How to correctly work with Math.Pow()?

Question

I would like to find all positive integer solutions to the equation a^3 + b^3 = c^3 + d^3 where a,b,c,d are integers between 1 to 1000;

for (int a = 1; a <= 1000; ++a)
{
    for (int b = 1; b <= 1000; ++b)
    {
        for (int c = 1; c <= 1000; ++c)
        {
            for (int d = 1; d <= 1000; ++d)
            {
                if (Math.Pow(a, 3) + Math.Pow(b, 3) == Math.Pow(c, 3) + Math.Pow(d, 3))
                {
                     Console.WriteLine("{0} {1} {2} {3}", a,b,c,d);
                }
            }
      }
}

I know that d = Math.Pow(a^3 + b^3 - c^3, 1/3), so

for (int a = 1; a <= 1000; ++a)
{
    for (int b = 1; b <= 1000; ++b)
    {
        for (int c = 1; c <= 1000; ++c)
        {
            int d = (int)Math.Pow(Math.Pow(a, 3) + Math.Pow(b, 3) - Math.Pow(c, 3), 1/3);

            if (Math.Pow(a, 3) + Math.Pow(b, 3) == Math.Pow(c, 3) + Math.Pow(d, 3))
            {
                Console.WriteLine("{0} {1} {2} {3}", a,b,c,d);
            }
        }
   }
}

But the second algorithm produce much smaller number of results. Where is the bug in my code?

I try (double)1/3 but the second algorithm still gives me smaller number of results then first one

1/3 is 0 because your doing integer division. So, `d` is always 1. Try `(double)1/3` — juharr, Jun 15 '16 at 19:19
I try `(double)1/3` but the second algorithm still gives me smaller number of results then first one — Anatoly, Jun 15 '16 at 19:31

Eric Lippert · Accepted Answer · 2016-06-15T21:33:36.907

Pretty much none of these answers are right.

Your question is:

Where is the bug in my code?

If you are using double-precision arithmetic to solve a problem in the integers, you are doing something wrong. Do not use Math.Pow at all, and particularly do not use it to extract cube roots and expect that you will get an exact integer answer.

So how should you actually solve this problem?

Let's be smarter about not doing unnecessary work. Your program discovers that 1³ + 12³ = 9³ + 10³, but also that 12³ + 1³ = 10³ + 9³, and so on. If you know the first one then you can know the second one pretty easily.

So what should we do to make this more efficient?

First, b must be larger than a. That way we never waste any time figuring out that 1³ + 12³ = 12³ + 1³.

Similarly, d must be larger than c.

Now, we can also say that c and d must be between a and b. Do you see why?

Once we put these restrictions in place:

    for (int a = 1; a <= 1000; ++a)
        for (int b = a + 1; b <= 1000; ++b)
            for (int c = a + 1; c < b; ++c)
                for (int d = c + 1; d < b; ++d)
                    if (a * a * a + b * b * b == c * c * c + d * d * d)
                        Console.WriteLine($"{a} {b} {c} {d}");

Your program becomes a lot faster.

Now, there are ways to make it faster still, if you're willing to trade more memory for less time. Can you think of some ways that this program is wasting time? How many times are the same computations done over and over again? How can we improve this situation?

We might notice for example that a * a * a is computed every time through the three inner loops!

    for (int a = 1; a <= 1000; ++a)
    {
        int a3 = a * a * a;
        for (int b = a + 1; b <= 1000; ++b)
        {
            int b3 = b * b * b;
            int sum = a3 + b3;
            for (int c = a + 1; c < b; ++c)
            {
                int c3 = c * c * c;
                int d3 = sum - c3;
                for (int d = c + 1; d < b; ++d)
                    if (d3 == d * d * d)
                        Console.WriteLine($"{a} {b} {c} {d}");
            }
        }
    }

But we could be even smarter than that. For example: what if we created a Dictionary<int, int> that maps cubes to their cube roots? There are only 1000 of them! Then we could say:

    for (int a = 1; a <= 1000; ++a)
    {
        int a3 = a * a * a;
        for (int b = a + 1; b <= 1000; ++b)
        {
            int b3 = b * b * b;
            int sum = a3 + b3;
            for (int c = a + 1; c < b; ++c)
            {
                int c3 = c * c * c;
                int d3 = sum - c3;
                if (dict.HasKey(d3))
                {
                    d = dict[d3];
                    Console.WriteLine($"{a} {b} {c} {d}");
                }
            }
        }
    }

Now you don't have to compute cubes or cube roots of d; you just look up whether it is a cube, and if it is, what its cube root is.

@Anatoly The dictionary would be pre-computed before we enter this loop, you would likely set it up in the static constructor of the class to a `static readonly IReadOnlyDictionary dict;` Dictionaries are thread safe for read only operations. — Scott Chamberlain, Jun 15 '16 at 21:21
@Anatoly: We never declare it either; I presume that the reader can fill in the details of how to initialize a dictionary containing a map from cubes to roots. Don't just copy-paste code off the internet; *think about how it works*. — Eric Lippert, Jun 15 '16 at 21:27
@EricLippert Sorry, I understand. What is the complexity of your algorithms? — Anatoly, Jun 15 '16 at 21:28
@Anatoly: I invite you to work that out for yourself. It builds character. — Eric Lippert, Jun 15 '16 at 21:28
once we have the map of all the `(c,d)` pairs and `c^3 + d^3`, we can just use that directly. We don't need to generate the (a,b) pairs. Each (a,b) will already be in map. — Anatoly, Jun 16 '16 at 15:12
@Anatoly: Yes; there are ways to make this algorithm even faster provided that you are willing to trade more space for less time. Also, this equation is well studied. We know, for example, that there is a pattern that lets you generate the *rational* solutions, but there is not yet a known pattern for generating the *integer* solutions. If this subject interests you, consider reading http://www.fq.math.ca/Scanned/38-3/chamberland.pdf — Eric Lippert, Jun 16 '16 at 15:57

score 3 · Answer 2 · answered Jun 15 '16 at 20:05

Even as a double, 1/3 cannot be represented exactly with any IEEE floating point types.

So, you aren't actually computing the cube root of a value, but rather raising the value to the power of 0.333333333333333333333333333334 or some slightly off version. This introduces rounding errors that result in increasing error factors.

Consider the case with a = 995, b = 990, c = 990: your first loop would generate a match with d = 995, but your computation produces d = 994 due to rounding errors, which fails to match your condition. This will rarely match, which explains the results you are seeing.

To fix the problem, you could work entirely with floating point, but that would get messy since you have to check ranges due to representational problems even then, and only after you find a possibly suitable range, you would then have to check the integer versions. Another solution is to approach this from the mathematics standpoint rather than brute-force, although that will likely require kernel methods and get very messy.

can you provide a code please? I don't understand how to fix the problem — Anatoly, Jun 15 '16 at 20:16
How are you planning to use this? Eliminating the error will not be easy, so knowing how you plan to use this will help in recommending a solution that will work; fixing your code is out of the scope of what I can do, because it would require two layers, the first to approximate the solution entirely in floating point, then a second to attempt to settle on an integer solution based on the floating point solution. However, even to provide advice will require more info, since most of the results are of the form A+B=B+A. — Matt Jordan, Jun 15 '16 at 20:21

telenachos · Answer 3 · 2016-06-15T21:11:28.633

1

Replace

int d = (int)Math.Pow(Math.Pow(a, 3) + Math.Pow(b, 3) - Math.Pow(c, 3), 1/3);

with

int d = (int)Math.Round(Math.Pow(Math.Pow(a, 3) + Math.Pow(b, 3) - Math.Pow(c, 3), 1.0/3), 0);
if(d > 1000) continue; // Restrict solutions as in brute force algorithm

The two errors fixed are,

1/3 is computed as an integer division resulting in 0. Using 1.0/3 forces a floating point division.
Because of floating point errors the calculation may return say 3.999999994 instead of 4. When cast to an integer this is truncated to 3, effectively removing a solution. Using Math.Round(3.999999994, 0) rounds this up to 4.0 which is cast to 4 as an integer.

edited Jun 15 '16 at 21:11

answered Jun 15 '16 at 20:05

telenachos

884
6
18

Your solution gives me +24 results then the first algorithm (n = 100) – Anatoly Jun 15 '16 at 20:28
This is because in your first algorithm d is restricted to d < n. But in the second algorithm d is calculated and so solutions where d > n are included as well. To get the same number of results just exclude cases when d > n. I've updated the answer to do this :) – telenachos Jun 15 '16 at 21:20

RomCoo · Answer 4 · 2016-06-15T20:18:07.700

0

Your problem is, that the result from this line: int d = (int)Math.Pow(Math.Pow(a, 3) + Math.Pow(b, 3) - Math.Pow(c, 3), 1/3); can be a non Integer value, which would not fullfill your requirement. But by casting it to int you get a larger amount of solutions.

You should change d to a double and then just check if d is an integer value between 1 and 1000:

double d = Math.Pow(Math.Pow(a, 3) + Math.Pow(b, 3) - Math.Pow(c, 3), 1.0/3.0);
double epsilon = 0.000000001;
double dint = Math.Round(d, 0);
if (dint<=d+epsilon && dint>=d-epsilon && dint>=1-epsilon && dint<=1000+epsilon)
{
   Console.WriteLine("{0} {1} {2} {3}", a,b,c,d);
}

Edit: I added an epsilon to make sure your double comparison will work

edited Jun 15 '16 at 20:18

answered Jun 15 '16 at 19:55

RomCoo

1,868
2
23
36

what about `(double) 1/3 `? – Anatoly Jun 15 '16 at 20:02
why d is double if you cast the result to int in the first line? – Anatoly Jun 15 '16 at 20:03
do you check third line? – Anatoly Jun 15 '16 at 20:12
yeah. perhaps you should use `Math.Round()` here. But you can leave it as double – RomCoo Jun 15 '16 at 20:19

How to correctly work with Math.Pow()?

4 Answers4