Congruent Numbers and the Lousy Algorithm
A while back I posted the exciting news that every congruent number less than one trillion had been identified using some impressive hardware hackery. One of the interesting things about this was that the team working on the problem had to multiply numbers that were too large to fit in system RAM.
Just to get a little more familiar with congruent numbers, I wrote a short program using what Keith Conrad refers to as a "lousy algorithm" to generate all congruent numbers.
The essence of this algorithm is that we can generate every congruent number, sooner or later, by running through all Pythagorean triples. The three integer sides of a right triangle can be formed by the values j^{2}-k^{2}, 2jk, j^{2}+k^{2}, where:
- j > k > 0
- The smallest common factor of j and k is 1.
- j and k cannot both be odd (nor even, from the previous condition)
Once you have the values of j and k, it is simple to determine the area of the resulting right triangle, which is ½ab. To reduce the value to an actual congruent number, the final step is to divide out all the squares. For example, if you had a right triangle with sides of 6, 8, and 10, it would have an area of 24. 24 is a product of 6 and 4. Casting out the 4 leaves a congruent number of 6, which is referred to as a squarefree value.
Why are we only interested in the squarefree values? You can trivially scale up congruent numbers by scaling the sides of the triangle. For example, the triangle with sides 3, 4, and 5 generates congruent number 6. Doubling all sides gives a triangle of sides 6, 8, and 10, and the congruent number of 24. 24 is the product of 6 and 4, or the 6 and the square of 2.
Determining the squarefree number like 6 is congruent can be tricky. But once you know that 6 is congruent, all the non-squarefree products like 24 and 54 are trivial, and not quite as interesting.
Implementation
Generating the Pythagorean triples is a breeze. Looping through all values of j and k, I just have to check to make sure they have a gcd of 1, which is not too expensive. I ran my values of j and k up to 4095, which has the interesting consequence of forcing the resulting area calculation to use a 64 bit number, which slows things down quite a bit. But the place where the program actually spends all its time is attempting to cast out squares. I didn't perform any optimization of this routine, and as a result the program actually took a couple of days to run.
When a program is only going to be run once, it is usually pretty easy to forget about optimization, but in this case a few hours spent optimizing the program might have been well spent. I encourage you to give it a try and see what you can come up with. The two key routines here are one that is called to test a pair for congruence, and the support routine that casts out squares: (my apologies for the poor code formatting the accursed blog software provides)
void congruent( int j,
int k,
std::map& numbers )
{
long a = j*j -k*k;
long b = 2*j*k;
long c = j*j + k*k;
long long area = (long long) a * b /2;
long long n = cast_out_squares( area );
long long divisor = area / n;
if ( numbers.find( n ) == numbers.end() ) {
data d = { a, b, c, area, divisor };
numbers[ n ] = d;
}
}
long long cast_out_squares( long long n )
{
for ( long long i = 2 ; ; i++ ) {
while ( (n % (i*i) ) == 0 )
n = n / (i*i);
if ( (i*i) > (n/2) )
break;
}
return n;
}
Results
Despite generating over 3 million congruent numbers, as predicted, I used a lousy algorithm. If we just consider the congruent numbers under 100, here are the ones I missed: 23, 29, 37, 47, 53, 61, 62, 79, 87, 94. That's roughly a third of the actual list. Still, it is interesting to see how good a lousy algorithm is compared to paper and pencil. It's hard to imagine how long it would take someone to find a congruent number like 71, considering that the right triangle that forms it has legs of size 1320/427 and 30317/660.
Show below are my list of numbers less than 100:
n | a | b | c | area | squares |
5 | 9 | 40 | 41 | 180 | 36 |
6 | 3 | 4 | 5 | 6 | 1 |
7 | 175 | 288 | 337 | 25200 | 3600 |
13 | 104329 | 23400 | 106921 | 1220649300 | 93896100 |
14 | 63 | 16 | 65 | 504 | 36 |
15 | 15 | 8 | 17 | 60 | 4 |
21 | 7 | 24 | 25 | 84 | 4 |
22 | 99 | 4900 | 4901 | 242550 | 11025 |
30 | 5 | 12 | 13 | 30 | 1 |
31 | 2553439 | 259200 | 2566561 | 330925694400 | 10675022400 |
34 | 17 | 144 | 145 | 1224 | 36 |
38 | 1478979 | 722500 | 1646021 | 534281163750 | 14060030625 |
39 | 25 | 312 | 313 | 3900 | 100 |
41 | 369 | 800 | 881 | 147600 | 3600 |
46 | 2783 | 7056 | 7585 | 9818424 | 213444 |
55 | 13689 | 11000 | 17561 | 75289500 | 1368900 |
65 | 65 | 72 | 97 | 2340 | 36 |
69 | 8303 | 64896 | 65425 | 269415744 | 3904576 |
70 | 45 | 28 | 53 | 630 | 9 |
71 | 12945359 | 871200 | 12974641 | 5638998380400 | 79422512400 |
77 | 39375 | 15820288 | 15820337 | 311461920000 | 4044960000 |
78 | 675 | 52 | 677 | 17550 | 225 |
85 | 5929 | 6120 | 8521 | 18142740 | 213444 |
86 | 111843 | 33124 | 116645 | 1852343766 | 21538881 |
93 | 2079511 | 216600 | 2090761 | 225211041300 | 2421624100 |
95 | 2082249 | 219640 | 2093801 | 228672585180 | 2407079844 |