OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgo/] [go/] [math/] [cbrt.go] - Blame information for rev 867

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 747 jeremybenn 
// Copyright 2009 The Go Authors. All rights reserved.
2 
// Use of this source code is governed by a BSD-style
3 
// license that can be found in the LICENSE file.
4 
 
5 
package math
6 
 
7 
/*
8 
        The algorithm is based in part on "Optimal Partitioning of
9 
        Newton's Method for Calculating Roots", by Gunter Meinardus
10 
        and G. D. Taylor, Mathematics of Computation © 1980 American
11 
        Mathematical Society.
12 
        (http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010)
13 
*/
14 
 
15 
// Cbrt returns the cube root of its argument.
16 
//
17 
// Special cases are:
18 
//      Cbrt(±0) = ±0
19 
//      Cbrt(±Inf) = ±Inf
20 
//      Cbrt(NaN) = NaN
21 
func Cbrt(x float64) float64 {
22 
        const (
23 
                A1 = 1.662848358e-01
24 
                A2 = 1.096040958e+00
25 
                A3 = 4.105032829e-01
26 
                A4 = 5.649335816e-01
27 
                B1 = 2.639607233e-01
28 
                B2 = 8.699282849e-01
29 
                B3 = 1.629083358e-01
30 
                B4 = 2.824667908e-01
31 
                C1 = 4.190115298e-01
32 
                C2 = 6.904625373e-01
33 
                C3 = 6.46502159e-02
34 
                C4 = 1.412333954e-01
35 
        )
36 
        // special cases
37 
        switch {
38 
        case x == 0 || IsNaN(x) || IsInf(x, 0):
39 
                return x
40 
        }
41 
        sign := false
42 
        if x < 0 {
43 
                x = -x
44 
                sign = true
45 
        }
46 
        // Reduce argument and estimate cube root
47 
        f, e := Frexp(x) // 0.5 <= f < 1.0
48 
        m := e % 3
49 
        if m > 0 {
50 
                m -= 3
51 
                e -= m // e is multiple of 3
52 
        }
53 
        switch m {
54 
        case 0: // 0.5 <= f < 1.0
55 
                f = A1*f + A2 - A3/(A4+f)
56 
        case -1:
57 
                f *= 0.5 // 0.25 <= f < 0.5
58 
                f = B1*f + B2 - B3/(B4+f)
59 
        default: // m == -2
60 
                f *= 0.25 // 0.125 <= f < 0.25
61 
                f = C1*f + C2 - C3/(C4+f)
62 
        }
63 
        y := Ldexp(f, e/3) // e/3 = exponent of cube root
64 
 
65 
        // Iterate
66 
        s := y * y * y
67 
        t := s + x
68 
        y *= (t + x) / (s + t)
69 
        // Reiterate
70 
        s = (y*y*y - x) / x
71 
        y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
72 
        if sign {
73 
                y = -y
74 
        }
75 
        return y
76 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.