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

Subversion Repositories openrisc

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

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 
func isOddInt(x float64) bool {
8 
        xi, xf := Modf(x)
9 
        return xf == 0 && int64(xi)&1 == 1
10 
}
11 
 
12 
// Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
13 
// updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
14 
 
15 
// Pow returns x**y, the base-x exponential of y.
16 
//
17 
// Special cases are (in order):
18 
//      Pow(x, ±0) = 1 for any x
19 
//      Pow(1, y) = 1 for any y
20 
//      Pow(x, 1) = x for any x
21 
//      Pow(NaN, y) = NaN
22 
//      Pow(x, NaN) = NaN
23 
//      Pow(±0, y) = ±Inf for y an odd integer < 0
24 
//      Pow(±0, -Inf) = +Inf
25 
//      Pow(±0, +Inf) = +0
26 
//      Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
27 
//      Pow(±0, y) = ±0 for y an odd integer > 0
28 
//      Pow(±0, y) = +0 for finite y > 0 and not an odd integer
29 
//      Pow(-1, ±Inf) = 1
30 
//      Pow(x, +Inf) = +Inf for |x| > 1
31 
//      Pow(x, -Inf) = +0 for |x| > 1
32 
//      Pow(x, +Inf) = +0 for |x| < 1
33 
//      Pow(x, -Inf) = +Inf for |x| < 1
34 
//      Pow(+Inf, y) = +Inf for y > 0
35 
//      Pow(+Inf, y) = +0 for y < 0
36 
//      Pow(-Inf, y) = Pow(-0, -y)
37 
//      Pow(x, y) = NaN for finite x < 0 and finite non-integer y
38 
func Pow(x, y float64) float64 {
39 
        switch {
40 
        case y == 0 || x == 1:
41 
                return 1
42 
        case y == 1:
43 
                return x
44 
        case y == 0.5:
45 
                return Sqrt(x)
46 
        case y == -0.5:
47 
                return 1 / Sqrt(x)
48 
        case IsNaN(x) || IsNaN(y):
49 
                return NaN()
50 
        case x == 0:
51 
                switch {
52 
                case y < 0:
53 
                        if isOddInt(y) {
54 
                                return Copysign(Inf(1), x)
55 
                        }
56 
                        return Inf(1)
57 
                case y > 0:
58 
                        if isOddInt(y) {
59 
                                return x
60 
                        }
61 
                        return 0
62 
                }
63 
        case IsInf(y, 0):
64 
                switch {
65 
                case x == -1:
66 
                        return 1
67 
                case (Abs(x) < 1) == IsInf(y, 1):
68 
                        return 0
69 
                default:
70 
                        return Inf(1)
71 
                }
72 
        case IsInf(x, 0):
73 
                if IsInf(x, -1) {
74 
                        return Pow(1/x, -y) // Pow(-0, -y)
75 
                }
76 
                switch {
77 
                case y < 0:
78 
                        return 0
79 
                case y > 0:
80 
                        return Inf(1)
81 
                }
82 
        }
83 
 
84 
        absy := y
85 
        flip := false
86 
        if absy < 0 {
87 
                absy = -absy
88 
                flip = true
89 
        }
90 
        yi, yf := Modf(absy)
91 
        if yf != 0 && x < 0 {
92 
                return NaN()
93 
        }
94 
        if yi >= 1<<63 {
95 
                return Exp(y * Log(x))
96 
        }
97 
 
98 
        // ans = a1 * 2**ae (= 1 for now).
99 
        a1 := 1.0
100 
        ae := 0
101 
 
102 
        // ans *= x**yf
103 
        if yf != 0 {
104 
                if yf > 0.5 {
105 
                        yf--
106 
                        yi++
107 
                }
108 
                a1 = Exp(yf * Log(x))
109 
        }
110 
 
111 
        // ans *= x**yi
112 
        // by multiplying in successive squarings
113 
        // of x according to bits of yi.
114 
        // accumulate powers of two into exp.
115 
        x1, xe := Frexp(x)
116 
        for i := int64(yi); i != 0; i >>= 1 {
117 
                if i&1 == 1 {
118 
                        a1 *= x1
119 
                        ae += xe
120 
                }
121 
                x1 *= x1
122 
                xe <<= 1
123 
                if x1 < .5 {
124 
                        x1 += x1
125 
                        xe--
126 
                }
127 
        }
128 
 
129 
        // ans = a1*2**ae
130 
        // if flip { ans = 1 / ans }
131 
        // but in the opposite order
132 
        if flip {
133 
                a1 = 1 / a1
134 
                ae = -ae
135 
        }
136 
        return Ldexp(a1, ae)
137 
}

powered by: WebSVN 2.1.0

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