macOS Numerical errors not seen under Debian

[cosmicassassin]

Just bought a Mac for the first time.
Running OS X 10.5.3 with g++ version 4.0.1 on an Intel iMac.
This code gives incorrect solutions when I run it on this computer but gives the correct solution on my Debian box.

#include <iostream>
#include <iomanip>
#include <math.h>

int main(){
  unsigned long int sum1, sum2, sum3, whereIam(6), lastHouse;
  int numberHouses(0);

  while(numberHouses != 10){
    sum1 = (whereIam*(whereIam-1))/2;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*whereIam*whereIam)); 
    sum3 = (whereIam*(whereIam+1))/2;
    sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    while(sum1 > sum2){
      lastHouse++;
      sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    }
    
    if(sum1 == sum2){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
    whereIam++;
  }
}

Just wondering what I'm missing on the Mac. Any help would be awesome.

 

[yeroen]

I compiled your code with g++ 4.0.1 on my intel mac pro and ran it. Here are the results I got:

         6         8
        35        49
       204       288
      1189      1681
      6930      9800
     66568     16512
     67382     22152
    262208      8192
    376663     16798
    468401     26689

Without attempting to grok exactly what it is you're trying to do with this code, I can say the 'unsigned long' declarations immediately raise a red flag and I'm betting that's the reason for the discrepancy. By varying the integral type of the unsigned long declarations, I get much different results (and always at the 6th iteration). Declaring these variables plain old int throws the code into an infinite loop, and an unsigned long long declaration yields different answers.

On your debian box, what is sizeof (unsigned long) ? Is it 8 bytes instead of 4?

For portability between the two, #include <stdint.h> and replace your unsigned long declarations with uint32_t or uint64_t, whichever one gives you the results you want.

 

[gnasher729]

Just bought a Mac for the first time.
Running OS X 10.5.3 with g++ version 4.0.1 on an Intel iMac.
This code gives incorrect solutions when I run it on this computer but gives the correct solution on my Debian box.

#include <iostream>
#include <iomanip>
#include <math.h>

int main(){
  unsigned long int sum1, sum2, sum3, whereIam(6), lastHouse;
  int numberHouses(0);

  while(numberHouses != 10){
    sum1 = (whereIam*(whereIam-1))/2;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*whereIam*whereIam)); 
    sum3 = (whereIam*(whereIam+1))/2;
    sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    while(sum1 > sum2){
      lastHouse++;
      sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    }
    
    if(sum1 == sum2){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
    whereIam++;
  }
}

Just wondering what I'm missing on the Mac. Any help would be awesome.

This looks awfully obfuscated and overly complicated. Lets first change to a for loop:

  for(whereIam = 6, numberHouses = 0; numberHouses < 10; whereIam++){
    sum1 = (whereIam*(whereIam-1))/2;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*whereIam*whereIam)); 
    sum3 = (whereIam*(whereIam+1))/2;
    sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    while(sum1 > sum2){
      lastHouse++;
      sum2 = lastHouse*(lastHouse+1)/2 - sum3;
    }
    
    if(sum1 == sum2){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
  }

Now we don't subtract sum3 from sum2 anymore; you will see why soon:

  for(whereIam = 6, numberHouses = 0; numberHouses < 10; whereIam++){
    sum1 = (whereIam*(whereIam-1))/2;
    sum3 = (whereIam*(whereIam+1))/2;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*whereIam*whereIam)); 
    sum2 = lastHouse*(lastHouse+1)/2;
    while(sum1 > sum2 - sum3){
      lastHouse++;
      sum2 = lastHouse*(lastHouse+1)/2;
    }
    
    if(sum1 == sum2 - sum3){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
  }

Now I move sum3 to the other side of the comparison and don't store sum2 at all:

  for(whereIam = 6, numberHouses = 0; numberHouses < 10; whereIam++){
    sum1 = (whereIam*(whereIam-1))/2;
    sum3 = (whereIam*(whereIam+1))/2;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*whereIam*whereIam)); 
    while(sum1 + sum3 > lastHouse*(lastHouse+1)/2)
      lastHouse++;
    
    if(sum1 + sum3 == lastHouse*(lastHouse+1)/2){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
  }

Well, what is sum1 + sum3? It is equal to whereIam * whereIam:

  for(whereIam = 6, numberHouses = 0; numberHouses < 10; whereIam++){
    sum = whereIam*whereIam;
    lastHouse = static_cast<int>(-0.5 + 0.5*sqrt(1 + 8*sum)); 
    while(sum > lastHouse*(lastHouse+1)/2)
      lastHouse++;
    
    if(sum == lastHouse*(lastHouse+1)/2){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
  }

Now wait a second... What are you doing here? You are checking in some rather complex way whether x (x+1) is twice a square. That would happen very rarely, and you are searching for ten solutions. I suspect that things go wrong rather soon when whereIam > 23170 and the expression in the square root exceeds 32 bits. I guess you are compiling for a 64 bit machine on the Debian box. As evidence, see the output by another poster: The second column should be roughly equal to the first column times 1.4142; after the fifth value it is nowhere near.

I'd first switch to 64 bit arithmetic, using unsigned long long instead of unsigned long. The static_cast should cast to unsigned long long as well. The expensive square root calculation is actually unnecessary: Just initialise lastHouse to 0 before the loop starts (because the square root will get bigger and bigger results). So:

  for(whereIam = 6, lastHouse = 0, numberHouses = 0; numberHouses < 10; whereIam++){
    sum = 2*whereIam*whereIam;
    while(sum > lastHouse*(lastHouse+1))
      lastHouse++;
    
    if(sum == lastHouse*(lastHouse+1)){
      numberHouses++;
      std::cout << std::setw(10) << whereIam << std::setw(10) << lastHouse << std::endl;
    }
  }

 

[gnasher729]

By the way, you can find all the solutions very very quickly with a bit of maths and without any searching. You want to find pairs (x, y) such that x * (x+1) = 2 * y^2. The following algorithm does it:

Start with a = b = 1. Then repeat the following as often as you like:
One solution is x = a^2, y = ab.
Let c = a + 2b, d = a + b.
Another solution is x = 2 * d^2, y = cd.
Let a = c + 2d, b = c + d and repeat.

You get the following solutions:
a = 1, b = 1, x = 1, y = 1. 1*2 = 2*1*1 = 2.
c = 3, d = 2, x = 8, y = 6. 8*9 = 2*6*6 = 72.
a = 7, b = 5, x = 49. y = 35. 49*50 = 2*35*35 = 2450.
c = 17, d = 12, x = 288, y = 204. 288*289 = 2*204*204 = 83232.
a = 41, b = 29, x = 1681, y = 1189
c = 99, d = 70, x = 9800, y = 6930
and so on. There are no other solutions.

The last solution you were looking for is x = 46,611,179 and y = 65,918,161. A rather large solution (18 solutions down the list) is a = 63,018,038,021 and b = 44,560,482,149. If you calculate x = a^2 (about 3.97 * 10^21), y = ab (about 2.81*10^21) you get x*(x+1) = 2*y*y exactly ; both products are about 1.58*10^43.

 

[Mac Player]

A little more info about how to solve second-degree diphantine equations : http://www.alpertron.com.ar/METHODS.HTM

 

[yeroen]

I'm impressed that you were able to parse the spaghetti for it's intent.