Athanasius has asked for the wisdom of the Perl Monks concerning the following question:
There are times when Perl’s inbuilt exponentiation operator, **, just doesn’t cut it:
16:46 >perl wE "my $cube_root = (8)**(1/3); say $cube_root;"
NaN
16:46 >perl v
This is perl 5, version 32, subversion 1 (v5.32.1) built for MSWin32x
+64multithread
I note in the documentation that ** “...is implemented using C's pow(3) function...”, and that’s where the limitation appears to lie. (I’m using Strawberry Perl running under Windows 8.1 64bit. I get the same results using Raku. I don’t know if C’s pow(3) function has the same limitations under Unix?)
I came up with the following, which works for my usecase but is really a kludge:
sub cube_root
{
my ($n) = @_;
return (abs( $n ) ** (1 / 3)) * ($n < 0 ? 1 : 1);
}
So, what’s the quickest/simplest/correct way to get Perl to act like the calculator app that comes with Windows, and give 2 as the cube root of 8? (I’m only interested in the real solution.) I had a quick look through CPAN, but nothing stood out. I’m probably overlooking the obvious.
Thanks,
Re: How to get better exponentiation?
by haukex (Archbishop) on Jan 22, 2022 at 08:57 UTC

I don’t know if C’s pow(3) function has the same limitations under Unix?
Yes, I get NaN on Linux too, and pow(3) says "The functions pow(x, y), powf(x, y), and powl(x, y) raise an invalid exception and return an NaN if x < 0 and y is not an integer." and similarly, on my system the manpage says "If x is a finite value less than 0, and y is a finite noninteger, a domain error occurs, and a NaN is returned." If you ask Wolfram Alpha, you have to explicitly tell it you want the realvalued root instead of the principal root, though I'm not enough of an expert to know whether that's relevant here (I did find out you can get all the principal roots via Math::Complex's root). Interestingly, Math::BigInt and Math::BigFloat's >broot(3) also return NaN, so throwing a use bignum; in the program unfortunately doesn't help in this case.
To get the cube root, you can apparently use cbrt(3), exposed via POSIX as of Perl v5.22:
$ perl wMstrict MPOSIX=cbrt le 'print cbrt(8)'
2
 [reply] [d/l] [select] 

 [reply] [d/l] [select] 

$ perl MMath::Prime::Util::GMP=rootreal le 'print rootreal(8, 3)'
2.000000000000000000000000000000000000000
Use 0+rootreal(...) to get just 2 in this case. powreal(8, 1/3) suffers a bit from the floatingpoint inaccuracies in 1/3. Note both functions support an optional third argument to specify the number of significant digits, e.g. powreal(8, 1/3, 5) is "2.0000".
Update: However, note the issues raised by vr and syphilis here!  [reply] [d/l] [select] 
Re: How to get better exponentiation?
by syphilis (Archbishop) on Jan 22, 2022 at 13:43 UTC

I came up with the following, which works for my usecase ...
Note that the floating point value 1/3 is not the same as the rational value 1/3, and while your subroutine works as you want for your use case, it won't work correctly for (eg) cube_root( 91197 ** 3 ), where it will return 91196.9999999999.
For general purposes, you therefore really need an implementation that performs cbrt(x), as opposed to pow(x, 1 / 3).
haukex has already shown that the excellent Math::Prime::Util::GMP module provides what you're after.
My own Math::MPFR module provides the same capability (courtesy of the mpfr C library):
C:\>perl MMath::MPFR=":mpfr" le "$op = Math::MPFR>new(9117 ** 3);
+Rmpfr_cbrt($op, $op, MPFR_RNDN); print $op;"
9.117e3
or, the more general:
C:\>perl MMath::MPFR=":mpfr" le "$op = Math::MPFR>new(9117 ** 3);
+Rmpfr_rootn_ui($op, $op, 3, MPFR_RNDN); print $op;"
9.117e3
Cheers, Rob
 [reply] [d/l] [select] 

> Note that the floating point value 1/3 is not the same as the rational value 1/3,
Good point, that might explain why it's considered undefined behavior by many.
And an approximation algorithm designed to work with floating point exponents won't be able to handle anything like 1/$odd differently.
Cheers Rolf
_{(addicted to the Perl Programming Language :)
Wikisyntax for the Monastery
}
 [reply] [d/l] 
Re: How to get better exponentiation?
by vr (Curate) on Jan 22, 2022 at 18:10 UTC

Calculator 10.2103.8.0
© 2021 Microsoft. All rights reserved
1 ^ ( 1 / 2 ) =
Invalid input
1 ^ ( 2 / 4 ) =
1
1 ^ ( 1 / 3 ) =
1
1 ^ ( 2 / 6 ) =
1
32 ^ ( 1 / 5 ) =
2
32 ^ ( 2 / 10 ) =
2
32 ^ 0.2 =
2
And so on. You see the pattern.
Math::Prime::Util::GMP
I'm not sure, maybe it's un(der)documented (about allowed inputs) or buggy/inconsistent w.r.t. functions mentioned:
perl MMath::Prime::Util::GMP=:all le 'print powreal(4, .5)'
2.000000000000000000000000000000000000000
perl MMath::Prime::Util::GMP=:all le "print rootreal(4, 2)"
# aborted, core dumped
I think if task is only to ever get real roots, then your "kludge" (and checking for integer parity and a sign of argument, then disallowing obvious combination) looks farfar less "kludgier" than, say, calling Math::Complex::root and then grepping list for (possibly) zero imaginary part, or something like that.
If, OTOH, what's required is my_real_pow with any (+/) real base and any (positive?) real exponent, then solution may be to use Math::BigRat::parts on exponent, and then disallowing "negative base and even exponent denominator", while numerator parity influences sign of result if base is negative. Just an idea to explore (what about Math::BigRat result of a reciprocal, will roundtrip call for our function be consistent?). Then real answer for (8)**(2/3) would be +4 same as Wolfram Alpha. Edit 23/01: Apparently, Math::BigRat just won't help with this approach (I should have checked what it returns for most simple 2/3 argument), if exponent already happens to be floating point approximation. CPAN doesn't seem to have solutions to convert such floats to reasonable rationals. For (nonPerl) point of reference, in J 8byte doubles, results of division such as 2%3 or (just randomly typed) 73364%294557 are converted to rationals with numerator/denominator exactly as shown.
 [reply] [d/l] [select] 

 [reply] [d/l] 
Re: How to get better exponentiation?
by salva (Canon) on Jan 22, 2022 at 16:22 UTC

Note that Perl converts 1/3 to a floating point number which is close but not exactly 1/3.
Then, 8 powered to that 1/3 aproximation is not a real number but a complex one (well, actually a very big set of roots).  [reply] 
Re: How to get better exponentiation? (undefined)
by LanX (Saint) on Jan 22, 2022 at 12:03 UTC

I had my doubts and did some quick research
German Wikipedia states that it's disputed, if the cube root of a negative real number was defined.°
I couldn't find it in English sources in the hurry, but I suppose it's the old problem of having functions which fulfill certain qualities like reversibility in order to play well with the whole mathematical model.
In the end you'll have to handle two definitions of odd roots.
Hence your approach to define your own sub cube_root() looks fine for me. :)
Cheers Rolf
_{(addicted to the Perl Programming Language :)
Wikisyntax for the Monastery
}
°) https://de.wikipedia.org/wiki/Wurzel_(Mathematik)#Wurzeln_aus_negativen_Zahlen
Bezüglich der ungeraden Wurzeln aus negativen Zahlen werden folgende Positionen vertreten:
 Wurzeln aus negativen Zahlen sind generell nicht definiert ...
 Wurzeln aus negativen Zahlen sind definiert ...
update
Here the google translate for the laziest of lazy ...
Regarding the odd roots of negative numbers, the following positions are held:
 Roots of negative numbers are generally undefined...
 Roots of negative numbers are defined...
 [reply] [d/l] 

German Wikipedia states that it's disputed, if the cube root of a negative real number was defined
I finally got around to taking a look at the C89, C99 and C11 standards.
There's no mention of "cbrt" in C89, but in C99 and C11 it's quite acceptable for cbrt/cbrtf/cbrtl to take a negative argument.
However, if pow/powf/powl are given a first argument that is negative && the second argument has a noninteger value, then a "domain error" occurs.
Cheers, Rob
 [reply] 

Yes, thanks.
But I was talking about mathematical definitions and these are CS standards.
For instance: pure math has no big notion of floating point numbers.
Personally I'm fine with allowing root($x,$o) with $x<0 and $o odd integer in a computer.
But I could imagine reasons in the realm of mathematical modeling of functions to consider them undefined.
> the second argument has a noninteger value
I'd say because there is no way to express 1/$o loss free as binary floating point number ('$o odd integer')
update
see also https://en.wikipedia.org/wiki/Cube_root#Complex_numbers
With this definition, the principal cube root of a negative number is a complex number, and for instance 3√−8 will not be −2, but rather 1 + i√3.
Cheers Rolf
_{(addicted to the Perl Programming Language :)
Wikisyntax for the Monastery
}
 [reply] [d/l] [select] 

Re: How to get better exponentiation?
by syphilis (Archbishop) on Jan 23, 2022 at 01:15 UTC

So, what’s the quickest/simplest/correct way to get Perl to act like the calculator app that comes with Windows, and give 2 as the cube root of 8?
Update: The below is based upon my observations of the Windows calculator on Windows 7.
Some further reading indicates that things might be different on Windows 10.
I think that the calculator app, when calculating x ** y for x < 0 proceeds roughly as follows:
if (x < 0) {
return x ** y if y is an integer;
calculate n = 1 / y;
return nth root of x if n is an odd integer;
die "Invalid input",
}
else {
return x ** y;
}
As you can see, in addition to needing an exponentiation function, that method also requires a function that calculates the nth root.
If you use that method and make use of POSIX::cbrt(), you should see that it DWYMs for the case 8 ** (1 / 3)
It might also be that the windows calculator is performing decimal arithmetic  and that could also produce results that differ slightly from perl in some cases.
Does anyone know for sure the number base that's used by the windows calculator ?
Cheers, Rob
 [reply] [d/l] [select] 
Re: How to get better exponentiation?
by rizzo (Curate) on Jan 23, 2022 at 01:20 UTC

So, what’s the quickest/simplest/correct way to get Perl to act like the calculator app that comes with Windows ...
With the emphasis on "calculator app" Gnu Octave could be an option. There is the Inline::Octave module that lets you Inline octave code into your perl.
I haven't tried it yet so I don't know if it's quick and easy and probably you'll find it to heavyweight for only calculating roots.
 [reply] 

