github.com/c9s/go@v0.0.0-20180120015821-984e81f64e0c/src/math/pow.go (about)

     1  // 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  
    40  func pow(x, y float64) float64 {
    41  	switch {
    42  	case y == 0 || x == 1:
    43  		return 1
    44  	case y == 1:
    45  		return x
    46  	case IsNaN(x) || IsNaN(y):
    47  		return NaN()
    48  	case x == 0:
    49  		switch {
    50  		case y < 0:
    51  			if isOddInt(y) {
    52  				return Copysign(Inf(1), x)
    53  			}
    54  			return Inf(1)
    55  		case y > 0:
    56  			if isOddInt(y) {
    57  				return x
    58  			}
    59  			return 0
    60  		}
    61  	case IsInf(y, 0):
    62  		switch {
    63  		case x == -1:
    64  			return 1
    65  		case (Abs(x) < 1) == IsInf(y, 1):
    66  			return 0
    67  		default:
    68  			return Inf(1)
    69  		}
    70  	case IsInf(x, 0):
    71  		if IsInf(x, -1) {
    72  			return Pow(1/x, -y) // Pow(-0, -y)
    73  		}
    74  		switch {
    75  		case y < 0:
    76  			return 0
    77  		case y > 0:
    78  			return Inf(1)
    79  		}
    80  	case y == 0.5:
    81  		return Sqrt(x)
    82  	case y == -0.5:
    83  		return 1 / Sqrt(x)
    84  	}
    85  
    86  	absy := y
    87  	flip := false
    88  	if absy < 0 {
    89  		absy = -absy
    90  		flip = true
    91  	}
    92  	yi, yf := Modf(absy)
    93  	if yf != 0 && x < 0 {
    94  		return NaN()
    95  	}
    96  	if yi >= 1<<63 {
    97  		// yi is a large even int that will lead to overflow (or underflow to 0)
    98  		// for all x except -1 (x == 1 was handled earlier)
    99  		switch {
   100  		case x == -1:
   101  			return 1
   102  		case (Abs(x) < 1) == (y > 0):
   103  			return 0
   104  		default:
   105  			return Inf(1)
   106  		}
   107  	}
   108  
   109  	// ans = a1 * 2**ae (= 1 for now).
   110  	a1 := 1.0
   111  	ae := 0
   112  
   113  	// ans *= x**yf
   114  	if yf != 0 {
   115  		if yf > 0.5 {
   116  			yf--
   117  			yi++
   118  		}
   119  		a1 = Exp(yf * Log(x))
   120  	}
   121  
   122  	// ans *= x**yi
   123  	// by multiplying in successive squarings
   124  	// of x according to bits of yi.
   125  	// accumulate powers of two into exp.
   126  	x1, xe := Frexp(x)
   127  	for i := int64(yi); i != 0; i >>= 1 {
   128  		if xe < -1<<12 || 1<<12 < xe {
   129  			// catch xe before it overflows the left shift below
   130  			// Since i !=0 it has at least one bit still set, so ae will accumulate xe
   131  			// on at least one more iteration, ae += xe is a lower bound on ae
   132  			// the lower bound on ae exceeds the size of a float64 exp
   133  			// so the final call to Ldexp will produce under/overflow (0/Inf)
   134  			ae += xe
   135  			break
   136  		}
   137  		if i&1 == 1 {
   138  			a1 *= x1
   139  			ae += xe
   140  		}
   141  		x1 *= x1
   142  		xe <<= 1
   143  		if x1 < .5 {
   144  			x1 += x1
   145  			xe--
   146  		}
   147  	}
   148  
   149  	// ans = a1*2**ae
   150  	// if flip { ans = 1 / ans }
   151  	// but in the opposite order
   152  	if flip {
   153  		a1 = 1 / a1
   154  		ae = -ae
   155  	}
   156  	return Ldexp(a1, ae)
   157  }