gitee.com/quant1x/num@v0.3.2/math32/pow.go (about)

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