Golang Float.Mul方法代码示例

// sincos iterates a sin or cos Taylor series.
func sincos(name string, c Context, index int, x *big.Float, z *big.Float, exp uint64, factorial *big.Float) *big.Float {
	term := newFloat(c).Set(floatOne)
	for j := 0; j < index; j++ {
		term.Mul(term, x)
	}
	xN := newFloat(c).Set(term)
	x2 := newFloat(c).Mul(x, x)
	n := newFloat(c)

	for loop := newLoop(c.Config(), name, x, 4); ; {
		// Invariant: factorial holds -1ⁿ*exponent!.
		factorial.Neg(factorial)
		term.Quo(term, factorial)
		z.Add(z, term)

		if loop.done(z) {
			break
		}
		// Advance x**index (multiply by x²).
		term.Mul(xN, x2)
		xN.Set(term)
		// Advance factorial.
		factorial.Mul(factorial, n.SetUint64(exp+1))
		factorial.Mul(factorial, n.SetUint64(exp+2))
		exp += 2
	}
	return z
}

// 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
}

func renderFloat(img *image.RGBA) {
	var yminF, ymaxMinF, heightF big.Float
	yminF.SetInt64(ymin)
	ymaxMinF.SetInt64(ymax - ymin)
	heightF.SetInt64(height)

	var xminF, xmaxMinF, widthF big.Float
	xminF.SetInt64(xmin)
	xmaxMinF.SetInt64(xmax - xmin)
	widthF.SetInt64(width)

	var y, x big.Float
	for py := int64(0); py < height; py++ {
		// y := float64(py)/height*(ymax-ymin) + ymin
		y.SetInt64(py)
		y.Quo(&y, &heightF)
		y.Mul(&y, &ymaxMinF)
		y.Add(&y, &yminF)

		for px := int64(0); px < width; px++ {
			// x := float64(px)/width*(xmax-xmin) + xmin
			x.SetInt64(px)
			x.Quo(&x, &widthF)
			x.Mul(&x, &xmaxMinF)
			x.Add(&x, &xminF)

			c := mandelbrotFloat(&x, &y)
			if c == nil {
				c = color.Black
			}
			img.Set(int(px), int(py), c)
		}
	}
}

// sincos iterates a sin or cos Taylor series.
func sincos(name string, index int, x, z, exponent, factorial *big.Float) *big.Float {
	plus := false
	term := newF().Set(floatOne)
	for j := 0; j < index; j++ {
		term.Mul(term, x)
	}
	xN := newF().Set(term)
	x2 := newF().Mul(x, x)

	loop := newLoop(name, x, 4)
	for {
		// Invariant: factorial holds exponent!.
		term.Quo(term, factorial)
		if plus {
			z.Add(z, term)
		} else {
			z.Sub(z, term)
		}
		plus = !plus

		if loop.terminate(z) {
			break
		}
		// Advance x**index (multiply by x²).
		term.Mul(xN, x2)
		xN.Set(term)
		// Advance exponent and factorial.
		exponent.Add(exponent, floatOne)
		factorial.Mul(factorial, exponent)
		exponent.Add(exponent, floatOne)
		factorial.Mul(factorial, exponent)
	}
	return z
}

// 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)
	}
}

// Multiply
func (a Scalar) Mul_S(b S) S {
	var x, y big.Float
	x = big.Float(a)
	y = big.Float(b.(Scalar))
	z := x.Mul(&x, &y)
	return (Scalar)(*z)

}

func abs(x *complexFloat) *big.Float {
	r := new(big.Float).SetPrec(prec)
	i := new(big.Float).SetPrec(prec)
	r.Copy(x.r)
	i.Copy(x.i)

	r.Mul(r, x.r) // r^2
	i.Mul(i, x.i) // i^2

	r.Add(r, i) // r^2 + i^2

	return sqrtFloat(r) // sqrt(r^2 + i^2)
}

// 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)

}

// sqrt for big.Float
func sqrt(given *big.Float) *big.Float {
	const prec = 200
	steps := int(math.Log2(prec))
	given.SetPrec(prec)
	half := new(big.Float).SetPrec(prec).SetFloat64(0.5)
	x := new(big.Float).SetPrec(prec).SetInt64(1)
	t := new(big.Float)
	for i := 0; i <= steps; i++ {
		t.Quo(given, x)
		t.Add(x, t)
		t.Mul(half, t)
	}
	return x
}

func mandelbrotFloat(a, b *big.Float) color.Color {
	var x, y, nx, ny, x2, y2, f2, f4, r2, tmp big.Float
	f2.SetInt64(2)
	f4.SetInt64(4)
	x.SetInt64(0)
	y.SetInt64(0)

	defer func() { recover() }()

	for n := uint8(0); n < iterations; n++ {
		// Not update x2 and y2
		// because they are already updated in the previous loop
		nx.Sub(&x2, &y2)
		nx.Add(&nx, a)

		tmp.Mul(&x, &y)
		ny.Mul(&f2, &tmp)
		ny.Add(&ny, b)

		x.Set(&nx)
		y.Set(&ny)

		x2.Mul(&x, &x)
		y2.Mul(&y, &y)
		r2.Add(&x2, &y2)

		if r2.Cmp(&f4) > 0 {
			return color.Gray{255 - contrast*n}
		}
	}
	return color.Black
}

// 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())
}

// hypot for big.Float
func hypot(p, q *big.Float) *big.Float {
	// special cases
	switch {
	case p.IsInf() || q.IsInf():
		return big.NewFloat(math.Inf(1))
	}
	p = p.Abs(p)
	q = q.Abs(q)
	if p.Cmp(p) < 0 {
		p, q = q, p
	}
	if p.Cmp(big.NewFloat(0)) == 0 {
		return big.NewFloat(0)
	}
	q = q.Quo(q, p)
	return sqrt(q.Mul(q, q).Add(q, big.NewFloat(1))).Mul(q, p)
}

func sqrtFloat(x *big.Float) *big.Float {
	t1 := new(big.Float).SetPrec(prec)
	t2 := new(big.Float).SetPrec(prec)
	t1.Copy(x)

	// Iterate.
	// x{n} = (x{n-1}+x{0}/x{n-1}) / 2
	for i := 0; i <= steps; i++ {
		if t1.Cmp(zero) == 0 || t1.IsInf() {
			return t1
		}
		t2.Quo(x, t1)
		t2.Add(t2, t1)
		t1.Mul(half, t2)
	}

	return t1
}

// 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
}

// This example shows how to use big.Float to compute the square root of 2 with
// a precision of 200 bits, and how to print the result as a decimal number.
func Example_sqrt2() {
	// We'll do computations with 200 bits of precision in the mantissa.
	const prec = 200

	// Compute the square root of 2 using Newton's Method. We start with
	// an initial estimate for sqrt(2), and then iterate:
	//     x_{n+1} = 1/2 * ( x_n + (2.0 / x_n) )

	// Since Newton's Method doubles the number of correct digits at each
	// iteration, we need at least log_2(prec) steps.
	steps := int(math.Log2(prec))

	// Initialize values we need for the computation.
	two := new(big.Float).SetPrec(prec).SetInt64(2)
	half := new(big.Float).SetPrec(prec).SetFloat64(0.5)

	// Use 1 as the initial estimate.
	x := new(big.Float).SetPrec(prec).SetInt64(1)

	// 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.
	t := new(big.Float)

	// Iterate.
	for i := 0; i <= steps; i++ {
		t.Quo(two, x)  // t = 2.0 / x_n
		t.Add(x, t)    // t = x_n + (2.0 / x_n)
		x.Mul(half, t) // x_{n+1} = 0.5 * t
	}

	// We can use the usual fmt.Printf verbs since big.Float implements fmt.Formatter
	fmt.Printf("sqrt(2) = %.50f\n", x)

	// Print the error between 2 and x*x.
	t.Mul(x, x) // t = x*x
	fmt.Printf("error = %e\n", t.Sub(two, t))

	// Output:
	// sqrt(2) = 1.41421356237309504880168872420969807856967187537695
	// error = 0.000000e+00
}