Golang Float.Prec方法代码示例

// Compute the square root of n using Newton's Method. We start with
// an initial estimate for sqrt(n), and then iterate
//     x_{i+1} = 1/2 * ( x_i + (n / x_i) )
// Result is returned in x
func (e *Pslq) Sqrt(n, x *big.Float) {
	if n == x {
		panic("need distinct input and output")
	}
	if n.Sign() == 0 {
		x.Set(n)
		return
	} else if n.Sign() < 0 {
		panic("Sqrt of negative number")
	}
	prec := n.Prec()

	// Use the floating point square root as initial estimate
	nFloat64, _ := n.Float64()
	x.SetPrec(prec).SetFloat64(math.Sqrt(nFloat64))

	// We use t as a temporary variable. There's no need to set its precision
	// since big.Float values with unset (== 0) precision automatically assume
	// the largest precision of the arguments when used as the result (receiver)
	// of a big.Float operation.
	var t big.Float

	// Iterate.
	for {
		t.Quo(n, x)        // t = n / x_i
		t.Add(x, &t)       // t = x_i + (n / x_i)
		t.Mul(&e.half, &t) // x_{i+1} = 0.5 * t
		if x.Cmp(&t) == 0 {
			// Exit loop if no change to result
			break
		}
		x.Set(&t)
	}
}

// Sqrt returns the square root n.
func Sqrt(n *big.Float) *big.Float {
	prec := n.Prec()

	x := new(big.Float).SetPrec(prec).SetInt64(1)
	z := new(big.Float).SetPrec(prec).SetInt64(1)

	half := new(big.Float).SetPrec(prec).SetFloat64(0.5)
	t := new(big.Float).SetPrec(prec)

	for {
		z.Copy(x)

		t.Mul(x, x)
		t.Sub(t, n)
		t.Quo(t, x)
		t.Mul(t, half)
		x.Sub(x, t)

		if x.Cmp(z) == 0 {
			break
		}
	}

	return x
}

// NearestInt set res to the nearest integer to x
func (e *Pslq) NearestInt(x *big.Float, res *big.Int) {
	prec := x.Prec()
	var tmp big.Float
	tmp.SetPrec(prec)
	if x.Sign() >= 0 {
		tmp.Add(x, &e.half)
	} else {
		tmp.Sub(x, &e.half)
	}
	tmp.Int(res)
}

// Sqrt returns a big.Float representation of the square root of
// z. Precision is the same as the one of the argument. The function
// panics if z is negative, returns ±0 when z = ±0, and +Inf when z =
// +Inf.
func Sqrt(z *big.Float) *big.Float {

	// panic on negative z
	if z.Sign() == -1 {
		panic("Sqrt: argument is negative")
	}

	// √±0 = ±0
	if z.Sign() == 0 {
		return big.NewFloat(float64(z.Sign()))
	}

	// √+Inf  = +Inf
	if z.IsInf() {
		return big.NewFloat(math.Inf(+1))
	}

	// Compute √(a·2**b) as
	//   √(a)·2**b/2       if b is even
	//   √(2a)·2**b/2      if b > 0 is odd
	//   √(0.5a)·2**b/2    if b < 0 is odd
	//
	// The difference in the odd exponent case is due to the fact that
	// exp/2 is rounded in different directions when exp is negative.
	mant := new(big.Float)
	exp := z.MantExp(mant)
	switch exp % 2 {
	case 1:
		mant.Mul(big.NewFloat(2), mant)
	case -1:
		mant.Mul(big.NewFloat(0.5), mant)
	}

	// Solving x² - z = 0 directly requires a Quo call, but it's
	// faster for small precisions.
	//
	// Solving 1/x² - z = 0 avoids the Quo call and is much faster for
	// high precisions.
	//
	// Use sqrtDirect for prec <= 128 and sqrtInverse for prec > 128.
	var x *big.Float
	if z.Prec() <= 128 {
		x = sqrtDirect(mant)
	} else {
		x = sqrtInverse(mant)
	}

	// re-attach the exponent and return
	return x.SetMantExp(x, exp/2)

}

// compute √z using newton to solve
// t² - z = 0 for t
func sqrtDirect(z *big.Float) *big.Float {
	// f(t)/f'(t) = 0.5(t² - z)/t
	half := big.NewFloat(0.5)
	f := func(t *big.Float) *big.Float {
		x := new(big.Float).Mul(t, t) // x = t²
		x.Sub(x, z)                   // x = t² - z
		x.Mul(half, x)                // x = 0.5(t² - z)
		return x.Quo(x, t)            // return x = 0.5(t² - z)/t
	}

	// initial guess
	zf, _ := z.Float64()
	guess := big.NewFloat(math.Sqrt(zf))

	return newton(f, guess, z.Prec())
}

func withinMandLimit(z *BigComplex, limit *big.Float) bool {
	// Approximate cmplx.Abs
	negLimit := MakeBigFloat(0.0, limit.Prec())
	negLimit.Neg(limit)

	r := z.Real()
	i := z.Imag()

	rLimCmp := r.Cmp(limit)
	rNegLimCmp := r.Cmp(&negLimit)
	iLimCmp := i.Cmp(limit)
	iNegLimCmp := i.Cmp(&negLimit)

	within := rLimCmp == -1 && rNegLimCmp == 1
	within = within && iLimCmp == -1 && iNegLimCmp == 1
	return within
}

// fast path for z**w when w is an integer
func powInt(z *big.Float, w int) *big.Float {

	// get mantissa and exponent of z
	mant := new(big.Float)
	exp := z.MantExp(mant)

	// result's exponent
	exp = exp * w

	// result's mantissa
	x := big.NewFloat(1).SetPrec(z.Prec())

	// Classic right-to-left binary exponentiation
	for w > 0 {
		if w%2 == 1 {
			x.Mul(x, mant)
		}
		w >>= 1
		mant.Mul(mant, mant)
	}

	return new(big.Float).SetMantExp(x, exp)
}

// compute √z using newton to solve
// 1/t² - z = 0 for x and then inverting.
func sqrtInverse(z *big.Float) *big.Float {
	// f(t)/f'(t) = -0.5t(1 - zt²)
	nhalf := big.NewFloat(-0.5)
	one := big.NewFloat(1)
	f := func(t *big.Float) *big.Float {
		u := new(big.Float)
		u.Mul(t, t)                     // u = t²
		u.Mul(u, z)                     // u = zt²
		u.Sub(one, u)                   // u = 1 - zt²
		u.Mul(u, nhalf)                 // u = -0.5(1 - zt²)
		return new(big.Float).Mul(t, u) // x = -0.5t(1 - zt²)
	}

	// initial guess
	zf, _ := z.Float64()
	guess := big.NewFloat(1 / math.Sqrt(zf))

	// There's another operation after newton,
	// so we need to force it to return at least
	// a few guard digits. Use 32.
	x := newton(f, guess, z.Prec()+32)
	return x.Mul(z, x).SetPrec(z.Prec())
}

func ExampleFloat_Add() {
	// Operating on numbers of different precision.
	var x, y, z big.Float
	x.SetInt64(1000)          // x is automatically set to 64bit precision
	y.SetFloat64(2.718281828) // y is automatically set to 53bit precision
	z.SetPrec(32)
	z.Add(&x, &y)
	fmt.Printf("x = %s (%s, prec = %d, acc = %s)\n", &x, x.Format('p', 0), x.Prec(), x.Acc())
	fmt.Printf("y = %s (%s, prec = %d, acc = %s)\n", &y, y.Format('p', 0), y.Prec(), y.Acc())
	fmt.Printf("z = %s (%s, prec = %d, acc = %s)\n", &z, z.Format('p', 0), z.Prec(), z.Acc())
	// Output:
	// x = 1000 (0x.fap10, prec = 64, acc = Exact)
	// y = 2.718281828 (0x.adf85458248cd8p2, prec = 53, acc = Exact)
	// z = 1002.718282 (0x.faadf854p10, prec = 32, acc = Below)
}

// Log returns a big.Float representation of the natural logarithm of
// z. Precision is the same as the one of the argument. The function
// panics if z is negative, returns -Inf when z = 0, and +Inf when z =
// +Inf
func Log(z *big.Float) *big.Float {

	// panic on negative z
	if z.Sign() == -1 {
		panic("Log: argument is negative")
	}

	// Log(0) = -Inf
	if z.Sign() == 0 {
		return big.NewFloat(math.Inf(-1)).SetPrec(z.Prec())
	}

	prec := z.Prec() + 64 // guard digits

	one := big.NewFloat(1).SetPrec(prec)
	two := big.NewFloat(2).SetPrec(prec)
	four := big.NewFloat(4).SetPrec(prec)

	// Log(1) = 0
	if z.Cmp(one) == 0 {
		return big.NewFloat(0).SetPrec(z.Prec())
	}

	// Log(+Inf) = +Inf
	if z.IsInf() {
		return big.NewFloat(math.Inf(+1)).SetPrec(z.Prec())
	}

	x := new(big.Float).SetPrec(prec)

	// if 0 < z < 1 we compute log(z) as -log(1/z)
	var neg bool
	if z.Cmp(one) < 0 {
		x.Quo(one, z)
		neg = true
	} else {
		x.Set(z)
	}

	// We scale up x until x >= 2**(prec/2), and then we'll be allowed
	// to use the AGM formula for Log(x).
	//
	// Double x until the condition is met, and keep track of the
	// number of doubling we did (needed to scale back later).

	lim := new(big.Float)
	lim.SetMantExp(two, int(prec/2))

	k := 0
	for x.Cmp(lim) < 0 {
		x.Mul(x, x)
		k++
	}

	// Compute the natural log of x using the fact that
	//     log(x) = π / (2 * AGM(1, 4/x))
	// if
	//     x >= 2**(prec/2),
	// where prec is the desired precision (in bits)
	pi := pi(prec)
	agm := agm(one, x.Quo(four, x)) // agm = AGM(1, 4/x)

	x.Quo(pi, x.Mul(two, agm)) // reuse x, we don't need it

	if neg {
		x.Neg(x)
	}

	// scale the result back multiplying by 2**-k
	// reuse lim to reduce allocations.
	x.Mul(x, lim.SetMantExp(one, -k))

	return x.SetPrec(z.Prec())
}

func runTest(encoderDecoder EncoderDecoder, n *big.Float) (nBytes uint64) {
	y := newBigFloat(0)

	buf := new(bytes.Buffer)
	err := encoderDecoder.Encode(buf, n)
	if err != nil {
		panic(err)
	}
	nBytes += uint64(buf.Len())
	buf = bytes.NewBuffer(buf.Bytes())

	err = encoderDecoder.Decode(buf, y)
	nBytes += uint64(buf.Len())
	if n.Cmp(y) != 0 {
		panic(fmt.Sprintf("write and read are not the same: %v, %v - %d - %d, %d", stringOfBigFloat(n), stringOfBigFloat(y), nBytes, n.Prec(), y.Prec()))
	}
	return
}

// Pow returns a big.Float representation of z**w. Precision is the same as the one
// of the first argument. The function panics when z is negative.
func Pow(z *big.Float, w *big.Float) *big.Float {

	if z.Sign() < 0 {
		panic("Pow: negative base")
	}

	// Pow(z, 0) = 1.0
	if w.Sign() == 0 {
		return big.NewFloat(1).SetPrec(z.Prec())
	}

	// Pow(z, 1) = z
	// Pow(+Inf, n) = +Inf
	if w.Cmp(big.NewFloat(1)) == 0 || z.IsInf() {
		return new(big.Float).Copy(z)
	}

	// Pow(z, -w) = 1 / Pow(z, w)
	if w.Sign() < 0 {
		x := new(big.Float)
		zExt := new(big.Float).Copy(z).SetPrec(z.Prec() + 64)
		wNeg := new(big.Float).Neg(w)
		return x.Quo(big.NewFloat(1), Pow(zExt, wNeg)).SetPrec(z.Prec())
	}

	// w integer fast path
	if w.IsInt() {
		wi, _ := w.Int64()
		return powInt(z, int(wi))
	}

	// compute w**z as exp(z log(w))
	x := new(big.Float).SetPrec(z.Prec() + 64)
	logZ := Log(new(big.Float).Copy(z).SetPrec(z.Prec() + 64))
	x.Mul(w, logZ)
	x = Exp(x)
	return x.SetPrec(z.Prec())

}

// Exp returns a big.Float representation of exp(z). Precision is
// the same as the one of the argument. The function returns +Inf
// when z = +Inf, and 0 when z = -Inf.
func Exp(z *big.Float) *big.Float {

	// exp(0) == 1
	if z.Sign() == 0 {
		return big.NewFloat(1).SetPrec(z.Prec())
	}

	// Exp(+Inf) = +Inf
	if z.IsInf() && z.Sign() > 0 {
		return big.NewFloat(math.Inf(+1)).SetPrec(z.Prec())
	}

	// Exp(-Inf) = 0
	if z.IsInf() && z.Sign() < 0 {
		return big.NewFloat(0).SetPrec(z.Prec())
	}

	guess := new(big.Float)

	// try to get initial estimate using IEEE-754 math
	zf, _ := z.Float64()
	if zfs := math.Exp(zf); zfs == math.Inf(+1) || zfs == 0 {
		// too big or too small for IEEE-754 math,
		// perform argument reduction using
		//     e^{2z} = (e^z)²
		halfZ := new(big.Float).Mul(z, big.NewFloat(0.5))
		halfExp := Exp(halfZ.SetPrec(z.Prec() + 64))
		return new(big.Float).Mul(halfExp, halfExp).SetPrec(z.Prec())
	} else {
		// we got a nice IEEE-754 estimate
		guess.SetFloat64(zfs)
	}

	// f(t)/f'(t) = t*(log(t) - z)
	f := func(t *big.Float) *big.Float {
		x := new(big.Float)
		x.Sub(Log(t), z)
		return x.Mul(x, t)
	}

	x := newton(f, guess, z.Prec())

	return x
}