Robert Gerbicz wrote:
For a more realictic view, we can ask even the next:
How many odd numbers less than 2^64 has only 4k+1 primes in its decomposition to divisors, if the largest divisor is less than 2^32 (2^40 or any other number, that we will set as a goal)
Why not simulate it? For 100000 random numbers gives 595 such numbers. So we can expect approx. 595/100000*2^64~10^17 numbers to test. It's a pretty huge. Only for one number on average you should check thousands of triangle pairs.
I'll recheck your result later.
Some notes to your estimate.
First, we no need to check g = prime 4k+1, because in this case g² has only 2 decompositions to a sum of squares. So we can skip these diagonals.
Second, primes p=4k+1, 2^64/25 < p < 2^64/5 produce diagonals with only 2 divisors: p and (5 or 13 or 17)
The square of these diagonals g^2 has only 4 decompositions to a sum of squares, at the same time one of them is trivial: g² = 0² + g²
Thus a huge amount of primes p=4k+1, 2^64/25 < p < 2^64/5 produce diagonal with only 1 triangle of pairs.
The chance that this 1 triangle of pairs will produce a perfect cuboid is miserable.
And all efforts will be spend not on diagonals checking, but simply to the bust of primes.
FYI. g = 5*p, where p = x²+y² produce the next 4 decompositions:
g² = 0² + (5p)²
g² = (y² - x²)² + (2xy)²
g² = (6xy - 4y² + 4x²)² + (8xy + 3y² - 3x²)²
g² = (8xy - 3y² + 3x²)² + (6xy + 4y² - 4x²)²
Maybe someone can prove that the last three are not a perfect cuboid pairs?
First decomposition we'll skip. Other three first must be ordered in way g² = k² + l², where k < l.
After that, these three must be ordered by k acsending.
Top 2 pairs make
g² = a² + f²
g² = b² + e²
The last equation contains c and d in unknown order (maybe c < d, maybe c > d).
Thus our the only task is to check, that, for example
a² + b² = d²
we can do this check in a different ways.
a, b, d are 64-bit integers, thus we can't calculate a sum of squares directly.
First, we can check, for example that
byte(byte(a)²) + byte(byte(b)²) = byte(byte(d)²)
byte returns the last 8 bit of number.
Thus if equation is correct for 64 bits, it will be correct for the last 8 bits too.
This routine is very fast and very effective.
If check is passed, after that I apply more strong condition
integer(integer(a)²) + integer(integer(b)²) = integer(integer(d)²)
integer returns the last 32 bits.
And only after passing and this check, I apply the check, which use sqrt:
a = sqrt(d - b) sqrt(d + b)
The argument of sqrt is a number <2^65.
sqrt is real class numbers function.
Especially to point that the result must be 64 bit I use in Delphi the next trick:
a = sqrt(d - b + 0.0) sqrt(d + b + 0.0)
Thus 64 bit extended multiplied by another 64 bit extended. The result is 64 bit extended too.
Thus it's very important to sqrt correctly figured at least 48 characters and multiplication of extended must have 64 bit precision.
It's one of the weak sides of my program.
Because I don't exactly know what is an actual precision of these functions from the math library, what Delphi uses.
They say that sqrt(extended) returns extended, but how accurate this extended is calculated, I don't know.
The same about an accurancy of multiply operation.
If someone could advise something better, I'll make it.
There are some other trivial checks before the sqrt is used but I don't want to spend your time to these details.