gistfile1.go

package main

import (
	"fmt"
	"math"
	"math/big"
	"time"

	"github.com/ALTree/bigfloat"
	"github.com/ericlagergren/decimal"
	dmath "github.com/ericlagergren/decimal/math"
)

func main() {
	z := new(big.Float).SetPrec(64).SetUint64(2)
	t := time.Now()
	Exp(z, z)
	fd := time.Now().Sub(t)

	xx := decimal.WithPrecision(17).SetUint64(2)
	t = time.Now()
	dmath.Exp(xx, xx)
	dd := time.Now().Sub(t)

	zz := new(big.Float).SetPrec(64).SetUint64(2)
	t = time.Now()
	zz = bigfloat.Exp(zz)
	bfd := time.Now().Sub(t)

	fmt.Printf("big.Float CF   - %g: %s\n", z, fd)
	fmt.Printf("decimal.Big CF - %s: %s\n", xx, dd)
	fmt.Printf("bigfloat       - %g: %s\n", zz, bfd)
	fmt.Println("stdlib        - ", math.Exp(2))
}

func Exp(z, x *big.Float) *big.Float {
	if x.IsInf() {
		if x.Signbit() {
			// e ** -Inf = 0
			return z.SetUint64(0)
		}
		// e ** +Inf = +Inf
		return z.SetInf(false)
	}
	if x.Sign() == 0 {
		// e ** 0 = 1
		return z.SetUint64(1)
	}

	prec := z.Prec()
	r := new(big.Float)
	k := x.MantExp(r)
	g := expg{
		z:   r,
		pow: new(big.Float).SetPrec(prec).Mul(r, r),
		A:   new(big.Float).SetPrec(prec),
		B:   new(big.Float).SetPrec(prec),
	}

	wallis(z, &g)
	pow(z, new(big.Float).Copy(z), 1<<uint(k))
	return z
}

func pow(z, x *big.Float, y int) *big.Float {
	for i := 0; i < y-1; i++ {
		z.Mul(z, x)
	}
	return z
}

func wallis(z *big.Float, g generator) *big.Float {
	if !g.Next() {
		return z
	}

	prec := z.Prec()
	a := new(big.Float).SetPrec(prec)
	a_1 := new(big.Float).SetPrec(prec)
	b := new(big.Float).SetPrec(prec)
	b_1 := new(big.Float).SetPrec(prec)
	p := new(big.Float).SetPrec(prec)
	eps := new(big.Float).SetPrec(prec).Quo(one, new(big.Float).SetUint64(10000000000000))

	A, B := g.Term()

	a_1.SetUint64(1)
	a.Copy(B)
	b.SetUint64(1)

	for i := 0; g.Next() && !p.IsInf(); i++ {
		A, B = g.Term()

		z.Copy(a)
		a_1.Mul(a_1, A)
		a.Mul(a, B)
		a.Add(a, a_1)
		a_1.Copy(z)

		z.Copy(b)
		b_1.Mul(b_1, A)
		b.Mul(b, B)
		b.Add(b, b_1)
		b_1.Copy(z)

		z.Quo(a, b)
		if p.Sub(z, p).Abs(p).Cmp(eps) <= 0 {
			break
		}
		p.Copy(z)
	}
	return z
}

type generator interface {
	Term() (A, B *big.Float)
	Next() bool
}

type expg struct {
	z    *big.Float
	pow  *big.Float // z * z
	m    uint64     // term number
	A, B *big.Float // terms
}

func (e *expg) Next() bool { return true }

var (
	negfour = new(big.Float).SetInt64(-4)
	one     = new(big.Float).SetUint64(1)
	two     = new(big.Float).SetUint64(2)
	six     = new(big.Float).SetUint64(6)
	sixteen = new(big.Float).SetUint64(16)
)

func (e *expg) Term() (A, B *big.Float) {
	switch e.m {
	// [0, 1]
	case 0:
		e.A.SetUint64(0)
		e.B.SetUint64(1)
	// [2z, 2-z]
	case 1:
		e.A.Mul(two, e.z)
		e.B.Sub(two, e.z)
	// [z^2/6, 1]
	case 2:
		e.A.Quo(e.pow, six)
		e.B.SetUint64(1)
	// [(1/(16((m-1)^2)-4))(z^2), 1]
	default:
		// maxM is the largest m value we can use to compute 4(2m - 3)(2m - 1)
		// using unsigned integers.
		const maxM = 1518500252

		// 4(2m - 3)(2m - 1) ≡ 16(m - 1)^2 - 4
		if e.m <= maxM {
			e.A.SetUint64(16*((e.m-1)*(e.m-1)) - 4)
		} else {
			e.A.SetUint64(e.m - 1)

			// (m-1)^2
			e.A.Mul(e.A, e.A)

			// 16 * (m-1)^2 - 4 = 16 * (m-1)^2 + (-4)
			// e.A.FMA(sixteen, e.t.A, negfour)
			e.A.Mul(sixteen, e.A)
			e.A.Add(e.A, negfour)
		}

		// 1 / (16 * (m-1)^2 - 4)
		e.A.Quo(one, e.A)

		// (1 / (16 * (m-1)^2 - 4)) * (z^2)
		e.A.Mul(e.A, e.pow)

		// e.t.B is set to 1 inside case 2.
	}
	e.m++
	return e.A, e.B
}