False Positives on an Integer Lattice Leaderboard

 User            Language      Score
 =========================================
 Ell             C++11         293,619,555
 feersum         C++11         100,993,667
 Ell             C++11          78,824,732
 Geobits         Java           27,817,255
 Ell             Python         27,797,402
 Peter Taylor    Java                2,468
 &lt;reference&gt;     Julia                 530


Background

When working on a 2-D grid of integer coordinates, you sometimes want to know whether two vectors (with integer components) have the same magnitude. Of course, in Euclidean geometry the magnitude of a vector (x,y) is given by

√(x² + y²)


So a naive implementation might compute this value for both vectors and compare the results. Not only does that incur an unnecessary square root computation, it also causes trouble with floating point inaccuracies, which might yield false positives: vectors whose magnitudes are different, but where the significant digits in the floating point representation are all identical.

For the purposes of this challenge, we define a false positive as a pair of coordinate pairs (a,b) and (c,d) for which:


Their squared magnitude is different when represented as 64-bit unsigned integers.
Their magnitude is identical when represented as 64-bit binary floating point number and computed via a  64-bit square root (as per IEEE 754).


As an example, using 16-bit representations (instead of 64), the smallest1 pair of vectors which yields a false positive is

(25,20) and (32,0)


Their squared squared magnitudes are 1025 and 1024. Taking the square root yields 

32.01562118716424 and 32.0


But in 16-bit floats both of these get truncated to 32.0.

Likewise, the smallest2 pair yielding a false positive for 32-bit representations is

(1659,1220) and (1951,659)



1"Smallest" as measured by their 16-bit floating point magnitude.
2"Smallest" as measured by their 32-bit floating point magnitude.


Lastly, here are a handful of valid 64-bit cases:

 (51594363,51594339) and (54792160,48184783)
 (54356775,54353746) and (54620742,54088476)
 (54197313,46971217) and (51758889,49645356)
 (67102042,  956863) and (67108864,       6) *


* The last case is one of several with the smallest possible magnitude for 64-bit false positives.

The Challenge

In less than 10,000 bytes of code, using a single thread, you're to find as many false positives for 64-bit (binary) floating point numbers  in the coordinate range 0 ≤ y ≤ x (that is, only within the first octant of the Euclidian plane) such that x² + y² ≤ 253 within 10 minutes. If two submissions tie for the same number of pairs, the tie breaker is the actual time taken to find the last of those pairs.

Your program must not use more than 4 GB of memory at any time (for practical reasons).

It must be possible to run your program in two modes: one which outputs every pair as it finds it, and one which only outputs the number of found pairs at the end. The first will be used to verify the validity of your pairs (by looking at some sample of outputs) and the second will be used for actually timing your submission. Note that the printing must be the only difference. In particular, the counting program may not hardcode the number of pairs it could find. It must still perform the same loop that would be used to print all the numbers, and only omit the printing itself!

I will test all submissions on my Windows 8 laptop, so please ask in the comments if you want to use some not-too-common language.

Note that pairs must not be counted twice under switching of the first and second coordinate pair.

Note also that I will run your process via a Ruby controller, which will kill your process if it hasn't finished after 10 minutes. Make sure to output the number of pairs found by then. You can either keep track of time yourself and print the result just before the 10 minutes elapse, or you just output the number of pairs found sporadically, and I will take the last such number as your score.
Python, 27,797,402
Just to set the bar a little higher...
from sys import argv
verbose = len(argv) &gt; 1
found = 0
for x in xrange(67108864, 94906266):
    found += 1
    if verbose:
        print &quot;(%d, 0) (%d, 1)&quot; % (x, x)
print found

It's easy to verify that for all 67,108,864 &lt;= x &lt;= 94,906,265 = floor(sqrt(253)) the pairs (x, 0) and (x, 1) are false positives.
Why it works: 67,108,864 = 226. Therfore, all numbers x in the above range are of the form 226 + x' for some 0 &lt;= x' &lt; 226.
For all positive e, (x + e)2 = x2 + 2xe + e2 = x2 + 227e + 2x'e + e2.
If we want to have  (x + e)2 = x2 + 1 we need at least 227e &lt;= 1, that is, e &lt;= 2-27
However, since the mantissa of double-precision floating point numbers is 52-bit wide, the smallest e such that x + e &gt; x is e = 226 - 52 = 2-26.
In other words, the smallest representable number greater than x is x + 2-26 while the result of sqrt(x2 + 1) is at most x + 2-27.
Since the default IEEE-754 rounding mode is round-to-nearest; ties-to-even, it will always round to x and never to x + 2-26 (where the tie-break is really only relevant for x = 67,108,864, if at all. Any larger number will round to x regardless).

C++, 75,000,000+
Recall that 32 + 42 = 52.
What this means is that the point (4, 3) lies on the circle of radius 5 centered around the origin.
Actually, for all integers n, (4n, 3n) lies on such circle of radius 5n.
For large enough n (namely, such that 5n &gt;= 226), we already know a false-positive for all points on this circle: (5n, 1).
Great! That's another 27,797,402 / 5 free false-positive pairs right there!
But why stop here? (3, 4, 5) is not the only such triplet.
This program looks for all positive integer triplets (a, b, c) such that a2 + b2 = c2, and counts false-positives in this fashion.
It gets to 70,000,000 false-positives pretty fast but then slows down considerably as the numbers grow.
Compile with g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3.
To verify the results, add -DVERIFY (this will be notably slower.)
Run with flspos. Any command-line argument for verbose mode.
#include &lt;cstdio&gt;
#include &lt;cmath&gt;
#include &lt;cfloat&gt;

using namespace std;

/* Make sure we actually work with 64-bit precision */
#if defined(VERIFY) &amp;&amp; FLT_EVAL_METHOD != 0 &amp;&amp; FLT_EVAL_METHOD != 1
#   error &quot;invalid FLT_EVAL_METHOD (did you forget `-msse2 -mfpmath=sse'?)&quot;
#endif

template &lt;typename T&gt; inline long long sqr(T x) { return 1ll * x * x; }
template &lt;typename T&gt;
T gcd(T a, T b) { while (b) { T t = a; a = b; b = t % b; } return a; }

int main(int argc, const char* argv[]) {
    const bool verbose = argc &gt; 1;

    const int min = 1 &lt;&lt; 26;
    const int max = sqrt(1ll &lt;&lt; 53);

    int found = 0;
    auto add = [=, &amp;found] (int x1, int y1, int x2, int y2) {
#ifdef VERIFY
        auto n1 = sqr(x1) + sqr(y1), n2 = sqr(x2) + sqr(y2);
        if (n1 == n2 || sqrt(n1) != sqrt(n2)) {
            fprintf(stderr, &quot;Wrong false-positive: (%d, %d) (%d, %d)\n&quot;,
                    x1, y1, x2, y2);
            return;
        }
#endif
        if (verbose) printf(&quot;(%d, %d) (%d, %d)\n&quot;, x1, x2, y1, y2);
        ++found;
    };

    for (int x = min; x &lt;= max; ++x) add(x, 0,    x, 1);

    for (int a = 1; a &lt; max; ++a) {
        auto a2b2 = sqr(a) + 1;
        for (int b = 1; b &lt;= a; a2b2 += 2 * b + 1, ++b) {
            int c = sqrt(a2b2);
            if (a2b2 == sqr(c) &amp;&amp; gcd(a, b) == 1) {
                int max_c = max / c;
                for (int n = (min + c - 1) / c; n &lt;= max_c; ++n)
                    add(n * a, n * b,    n * c, 1);
            }
        }

        if (a % 512 == 0) printf(&quot;%d\n&quot;, found), fflush(stdout);
    }

    printf(&quot;%d\n&quot;, found);
}