Hi all! Just noticed that fractional powers of negative numbers seem to return exclusively nan. Is this a feature, a bug or am I doing something horribly wrong?
-1.pow(1/2) // nan, fair enough
-1.pow(1/3) // also nan !?!
I’d say it looks like a bug to me, especially since the documentation does not mention any explicit assumptions about the domain of pow.
If you try it with a complex number, you get back one of the non-real solutions:
Complex(-1,0).pow(1/3)
-> Complex( 0.5, 0.86602540378444 )
so that doesn’t really work as a satisfying workaround either.
I believe we are dependent on the C math library here. Do we have an obligation to do math slower, but better, than C?
hjh
No obligation to do anything as far as I’m concerned, but documenting the assumptions/limitations wouldn’t hurt.
Agreed.
In the language, I trace it to PyrMathOps.cpp line 1414 (other variants for int elsewhere):
case opPow:
SetRaw(a, pow(slotRawFloat(a), (double)slotRawInt(b)));
Assuming pow is from C standard libs, that would be the source of the limitation.
hjh
DISCLAIMER: I am not a mathematician.
Which answer would you expect?
With real numbers, only certain exponents seem to be possible:
You can raise a negative number to some fractional powers and get a real number answer, but only if the denominator of the fraction (in its lowest terms) is odd.
That’s the theory. In practice, such an exponent can’t be expressed accurately with floating point numbers, so it’s not really possible - unless the language has first class support for rational numbers.
That’s an excellent point. In base-two floating-point, the denominator must be even.
Consider decimals. If you have a finite number of digits after the decimal point, then the denominator may consist of any combination of prime factors of the base. For decimal, those prime factors are 2 and 5: 0.333 = 333/1000 = 333 / (2*2*2*5*5*5) but there is no finite-precision decimal where the denominator includes a factor of 3 or 7 or 11.
In binary, the base has only one prime factor: 2. And of course floating point can’t represent infinite precision (double-precision allows 53 bits in the mantissa).
a = 1/3;
a.high32Bits.asBinaryString(32)
a.low32Bits.asBinaryString(32)
-> 00111111110101010101010101010101
-> 01010101010101010101010101010101
=
sign bit = 0
exponent = 01111111101
mantissa = implicit 1. then 0101010101010101010101010101010101010101010101010101
So the denominator of this fraction 1.0101010101… will be a very large power of two, but it is a power of two, hence even, hence -1.pow(1/3) cannot be computed using a system of finite-precision binary fractions.
TL;DR -1.pow(0.33333333333333) = -1.pow(33333333333333 / 100000000000000) which is likewise an even denominator.
TL;DR It’s hard to argue with the IEEE; you’ll probably lose 
hjh
Hello,
Dynamic typing!
Scheme:
Chez Scheme Version 9.5.4
> (expt -1 (/ 1 2))
6.123233995736766e-17+1.0i
> (expt -1 (/ 1 3))
0.5000000000000001+0.8660254037844386i
(also python:)
Python 3.9.2 (default, Feb 28 2021, 17:03:44) [GCC 10.2.1]
>>> (-1) ** (1/2)
(6.123233995736766e-17+1j)
>>> (-1) ** (1/3)
(0.5000000000000001+0.8660254037844386j)
Squeak:
-1 raisedTo: 1/2. "-> ArithmeticError: Negative numbers don't have even roots. c.f. nthRoot:"
-1 raisedTo: 1/3. "-> -1"
Static typing…
GHCi, version 8.8.4
> (-1) ** (1/2)
NaN
> (-1) ** (1/3)
NaN
> import Data.Complex
> ((-1) :+ 0) ** (1/3)
0.5000000000000001 :+ 0.8660254037844386
> ((-1) :+ 0) ** (1/2)
6.123233995736766e-17 :+ 1.0
Best,
Rohan
Ps.
-1 raisedTo: 0.5. "-> ArithmeticError: Negative numbers can't be raised to float powers."
That’s interesting!
I find it funny that none of them yields 0.0 + 1j, which would be the precise correct answer. Probably some precision problem again.