Math.Pow() is broken

Question

And no, this does not (to my understanding) involve integer division or floating-point rounding issues.

My exact code is:

    static void Main(string[] args)
    {
        double power = (double)1.0 / (double)7.0;
        double expBase = -128.0;
        System.Console.WriteLine("sanity check: expected: -128 ^ 0.142857142857143 = -2.    actual: " + expBase + " ^ " + power + " = " + Math.Pow(expBase, power));
        System.Console.ReadLine();
    }

The output is:

sanity check: expected: -128 ^ 0.142857142857143 = -2.    actual: -128 ^ 0.14285
7142857143 = NaN

The Target Framework for this code is (according to solution properties) .NET Framework 4.0 Client Profile.

Strangely I haven't found any mention of this anywhere on the Web. Am I taking crazy pills here!?

I don't think the answer to this is what you expect: https://www.wolframalpha.com/input/?i=%28-128%29%5E%281%2F7%29 — Greg Hewgill, Jan 16 '13 at 20:44
No, the problem is one of precedence. The statement -128 ^ (1/7) could be interpreted as -(128^(1/7)) or as (-128)^(1/7). In this case, the latter is what is in use. — Kenogu Labz, Jan 16 '13 at 20:48
This is not about precedence. (-128)^(1/7) _is_ -2, but 0.14285 7142857143 != 1/7. — JLRishe, Jan 16 '13 at 20:57
@JLRishe: This is neither about precedence nor about rounding. See my answer for the mathematical reasons. — Greg Hewgill, Jan 16 '13 at 21:10
Obviously this type of calculation is not supported by System.Math. Any suggestions on how to calculate the answer? — cowlinator, Jan 16 '13 at 22:29
Does anyone know of a math library with an alternative to Math.Pow()? — cowlinator, Jan 22 '13 at 20:18
Think about this in terms of mathematical definitions. What you're trying to find is the 7th root of -128. However, this is exactly equivalent to the negated 7th root of 128. Using a positive radicand is guaranteed to give a valid value for all 'standard' roots. Instead of trying to find a different math library, just restructure the problem to make it work with your code: in this case, `-Math.Pow(-expBase, power)`. Do note that this only works in this case because your radix is 7, an odd number. An even-radix root would produce a complex number, which has no standard form in C#, AFAIK. — Kenogu Labz, Jan 30 '13 at 18:16

score 18 · Accepted Answer · edited Jun 20 '20 at 09:12

18

Seems to be exactly as specified; from the Math.Pow() remarks section on Pow(x,y);

Parameters
x < 0 but not NegativeInfinity; y is not an integer, NegativeInfinity, or PositiveInfinity.

Result
NaN

edited Jun 20 '20 at 09:12

Community

1
1

answered Jan 16 '13 at 20:46

Joachim Isaksson

176,943
25
281
294

Stephen Canon · Answer 2 · 2013-01-16T20:52:48.080

13

Joachim's answer explains that pow is behaving according to its specification.

Why is pow( ) specified that way? Because 1.0/7.0 is not equal to 1/7. You are asking for the 0.14285714285714285 power of -128.0, and there is no real number with that property, so the result is correctly NaN. For all odd n != 1, 1.0/(double)n is not exactly representable, so you can't compute the nth root of x by using pow(x, 1.0/(double)n). Therefore pow(x, y) is specified to return NaN for negative x and non-integer y -- there is no appropriate real result for any of those cases.

edited Jan 16 '13 at 20:52

answered Jan 16 '13 at 20:47

Stephen Canon

103,815
19
183
269

It should return the representable value closest to (-128)**(.1428…). :-) – Eric Postpischil Jan 16 '13 at 21:11
@EricPostpischil: Somehow I doubt the questioner expects the answer `1.8019377358048383`; it's really not very close at all. – Stephen Canon Jan 16 '13 at 21:15
2

But think how much it would add to the questions on Stack Overflow. Instead of explaining how .1+.2 gives .30…01 due to correctly rounded arithmetic, we could explain how 1.8019… is the correctly rounded result of `pow(-128, 1/7.)`. – Eric Postpischil Jan 16 '13 at 21:45
I guess I fail to understand how it has anything to do with 1/7 not being exactly representable, when 1/7 is still not representable if you use a positive base i.e. `double power = (double)1.0 / (double)7.0;' 'double expBase = 128.0;` results in exactly 2.0 – cowlinator Jan 16 '13 at 22:06
@user1698736: the difference is that `pow(128.0, 1.0/7.0)` is a real number; specifically, it's `1.999999999999999923045204...`, which happens to round to `2.0` in double precision. `pow(-128.0, 1.0/7.0)`, by contrast, is not real-valued. Just like `sqrt(-1)`, which is also not real-valued, it returns NaN. – Stephen Canon Jan 17 '13 at 10:59
@user1698736: (basically, I agree with your comment on Greg's answer--if there were a simple real number result, it would be useful for `pow` produce it, even though it wouldn't be the principle value. However, because of the fact that 1/7 is not representable there is no simple real number result; all of the solutions have non-zero imaginary part). – Stephen Canon Jan 17 '13 at 11:25
I think the fundamental failing is the lack of an `NthRoot(double,int)` function. Even when calculating roots using logarithms, computing `2^(lg(128)/7)` is more semantically precise and numerically accurate than `2^(lg(128)*(1/7)`. – supercat May 07 '14 at 19:11
1

@supercat: IEEE-754 recommended that languages provide the `rootn` function for this purpose in 2008; it will take a few years to find its way into C, and a few years after that to find its way into most other languages. – Stephen Canon May 07 '14 at 23:32
@StephenCanon: With or without such a recommendation, it would seem like an obvious thing for a framework designer to include. Validate `n` and the argument, look for "nice" cases, and compute `x^(1/n)`. If one is feeling really fancy, check whether the result is too high or low and nudge it accordingly (I think it's possible to compute a maximum error term, and multiplies should be cheap compared to the cost of computing a log and exponentiation. – supercat May 08 '14 at 00:14
@supercat: It’s nontrivial to implement correctly, if on cares about the little details like getting the inexact flag right. But yes, it’s a nice feature to have. – Stephen Canon May 08 '14 at 00:32
@StephenCanon: Implementing it to bit-level precision may be difficult, but if the API contract specified it as being at least as precise as x^(1/n), code which used the API rather than x^(1/n) would benefit if the API implementation was later improved. By contrast, code written with x^(1/n) will be forever limited to any imprecision of that approach. Likewise, if code wants to compute the sine of an angle that's some number of full cycles, using an API for that would be cleaner than multiplying the angle by 2*pi, though the latter case has a certain irony. – supercat May 08 '14 at 12:59
@supercat: I absolutely agree; that's why the committee added it in 2008 (along with sinpi). Library implementors are a conservative bunch, though, and tend not to add new features until they're reasonably certain that they will be standardized (and in what form they will be). It's a slow process. – Stephen Canon May 08 '14 at 13:05
@StephenCanon: Given that the IEEE-754 goes back decades, and people have been computing things like `signal = sin(time*freq*(2*pi))` for a *looong* time, the pace seems positively glacial. Even if it took awhile for implementations to catch up to recommendations, I would think the need for things like like whole-circle-based trig, roots, and transitive relational operators would have become apparent in the 1980s. What has the Committee been doing in the last 20 years? – supercat May 08 '14 at 13:09

Greg Hewgill · Answer 3 · 2013-01-16T21:13:34.553

3

The problem here is the mathematical definition of "seventh root" is a multivalued function. While it is true that

(-2)⁷ = -128

this does not mean that -2 is the only answer to (-128)^1/7. In the complex plane, the seventh-root function is multivalued and has many possible answers (just as the square root of 4 can be considered to be either +2 or -2, but +2 is the usual answer).

In order to simplify the mathematical handling of such expressions, the principal value is chosen by convention for the function in question so that the function becomes single-valued. In the case of seventh-root, the principal value is that given by Wolfram Alpha for (-128)^1/7.

The Math.Pow() function in C# attempts to return the principal value for the pow function. When the principal value of the result would be a complex number, it returns NaN.

edited Jan 16 '13 at 21:13

answered Jan 16 '13 at 21:08

Greg Hewgill

951,095
183
1,149
1,285

That is not the problem. sqrt is in generally multivalued, but we define `sqrt` to return the positive root. `pow` could be similarly defined. The problem in this specific case is that `pow` is not being passed 1/7 for the power. It is being passed a `double` that is near 1/7, and it has no way to know that 1/7 was intended. – Eric Postpischil Jan 16 '13 at 21:14
@EricPostpischil: There is no reason, even with exact calculations, why (-128)^(1/7) should return -2. That answer would be on a branch that wouldn't be considered a principal branch. Wolfram Alpha agrees and shows a result that is nowhere near -2. – Greg Hewgill Jan 16 '13 at 21:23
@EricPostpischil Not a difference in principle, `1.0/7.0` is still a rational number, so the power is a finitely branched covering of the plane. – Daniel Fischer Jan 16 '13 at 21:24
None of which get the intended answer. I.e., the problem here is not that the function is multivalued. The problem is that information has been lost before `pow` has been called. There is **no** definition of a function that gives the correct answer given incorrect inputs (in general). – Eric Postpischil Jan 16 '13 at 21:39
Math.Sqrt() and Math.Pow() often should have multiple values, but only return one. I figured that the solution was chosen on what was a simple, real number. The standard windows calculator seems to be able to choose a human-readable solution to -128^(1/7)... I'm surprised Math.Pow() does not. – cowlinator Jan 16 '13 at 22:26
@user1698736: Ask Windows calculator for `(-128)^(1/7)`, not `-128^(1/7)`. – Eric Postpischil Jan 17 '13 at 14:50
Interesting. The fact is that for `n` a positive **odd** number, there are two different definitions of the `n`-th root of `x` or `Pow(x, 1.0 / n)`. When working with complex numbers, it is natural to take the principal branch defined by convention, and this is like Greg Hewgill describes. But when working exclusively with real numbers, it is common to take the negative real solution. See [cube root on Wikipedia](http://en.wikipedia.org/wiki/Cube_root) where the section _Formal definition_ gives both conventions. @EricPostpischil is right, Windows calculator does `(-128)^(1/7)` as `-2`. – Jeppe Stig Nielsen Jan 19 '13 at 19:27
_(continued)_ Now, a `double` represents only real numbers, not complex. But the problem for `Pow` is to determine if its second argument is a rational number with **odd** denominator, or not. That's quite impossible. For any odd `n>1`, the binary expansion of `1.0 / n` is periodic and non-terminating. So `Pow` simply says `NaN` for any `Pow(x, q)` where `x` is negative and `q` is not whole. Many calculators try to see if `q` looks like a rational number with a (small) odd denominator. – Jeppe Stig Nielsen Jan 19 '13 at 19:33
@JeppeStigNielsen: I would suggest that taking integer (n) powers and roots of some number (x) should be considered as being totally different operations from the raising a number to a real-valued exponent. Raising *anything* to the *integer*-zero power yields the multiplicative identity. Otherwise, x*x^n = x^(n+1). The nth root y of some value x is an object in the same domain of x such that y^n=x. The fact that when x is positive, x^n = exp(log(x)*n) and the nth root of x is exp(log(x)/n) provide convenient ways of calculating the nth root in cases where they work, but don't *define* it. – supercat May 07 '14 at 19:08

score 1 · Answer 4 · edited Nov 03 '22 at 19:21

1

A fractional power of a negative real number is a complex number (see Math Forum for a detailed explanation).

edited Nov 03 '22 at 19:21

Glorfindel

21,988
13
81
109

answered Jan 16 '13 at 20:49

Simon

10,679
1
30
44

Not always. Given an odd whole number y and a real number x, x^(1/y) will always be real, even if x is negative. – JLRishe Jan 16 '13 at 21:05
@JLRishe Not if you interpret a power function `\z -> z^e` as an analytic function. The analytic continuations of that from the positive half-line take nonreal values on the negative half-line for `e = 1/y`, `y > 1` integer. – Daniel Fischer Jan 16 '13 at 21:20

cowlinator · Answer 5 · 2013-01-24T22:29:32.713

I fixed Math.Pow().

It now has a larger accepted domain (i.e. for Parameters: x < 0 but not NegativeInfinity; y is a fraction with a numerator of 1 and an odd denominator), and returns a real number result for the new domain area.

In other words, (-128)^(1/7) returns -2.

Note: due to double-float precision limitations, it will work for most, but not all, fractional exponents.

Below is the code for the wrapper for Math.Pow() I wrote.

public class MathUtil
{
    /// <summary>
    /// Wrapper for Math.Pow()
    /// Can handle cases like (-8)^(1/3) or  (-1/64)^(1/3)
    /// </summary>
    public static double Pow(double expBase, double power)
    {
        bool sign = (expBase < 0);
        if (sign && HasEvenDenominator(power)) 
            return double.NaN;  //sqrt(-1) = i
        else
        {
            if (sign && HasOddDenominator(power))
                return -1 * Math.Pow(Math.Abs(expBase), power);
            else
                return Math.Pow(expBase, power);
        }
    }

    private static bool HasEvenDenominator(double input)
    {
        if(input == 0)
            return false;
        else if (input % 1 == 0)
            return false;

        double inverse = 1 / input;
        if (inverse % 2 < double.Epsilon)
            return true;
        else
            return false;
    }

    private static bool HasOddDenominator(double input)
    {
        if (input == 0)
            return false;
        else if (input % 1 == 0)
            return false;

        double inverse = 1 / input;
        if ((inverse + 1) % 2 < double.Epsilon)
            return true;
        else
            return false;
    }
}

Why do you use HasOddDenominator. Cant you just say: !HasEvenDenominator ? — Eystein Bye, Oct 07 '13 at 11:55
Thx,but the code cant be use `x=-10; x ^ (2/3)` and the result is positive. — lindexi, Apr 21 '17 at 10:11
(-10)^(2/3) is a multivalued function, and one of the results is positive (in RL). So what is the problem? — cowlinator, Oct 18 '17 at 20:23

score 0 · Answer 6 · edited Aug 31 '23 at 22:00

For 1.0/3 != 1/3 ,I use Rational that can Accurately represent 1/3 in Microsoft.SolverFoundation.Common . See https://msdn.microsoft.com/en-us/library/microsoft.solverfoundation.common.rational%28v=vs.93%29.aspx?f=255&MSPPError=-2147217396

And I can catch the odd root link 1/3 for I can get the Denominator.

If I get a root is ax that I use the code get the Numerator and Denominator.

            var at = (double)ax.Numerator;
            var down = (double)ax.Denominator;

Rational can make the 2/6=1/3.

But the Rational.Pow cant calculate powerBase isnt positive.

I find the powerBase isnt positive and the Denominator is even and Numerator is odd.

             if (at % 2 == 1 && down % 2 == 0)
            {
                return Double.NaN;
            }

If the Denominator is odd ,I use x = x * -1

            if (at % 2 == 1 && down % 2 == 1)
            {
                x = Math.Pow(x, (int)at);
                x = x * -1;
                return -1 * Math.Pow(x, 1.0 / (int)down);
            }

If the Numerator is even ,the pow of Numerator make powerBase is positive.

Like pow(x,2/3), if x isnt positive ,it will be positive when use pow(x,2)

            x = Math.Pow(x, (int)at);
            return Math.Pow(x, 1.0 / (int)down);

The code you can use

       if (x < 0)
        {                
            var at = (double)ax.Numerator;
            var down = (double)ax.Denominator;

            if (at % 2 == 1 && down % 2 == 0)
            {
                return Double.NaN;
            }
            if (at % 2 == 1 && down % 2 == 1)
            {
                x = Math.Pow(x, (int)at);
                x = x * -1;
                return -1 * Math.Pow(x, 1.0 / (int)down);
            }
            x = Math.Pow(x, (int)at);
            return Math.Pow(x, 1.0 / (int)down);
        }
        else
        {
            return Math.Pow(x, a);
        }

Math.Pow() is broken

6 Answers6

Linked