Golang math.Cbrt函数代码示例

func test_cbrt(x float64, t *testing.T) {
	if (math.IsNaN(math.Cbrt(x)) && !math.IsNaN(Cbrt(x))) ||
		(math.IsNaN(Cbrt(x)) && !math.IsNaN(math.Cbrt(x))) ||
		math.Abs(Cbrt(x)-math.Cbrt(x)) > math.Abs(math.Cbrt(x))/1e-15 {
		t.Errorf("Cbrt(%v) is %v, not %v\n", x, Cbrt(x), math.Cbrt(x))
	}
}

func handleConn(conn net.Conn) {
	defer conn.Close()
	for {
		conn.SetReadDeadline(time.Now().Add(10 * time.Second))
		strReq, err := read(conn)
		if err != nil {
			if err == io.EOF {
				fmt.Println("The connection is closed by another side. (Server)")
			} else {
				fmt.Printf("Read Error: %s (Server)\n", err)
			}
			break
		}
		fmt.Printf("Received request: %s (Server)\n", strReq)

		i32Req, err := strconv.Atoi(strReq)
		if err != nil {
			n, err := write(conn, err.Error())
			if err != nil {
				fmt.Printf("Write Eoor (writen %d bytes:) %s (Server)\n", err)
			}
			fmt.Printf("Sent response (written %d bytes %s (server)\n", n, err)
			continue
		}

		f64Resp := math.Cbrt(float64(i32Req))
		respMsg := fmt.Sprintf("The cube root of %d is %f.", i32Req, f64Resp)
		n, err := write(conn, respMsg)
		if err != nil {
			fmt.Printf("Write Error: %s (Server)\n", err)
		}
		fmt.Printf("Sent response (writtn %d bytes:) %s (Server)\n", n, respMsg)
	}
}

// AnomalyDistance returns true anomaly and distance for near-parabolic orbits.
//
// Distance r returned in AU.
// An error is returned if the algorithm fails to converge.
func (e *Elements) AnomalyDistance(jde float64) (ν unit.Angle, r float64, err error) {
	// fairly literal translation of code on p. 246
	q1 := base.K * math.Sqrt((1+e.Ecc)/e.PDis) / (2 * e.PDis) // line 20
	g := (1 - e.Ecc) / (1 + e.Ecc)                            // line 20

	t := jde - e.TimeP // line 22
	if t == 0 {        // line 24
		return 0, e.PDis, nil
	}
	d1, d := 10000., 1e-9        // line 14
	q2 := q1 * t                 // line 28
	s := 2. / (3 * math.Abs(q2)) // line 30
	s = 2 / math.Tan(2*math.Atan(math.Cbrt(math.Tan(math.Atan(s)/2))))
	if t < 0 { // line 34
		s = -s
	}
	if e.Ecc != 1 { // line 36
		l := 0 // line 38
		for {
			s0 := s // line 40
			z := 1.
			y := s * s
			g1 := -y * s
			q3 := q2 + 2*g*s*y/3 // line 42
			for {
				z += 1                          // line 44
				g1 = -g1 * g * y                // line 46
				z1 := (z - (z+1)*g) / (2*z + 1) // line 48
				f := z1 * g1                    // line 50
				q3 += f                         // line 52
				if z > 50 || math.Abs(f) > d1 { // line 54
					return 0, 0, errors.New("No convergence")
				}
				if math.Abs(f) <= d { // line 56
					break
				}
			}
			l++ // line 58
			if l > 50 {
				return 0, 0, errors.New("No convergence")
			}
			for {
				s1 := s // line 60
				s = (2*s*s*s/3 + q3) / (s*s + 1)
				if math.Abs(s-s1) <= d { // line 62
					break
				}
			}
			if math.Abs(s-s0) <= d { // line 64
				break
			}
		}
	}
	ν = unit.Angle(2 * math.Atan(s))               // line 66
	r = e.PDis * (1 + e.Ecc) / (1 + e.Ecc*ν.Cos()) // line 68
	if ν < 0 {                                     // line 70
		ν += 2 * math.Pi
	}
	return
}

func BenchmarkCbrtFloat64ToUint64(b *testing.B) {
	k := uint(32)
	for i := 0; i < b.N; i++ {
		for j := uint64(0x10000); j < 0x4000000; j += 0x10000 {
			_ = uint64(math.Cbrt(float64(j << k)))
		}
	}
}

// AnomalyDistance returns true anomaly and distance of a body in a parabolic orbit of the Sun.
//
// True anomaly ν returned in radians.
// Distance r returned in AU.
func (e *Elements) AnomalyDistance(jde float64) (ν, r float64) {
	W := 3 * base.K / math.Sqrt2 * (jde - e.TimeP) / e.PDis / math.Sqrt(e.PDis)
	G := W * .5
	Y := math.Cbrt(G + math.Sqrt(G*G+1))
	s := Y - 1/Y
	ν = 2 * math.Atan(s)
	r = e.PDis * (1 + s*s)
	return
}

func (*lstarCompander) Compand(p Point) Point {
	l := p[0]
	if l < 0.0 {
		l = -l
	}
	if l <= 216.0/24389.0 {
		l = l * 24389.0 / 2700.0
	} else {
		l = 1.16*math.Cbrt(l) - 0.16
	}
	if p[0] < 0.0 {
		l = -l
	}
	p[0] = l

	l = p[1]
	if l < 0.0 {
		l = -l
	}
	if l <= 216.0/24389.0 {
		l = l * 24389.0 / 2700.0
	} else {
		l = 1.16*math.Cbrt(l) - 0.16
	}
	if p[1] < 0.0 {
		l = -l
	}
	p[1] = l

	l = p[2]
	if l < 0.0 {
		l = -l
	}
	if l <= 216.0/24389.0 {
		l = l * 24389.0 / 2700.0
	} else {
		l = 1.16*math.Cbrt(l) - 0.16
	}
	if p[2] < 0.0 {
		l = -l
	}
	p[2] = l
	return p
}

// TestIRoot16BitCbrt tests IRoot() with all 16-bit numbers and 3.
func TestIRoot16BitCbrt(t *testing.T) {
	for n := 0; n <= math.MaxInt16; n++ {
		expectedR := int64(math.Cbrt(float64(n)))
		r := IRoot(big.NewInt(int64(n)), 3)
		if r.Cmp(big.NewInt(expectedR)) != 0 {
			t.Errorf("For n=%d and p=3, expected r=%d, got r=%s",
				&n, expectedR, r)
		}
	}
}

func BenchmarkMatmult(b *testing.B) {
	b.StopTimer()
	// matmult is O(N**3) but testing expects O(b.N),
	// so we need to take cube root of b.N
	n := int(math.Cbrt(float64(b.N))) + 1
	A := makeMatrix(n)
	B := makeMatrix(n)
	C := makeMatrix(n)
	b.StartTimer()
	matmult(nil, A, B, C, 0, n, 0, n, 0, n, 8)
}

func XyzToLuvWhiteRef(x, y, z float64, wref [3]float64) (l, u, v float64) {
	if y/wref[1] <= 6.0/29.0*6.0/29.0*6.0/29.0 {
		l = y / wref[1] * 29.0 / 3.0 * 29.0 / 3.0 * 29.0 / 3.0
	} else {
		l = 1.16*math.Cbrt(y/wref[1]) - 0.16
	}
	ubis, vbis := xyz_to_uv(x, y, z)
	un, vn := xyz_to_uv(wref[0], wref[1], wref[2])
	u = 13.0 * l * (ubis - un)
	v = 13.0 * l * (vbis - vn)
	return
}

// Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
// This solution is adapted from Geocentric::Reverse.
func astroid(x, y float64) float64 {
	p, q := x*x, y*y
	r := (p + q - 1) / 6
	if q == 0 && r <= 0 {
		return 0
	}

	// Avoid possible division by zero when r = 0 by multiplying equations
	// for s and t by r^3 and r, resp.
	S := p * q / 4 // S = r^3 * s
	r2 := r * r
	r3 := r * r2
	// The discrimant of the quadratic equation for T3.  This is zero on
	// the evolute curve p^(1/3)+q^(1/3) = 1
	disc := S * (S + 2*r3)
	u := r
	if disc >= 0 {
		T3 := S + r3
		// Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
		// of precision due to cancellation.  The result is unchanged because
		// of the way the T is used in definition of u.
		if T3 < 0 {
			T3 += -math.Sqrt(disc)
		} else {
			T3 += math.Sqrt(disc) // T3 = (r * t)^3
		}
		// N.B. cbrt always returns the real root.  cbrt(-8) = -2.
		T := math.Cbrt(T3) // T = r * t
		// T can be zero; but then r2 / T -> 0.
		if T != 0 {
			u += T + r2/T
		}
	} else {
		// T is complex, but the way u is defined the result is real.
		ang := math.Atan2(math.Sqrt(-disc), -(S + r3))
		// There are three possible cube roots.  We choose the root which
		// avoids cancellation.  Note that disc < 0 implies that r < 0.
		u += 2 * r * math.Cos(ang/3)
	}

	v := math.Sqrt(u*u + q) // guaranteed positive
	// Avoid loss of accuracy when u < 0.
	var uv float64
	if u < 0 {
		uv = q / (v - u)
	} else {
		uv = u + v // u+v, guaranteed positive
	}
	w := (uv - q) / (2 * v) // positive?
	// Rearrange expression for k to avoid loss of accuracy due to
	// subtraction.  Division by 0 not possible because uv > 0, w >= 0.
	return uv / (math.Sqrt(uv+w*w) + w) // guaranteed positive
}

func BenchmarkCbrtFloat64ToUint(b *testing.B) {
	k := uint(1)
	if k<<32 != 0 {
		k = 32
	} else {
		k = 0
	}
	for i := 0; i < b.N; i++ {
		for j := uint(0x10000); j < 0x4000000; j += 0x10000 {
			_ = uint(math.Cbrt(float64(j << k)))
		}
	}
}

// Invert converts an XYZ point to LAB
func (t *LabTransformer) Invert(p XYZ) Lab {
	xr, yr, zr := p.X()/t.refWp.X(), p.Y()/t.refWp.Y(), p.Z()/t.refWp.Z()

	var fx, fy, fz float64
	if xr > CIEEps {
		fx = math.Cbrt(xr)
	} else {
		fx = (CIEKappa*xr + 16.0) / 116.0
	}
	if yr > CIEEps {
		fy = math.Cbrt(yr)
	} else {
		fy = (CIEKappa*yr + 16.0) / 116.0
	}
	if zr > CIEEps {
		fz = math.Cbrt(zr)
	} else {
		fz = (CIEKappa*zr + 16.0) / 116.0
	}

	return Lab{116*fy - 16.0, 500.0 * (fx - fy), 200.0 * (fy - fz)}
}

func Root(param float64, input *oproto.ValueStream) *oproto.ValueStream {
	output := &oproto.ValueStream{}
	for _, v := range input.Value {
		newv := &*v
		if param == 2 {
			newv.DoubleValue = math.Sqrt(v.DoubleValue)
		} else if param == 3 {
			newv.DoubleValue = math.Cbrt(v.DoubleValue)
		} else {
			newv.DoubleValue = math.Pow(v.DoubleValue, 1.0/param)
		}
		output.Value = append(output.Value, newv)
	}
	return output
}

func (t *LuvTransformer) Invert(p XYZ) Luv {
	d := p.X() + 15.0*p.Y() + 3.0*p.Z()

	var l, up, vp float64

	if d > 0 {
		up = 4.0 * p.X() / d
		vp = 9.0 * p.Y() / d
	}
	yr := p.Y() / t.refWp.Y()

	if yr > CIEEps {
		l = 116.0*math.Cbrt(yr) - 16.0
	} else {
		l = CIEKappa * yr
	}
	return Luv{l, 13.0 * l * (up - t.u0), 13.0 * l * (vp - t.v0)}
}

func angleContribution(path *geo.Path, index int, scale float64) *geo.Point {
	angle := geo.NewPoint(0, 0)

	if scale != 0.0 {
		n1 := path.GetAt(index - 1).Clone().Subtract(path.GetAt(index))
		n2 := path.GetAt(index + 1).Clone().Subtract(path.GetAt(index))

		len1 := n1.DistanceFrom(geo.NewPoint(0, 0))
		len2 := n2.DistanceFrom(geo.NewPoint(0, 0))

		n1.Normalize()
		n2.Normalize()

		// cbrt
		factor := math.Cbrt(n1.Dot(n2)) + 1
		angle = n1.Add(n2).Normalize().Scale(math.Min(len1, len2) * scale * factor)
	}

	return angle
}