github.com/gnolang/gno@v0.0.0-20240520182011-228e9d0192ce/gnovm/stdlibs/math/pow.gno (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  	if Abs(x) >= (1 << 53) {
     9  		// 1 << 53 is the largest exact integer in the float64 format.
    10  		// Any number outside this range will be truncated before the decimal point and therefore will always be
    11  		// an even integer.
    12  		// Without this check and if x overflows int64 the int64(xi) conversion below may produce incorrect results
    13  		// on some architectures (and does so on arm64). See issue #57465.
    14  		return false
    15  	}
    16  
    17  	xi, xf := Modf(x)
    18  	return xf == 0 && int64(xi)&1 == 1
    19  }
    20  
    21  // Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
    22  // updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
    23  
    24  // Pow returns x**y, the base-x exponential of y.
    25  //
    26  // Special cases are (in order):
    27  //
    28  //	Pow(x, ±0) = 1 for any x
    29  //	Pow(1, y) = 1 for any y
    30  //	Pow(x, 1) = x for any x
    31  //	Pow(NaN, y) = NaN
    32  //	Pow(x, NaN) = NaN
    33  //	Pow(±0, y) = ±Inf for y an odd integer < 0
    34  //	Pow(±0, -Inf) = +Inf
    35  //	Pow(±0, +Inf) = +0
    36  //	Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
    37  //	Pow(±0, y) = ±0 for y an odd integer > 0
    38  //	Pow(±0, y) = +0 for finite y > 0 and not an odd integer
    39  //	Pow(-1, ±Inf) = 1
    40  //	Pow(x, +Inf) = +Inf for |x| > 1
    41  //	Pow(x, -Inf) = +0 for |x| > 1
    42  //	Pow(x, +Inf) = +0 for |x| < 1
    43  //	Pow(x, -Inf) = +Inf for |x| < 1
    44  //	Pow(+Inf, y) = +Inf for y > 0
    45  //	Pow(+Inf, y) = +0 for y < 0
    46  //	Pow(-Inf, y) = Pow(-0, -y)
    47  //	Pow(x, y) = NaN for finite x < 0 and finite non-integer y
    48  func Pow(x, y float64) float64 {
    49  	return pow(x, y)
    50  }
    51  
    52  func pow(x, y float64) float64 {
    53  	switch {
    54  	case y == 0 || x == 1:
    55  		return 1
    56  	case y == 1:
    57  		return x
    58  	case IsNaN(x) || IsNaN(y):
    59  		return NaN()
    60  	case x == 0:
    61  		switch {
    62  		case y < 0:
    63  			if Signbit(x) && isOddInt(y) {
    64  				return Inf(-1)
    65  			}
    66  			return Inf(1)
    67  		case y > 0:
    68  			if Signbit(x) && isOddInt(y) {
    69  				return x
    70  			}
    71  			return 0
    72  		}
    73  	case IsInf(y, 0):
    74  		switch {
    75  		case x == -1:
    76  			return 1
    77  		case (Abs(x) < 1) == IsInf(y, 1):
    78  			return 0
    79  		default:
    80  			return Inf(1)
    81  		}
    82  	case IsInf(x, 0):
    83  		if IsInf(x, -1) {
    84  			return Pow(1/x, -y) // Pow(-0, -y)
    85  		}
    86  		switch {
    87  		case y < 0:
    88  			return 0
    89  		case y > 0:
    90  			return Inf(1)
    91  		}
    92  	case y == 0.5:
    93  		return Sqrt(x)
    94  	case y == -0.5:
    95  		return 1 / Sqrt(x)
    96  	}
    97  
    98  	yi, yf := Modf(Abs(y))
    99  	if yf != 0 && x < 0 {
   100  		return NaN()
   101  	}
   102  	if yi >= 1<<63 {
   103  		// yi is a large even int that will lead to overflow (or underflow to 0)
   104  		// for all x except -1 (x == 1 was handled earlier)
   105  		switch {
   106  		case x == -1:
   107  			return 1
   108  		case (Abs(x) < 1) == (y > 0):
   109  			return 0
   110  		default:
   111  			return Inf(1)
   112  		}
   113  	}
   114  
   115  	// ans = a1 * 2**ae (= 1 for now).
   116  	a1 := 1.0
   117  	ae := 0
   118  
   119  	// ans *= x**yf
   120  	if yf != 0 {
   121  		if yf > 0.5 {
   122  			yf--
   123  			yi++
   124  		}
   125  		a1 = Exp(yf * Log(x))
   126  	}
   127  
   128  	// ans *= x**yi
   129  	// by multiplying in successive squarings
   130  	// of x according to bits of yi.
   131  	// accumulate powers of two into exp.
   132  	x1, xe := Frexp(x)
   133  	for i := int64(yi); i != 0; i >>= 1 {
   134  		if xe < -1<<12 || 1<<12 < xe {
   135  			// catch xe before it overflows the left shift below
   136  			// Since i !=0 it has at least one bit still set, so ae will accumulate xe
   137  			// on at least one more iteration, ae += xe is a lower bound on ae
   138  			// the lower bound on ae exceeds the size of a float64 exp
   139  			// so the final call to Ldexp will produce under/overflow (0/Inf)
   140  			ae += xe
   141  			break
   142  		}
   143  		if i&1 == 1 {
   144  			a1 *= x1
   145  			ae += xe
   146  		}
   147  		x1 *= x1
   148  		xe <<= 1
   149  		if x1 < .5 {
   150  			x1 += x1
   151  			xe--
   152  		}
   153  	}
   154  
   155  	// ans = a1*2**ae
   156  	// if y < 0 { ans = 1 / ans }
   157  	// but in the opposite order
   158  	if y < 0 {
   159  		a1 = 1 / a1
   160  		ae = -ae
   161  	}
   162  	return Ldexp(a1, ae)
   163  }