Tabla de clasificación

 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
 <reference>     Julia                 530

Antecedentes

Cuando se trabaja en una cuadrícula 2D de coordenadas enteras, a veces desea saber si dos vectores (con componentes enteros) tienen la misma magnitud. Por supuesto, en geometría euclidiana, la magnitud de un vector (x,y)viene dada por

√(x² + y²)

Entonces, una implementación ingenua podría calcular este valor para ambos vectores y comparar los resultados. Esto no solo genera un cálculo innecesario de raíz cuadrada, sino que también causa problemas con las imprecisiones de coma flotante, lo que puede generar falsos positivos: vectores cuyas magnitudes son diferentes, pero donde los dígitos significativos en la representación de coma flotante son todos idénticos.

Para los propósitos de este desafío, definimos un falso positivo como un par de pares de coordenadas (a,b)y (c,d)para el cual:

• Su magnitud al cuadrado es diferente cuando se representan como enteros sin signo de 64 bits.
• Su magnitud es idéntica cuando se representa como un número de punto flotante binario de 64 bits y se calcula a través de una raíz cuadrada de 64 bits (según IEEE 754 ).

Como ejemplo, usando representaciones de 16 bits (en lugar de 64), el ¹ par más pequeño de vectores que produce un falso positivo es

(25,20) and (32,0)

Sus magnitudes al cuadrado al cuadrado son 1025y 1024. Tomando los rendimientos de raíz cuadrada

32.01562118716424 and 32.0

Pero en flotantes de 16 bits, ambos se truncan 32.0.

Del mismo modo, el par más pequeño de ^{2 que} produce un falso positivo para representaciones de 32 bits es

(1659,1220) and (1951,659)

_{¹ "Más pequeño" medido por su magnitud de coma flotante de 16 bits.
² "Más pequeño" medido por su magnitud de coma flotante de 32 bits.}

Por último, aquí hay un puñado de casos válidos de 64 bits:

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

* El último caso es uno de varios con la menor magnitud posible para falsos positivos de 64 bits.

El reto

En menos de 10,000 bytes de código, usando un solo hilo, encontrará tantos falsos positivos para números de punto flotante (binario) de 64 bits en el rango de coordenadas 0 ≤ y ≤ x(es decir, solo dentro del primer octante del plano euclidiano) tal que dentro de los 10 minutos . Si dos envíos empatan para la misma cantidad de pares, el desempate es el tiempo real necesario para encontrar el último de esos pares.x² + y² ≤ 253

Su programa no debe usar más de 4 GB de memoria en ningún momento (por razones prácticas).

Debe ser posible ejecutar su programa en dos modos: uno que genera cada par a medida que lo encuentra, y otro que solo genera el número de pares encontrados al final. El primero se usará para verificar la validez de sus pares (observando alguna muestra de resultados) y el segundo se usará para realmente sincronizar su envío. Tenga en cuenta que la impresión debe ser la única diferencia. En particular, el programa de conteo puede no codificar el número de pares que podría encontrar. ¡Todavía debe realizar el mismo bucle que se usaría para imprimir todos los números, y solo omitir la impresión misma!

Probaré todas las presentaciones en mi computadora portátil con Windows 8, así que pregunte en los comentarios si desea usar un lenguaje no muy común.

Tenga en cuenta que los pares no se deben contar dos veces al cambiar el primer y el segundo par de coordenadas.

Tenga en cuenta también que ejecutaré su proceso a través de un controlador Ruby, que matará su proceso si no ha terminado después de 10 minutos. Asegúrese de generar el número de pares encontrados para entonces. Puede realizar un seguimiento del tiempo usted mismo e imprimir el resultado justo antes de que transcurran los 10 minutos, o simplemente generar el número de pares encontrados esporádicamente, y tomaré el último número como su puntaje.

C ++, 275,000,000+

Nos referiremos a pares cuya magnitud es representable con precisión, como (x, 0) , como pares honestos y a todos los demás pares como pares deshonestos de magnitud m , donde m es la magnitud informada erróneamente del par. El primer programa en la publicación anterior usó un conjunto de parejas estrechamente relacionadas de pares honestos y deshonestos:
(x, 0) y (x, 1) , respectivamente, para x suficientemente grande. El segundo programa usó el mismo conjunto de pares deshonestos pero extendió el conjunto de pares honestos al buscar todos los pares honestos de magnitud integral. El programa no finaliza en diez minutos, pero encuentra la gran mayoría de sus resultados muy pronto, lo que significa que la mayor parte del tiempo de ejecución se desperdicia. En lugar de seguir buscando pares honestos cada vez menos frecuentes, este programa usa el tiempo libre para hacer la siguiente cosa lógica: extender el conjunto de pares deshonestos .

De la publicación anterior, sabemos que para todos los enteros suficientemente grandes r , sqrt (r ² + 1) = r , donde sqrt es la función de raíz cuadrada de punto flotante. Nuestro plan de ataque es encontrar pares P = (x, y) de modo que x ² + y ² = r ² + 1 para algún número entero suficientemente grande r . Eso es lo suficientemente simple como para hacerlo, pero la ingenua búsqueda de pares individuales es demasiado lenta para ser interesante. Queremos encontrar estos pares a granel, tal como lo hicimos para los pares honestos en el programa anterior.

Sea { v , w } un par de vectores ortonormales. Para todos los escalares reales r , || r v + w || ² = r ² + 1 . En ℝ ² , este es un resultado directo del teorema de Pitágoras:

Estamos buscando los vectores v y w de modo que exista un número entero r para el cual x e y también sean números enteros. Como nota al margen, tenga en cuenta que el conjunto de pares deshonestos que utilizamos en los dos programas anteriores fue simplemente un caso especial de esto, donde { v , w } era la base estándar de ℝ ² ; Esta vez deseamos encontrar una solución más general. Aquí es donde los trillizos pitagóricos (trillizos enteros (a, b, c) que satisfacen a ² + b ² = c ², que utilizamos en el programa anterior) hacen su regreso.

Sea (a, b, c) un triplete pitagórico. Los vectores v = (b / c, a / c) y w = (-a / c, b / c) (y también
w = (a / c, -b / c) ) son ortonormales, como es fácil de verificar . Resulta que, para cualquier elección del triplete pitagórico, existe un número entero r tal que x e y son números enteros. Para probar esto, y para encontrar efectivamente r y P , necesitamos una pequeña teoría de números / grupos; Voy a ahorrar los detalles. De cualquier manera, supongamos que tenemos nuestra integral r , x e y . Todavía nos faltan algunas cosas: necesitamos rser lo suficientemente grande y queremos un método rápido para derivar muchos más pares similares de este. Afortunadamente, hay una manera simple de lograr esto.

Tenga en cuenta que la proyección de P en v es r v , por lo tanto, r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , todo esto para decir que xb + ya = rc . Como resultado, para todos los enteros n , (x + bn) ² + (y + an) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1. En otras palabras, la magnitud al cuadrado de los pares de la forma
(x + bn, y + an) es (r + cn) ² + 1 , ¡que es exactamente el tipo de pares que estamos buscando! Para n suficientemente grande , estos son pares deshonestos de magnitud r + cn .

Siempre es bueno mirar un ejemplo concreto. Si tomamos el triplete pitagórico (3, 4, 5) , entonces en r = 2 tenemos P = (1, 2) (puede verificar que (1, 2) · (4/5, 3/5) = 2 y, claramente, 1 ² + 2 ² = 2 ² + 1. ) Sumar 5 a r y (4, 3) a P nos lleva a r '= 2 + 5 = 7 y P' = (1 + 4, 2 + 3) = (5, 5) . Lo y he aquí, 5 ² + 5 ² = 7 ² + 1. Las siguientes coordenadas son r '' = 12 y P '' = (9, 8) , y nuevamente, 9 ² + 8 ² = 12 ² + 1 , y así sucesivamente ...

Una vez que r es lo suficientemente grande, comenzamos a obtener pares deshonestos con incrementos de magnitud de 5 . Eso es aproximadamente 27,797,402 / 5 pares deshonestos.

Así que ahora tenemos muchos pares deshonestos de magnitud integral. Podemos unirlos fácilmente con los pares honestos del primer programa para formar falsos positivos, y con el debido cuidado también podemos usar los pares honestos del segundo programa. Esto es básicamente lo que hace este programa. Al igual que el programa anterior, también encuentra la mayoría de sus resultados muy pronto, llega a 200,000,000 de falsos positivos en unos pocos segundos, y luego se ralentiza considerablemente.

Compilar con g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Para verificar los resultados, agregue -DVERIFY(esto será notablemente más lento).

Corre con flspos. Cualquier argumento de línea de comandos para el modo detallado.

#include <cstdio>
#define _USE_MATH_DEFINES
#undef __STRICT_ANSI__
#include <cmath>
#include <cfloat>
#include <vector>
#include <iterator>
#include <algorithm>

using namespace std;

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

template <typename T> struct widen;
template <> struct widen<int> { typedef long long type; };

template <typename T>
inline typename widen<T>::type mul(T x, T y) {
    return typename widen<T>::type(x) * typename widen<T>::type(y);
}
template <typename T> inline T div_ceil(T a, T b) { return (a + b - 1) / b; }
template <typename T> inline typename widen<T>::type sq(T x) { return mul(x, x); }
template <typename T>
T gcd(T a, T b) { while (b) { T t = a; a = b; b = t % b; } return a; }
template <typename T>
inline typename widen<T>::type lcm(T a, T b) { return mul(a, b) / gcd(a, b); }
template <typename T>
T div_mod_n(T a, T b, T n) {
    if (b == 0) return a == 0 ? 0 : -1;
    const T n_over_b = n / b, n_mod_b = n % b;
    for (T m = 0; m < n; m += n_over_b + 1) {
        if (a % b == 0) return m + a / b;
        a -= b - n_mod_b;
        if (a < 0) a += n;
    }
    return -1;
}

template <typename T> struct pythagorean_triplet { T a, b, c; };
template <typename T>
struct pythagorean_triplet_generator {
    typedef pythagorean_triplet<T> result_type;
private:
    typedef typename widen<T>::type WT;
    result_type p_triplet;
    WT p_c2b2;
public:
    pythagorean_triplet_generator(const result_type& triplet = {3, 4, 5}) :
        p_triplet(triplet), p_c2b2(sq(triplet.c) - sq(triplet.b))
    {}
    const result_type& operator*() const { return p_triplet; }
    const result_type* operator->() const { return &p_triplet; }
    pythagorean_triplet_generator& operator++() {
        do {
            if (++p_triplet.b == p_triplet.c) {
                ++p_triplet.c;
                p_triplet.b = ceil(p_triplet.c * M_SQRT1_2);
                p_c2b2 = sq(p_triplet.c) - sq(p_triplet.b);
            } else
                p_c2b2 -= 2 * p_triplet.b - 1;
            p_triplet.a = sqrt(p_c2b2);
        } while (sq(p_triplet.a) != p_c2b2 || gcd(p_triplet.b, p_triplet.a) != 1);
        return *this;
    }
    result_type operator()() { result_type t = **this; ++*this; return t; }
};

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

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

    const size_t small_triplet_count = 1000;
    vector<pythagorean_triplet<int>> small_triplets;
    small_triplets.reserve(small_triplet_count);
    generate_n(
        back_inserter(small_triplets),
        small_triplet_count,
        pythagorean_triplet_generator<int>()
    );

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

    int output_counter = 0;
    for (int x = min; x <= max; ++x) add(x, 0,    x, 1);
    for (pythagorean_triplet_generator<int> i; i->c <= max; ++i) {
        const auto& t1 = *i;

        for (int n = div_ceil(min, t1.c); n <= max / t1.c; ++n)
            add(n * t1.b, n * t1.a,    n * t1.c, 1);

        auto find_false_positives = [&] (int r, int x, int y) {
            {
                int n = div_ceil(min - r, t1.c);
                int min_r = r + n * t1.c;
                int max_n = n + (max - min_r) / t1.c;
                for (; n <= max_n; ++n)
                    add(r + n * t1.c, 0,    x + n * t1.b, y + n * t1.a);
            }
            for (const auto t2 : small_triplets) {
                int m = div_mod_n((t2.c - r % t2.c) % t2.c, t1.c % t2.c, t2.c);
                if (m < 0) continue;
                int sr = r + m * t1.c;
                int c = lcm(t1.c, t2.c);
                int min_n = div_ceil(min - sr, c);
                int min_r = sr + min_n * c;
                if (min_r > max) continue;
                int x1 = x + m * t1.b, y1 = y + m * t1.a;
                int x2 = t2.b * (sr / t2.c), y2 = t2.a * (sr / t2.c);
                int a1 = t1.a * (c / t1.c), b1 = t1.b * (c / t1.c);
                int a2 = t2.a * (c / t2.c), b2 = t2.b * (c / t2.c);
                int max_n = min_n + (max - min_r) / c;
                int max_r = sr + max_n * c;
                for (int n = min_n; n <= max_n; ++n) {
                    add(
                        x2 + n * b2, y2 + n * a2,
                        x1 + n * b1, y1 + n * a1
                    );
                }
            }
        };
        {
            int m = div_mod_n((t1.a - t1.c % t1.a) % t1.a, t1.b % t1.a, t1.a);
            find_false_positives(
                /* r = */ (mul(m, t1.c) + t1.b) / t1.a,
                /* x = */ (mul(m, t1.b) + t1.c) / t1.a,
                /* y = */ m
            );
        } {
            int m = div_mod_n((t1.b - t1.c % t1.b) % t1.b, t1.a, t1.b);
            find_false_positives(
                /* r = */ (mul(m, t1.c) + t1.a) / t1.b,
                /* x = */ m,
                /* y = */ (mul(m, t1.a) + t1.c) / t1.b
            );
        }

        if (output_counter++ % 50 == 0)
            printf("%d\n", found), fflush(stdout);
    }
    printf("%d\n", found);
}

Pitón, 27,797,402

Solo para poner el listón un poco más alto ...

from sys import argv
verbose = len(argv) > 1
found = 0
for x in xrange(67108864, 94906266):
    found += 1
    if verbose:
        print "(%d, 0) (%d, 1)" % (x, x)
print found

Es fácil verificar que para todos los 67,108,864 <= x <= 94,906,265 = floor (sqrt (2 ⁵³ )) los pares (x, 0) y (x, 1) son falsos positivos.

Por qué funciona : 67,108,864 = 2 ²⁶ . Por lo tanto, todos los números x en el rango anterior tienen la forma 2 ²⁶ + x ' para algunos 0 <= x' <2 ²⁶ . Para todos los positivos e , (x + e) ² = x ² + 2xe + e ² = x ² + 2 ²⁷ e + 2x'e + e ² . Si queremos tener
(x + e) ² = x ² + 1 necesitamos al menos 2 ²⁷ e <= 1 , es decir, e <= 2 ^-27 Sin embargo, dado que la mantisa de los números de coma flotante de doble precisión es de 52 bits de ancho, la más pequeña e tal que x + e> x es e = 2 ^{26 - 52} = 2 ^-26 . En otras palabras, el número representable más pequeño mayor que x es x + 2 ^-26, mientras que el resultado de sqrt (x ² + 1) es como máximo x + 2 ^-27 . Dado que el modo de redondeo IEEE-754 predeterminado es redondear al más cercano; lazos a pares, siempre se redondeará a x y nunca a x + 2 ^-26 (donde el desempate solo es relevante para x = 67,108,864, como mucho. Cualquier número mayor se redondeará a x independientemente).

C ++, 75,000,000+

Recuerde que 3 ² + 4 ² = 5 ² . Lo que esto significa es que el punto (4, 3) se encuentra en el círculo de radio 5 centrado alrededor del origen. En realidad, para todos los enteros n , (4n, 3n) se encuentra en dicho círculo de radio 5n . Para n suficientemente grande (es decir, tal que 5n> = 2 ²⁶ ), ya conocemos un falso positivo para todos los puntos en este círculo: (5n, 1) . ¡Excelente! ¡Eso es otros 27,797,402 / 5 pares de falsos positivos gratis! ¿Pero por qué parar aquí? (3, 4, 5) no es el único triplete.

Este programa busca todos los tripletes enteros positivos (a, b, c) de modo que a ² + b ² = c ² , y cuenta los falsos positivos de esta manera. Llega a 70,000,000 falsos positivos bastante rápido, pero luego disminuye considerablemente a medida que aumentan los números.

Compilar con g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Para verificar los resultados, agregue -DVERIFY(esto será notablemente más lento).

Corre con flspos. Cualquier argumento de línea de comandos para el modo detallado.

#include <cstdio>
#include <cmath>
#include <cfloat>

using namespace std;

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

template <typename T> inline long long sqr(T x) { return 1ll * x * x; }
template <typename T>
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 > 1;

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

    int found = 0;
    auto add = [=, &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, "Wrong false-positive: (%d, %d) (%d, %d)\n",
                    x1, y1, x2, y2);
            return;
        }
#endif
        if (verbose) printf("(%d, %d) (%d, %d)\n", x1, x2, y1, y2);
        ++found;
    };

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

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

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

    printf("%d\n", found);
}

C ++ 11 - 100,993,667

EDITAR: Nuevo programa.

El viejo usaba demasiada memoria. Este reduce a la mitad el uso de memoria mediante el uso de una matriz de vectores gigantes en lugar de una tabla hash. También elimina el hilo aleatorio cruft.

   /* by feersum  2014/9
   http://codegolf.stackexchange.com/questions/37627/false-positives-on-an-integer-lattice */
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <vector>
using namespace std;
#define ul unsigned long long

#define K const



#define INS(A)   { bool already = false; \
    for(auto e = res[A.p[0][0]].end(), it = res[A.p[0][0]].begin(); it != e; ++it) \
        if(A.p[0][1] == it->y1 && A.p[1][0] == it->x2 && A.p[1][1] == it->y2) { \
            already = true; \
            break; } \
    if(!already) res[A.p[0][0]].push_back( {A.p[0][1], A.p[1][0], A.p[1][1]} ), ++n; }

#define XMAXMIN (1<<26)

struct ints3 {
    int y1, x2, y2;
};


struct pparm {
    int a,b,c,d;
    int E[4];
    pparm(int A,int B,int C, int D):
        E{B*B+D*D,A*B+C*D,A*A+C*C+2*(B+D),A+C}
    {
        a=A;b=B;c=C;d=D;
    }

};

struct ans {
    int p[2][2];

};
ostream&operator<<(ostream&o, ans&a)
{
    o<<'('<<a.p[0][0]<<','<<a.p[0][1]<<"),("<<a.p[1][0]<<','<<a.p[1][1]<<')'<<endl;
    return o;
}



vector<ints3> res[XMAXMIN];

bool print;
int n;

void gen(K pparm&p1, K pparm&p2)
{
#ifdef DBUG
    for(int i=0;i<2;i++){
    K pparm&p=i?p2:p1;
    cout<<' '<<p.a<<' '<<p.b<<' '<<p.c<<' '<<p.d<<' ';}
    cout<<endl;
#endif

    for(ul x = 0; ; ++x) {
        ans a;
        ul s[2];
        for(int i = 0; i < 2; i++) {
            K pparm &p = i?p2:p1;
            int *pt = a.p[i];
            pt[0] = p.b+x*(p.a+x);
            pt[1] = p.d+x*(p.c+x);
            s[i] = (ul)pt[0]*pt[0] + (ul)pt[1]*pt[1];
        }
        if(*s >> 53)
            break;
if(s[1] - s[0] != 1)
exit(4);

        if(sqrt(s[0]) == sqrt(s[1])) {
             for(int i = 0; i < 2; i++)
                if(a.p[i][0] > a.p[i][1])
                    swap(a.p[i][0], a.p[i][1]);
            if(a.p[0][0] > a.p[0][1])
                for(int i = 0; i < 2; i++)
                    swap(a.p[0][i], a.p[1][i]);
            INS(a)
        }
    }
}



int main(int ac, char**av)
{
    for(int i = 1; i < ac; i++) {
        print |= !strcmp(av[1], "-P");
    }


    #define JMAX 43000000
    for(ul j = 0; j < JMAX; j++) {
        pparm p1(-~j,j,~-j,0),p2(j,1,j,j);
        gen(p1,p2);
        if(!print && !(j%1024))
#ifdef DBUG
            cout<<j<<' ',
#endif
            cout<<n<<endl;

    }
    if(print) 
        for(vector<ints3>& v: res)
            for(ints3& i: v)
                printf("(%d,%d),(%d,%d)\n", &v - res, i.y1, i.x2, i.y2);

    return 0;
}

Ejecute con un -Pargumento para imprimir los puntos en lugar del número de los mismos.

Para mí, lleva menos de 2 minutos en el modo de conteo y aproximadamente 5 minutos con la impresión dirigida a un archivo (~ 4 GB), por lo que no se limitó a E / S.

Mi programa original estaba ordenado, pero dejé caer la mayor parte ya que solo podía producir en el orden de 10 ^ 5 resultados. Lo que hizo fue buscar parametrizaciones de la forma (x ^ 2 + Ax + B, x ^ 2 + Cx + D), (x ^ 2 + ax + b, x ^ 2 + cx + d) de modo que para cualquier x, (x ^ 2 + Ax + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Cuando encontró tal conjunto de parámetros {a, b, c, d, A, B, C, D}, procedió a verificar todos los valores de x por debajo del máximo. Mientras miraba mi salida de depuración de este programa, noté una cierta parametrización de la parametrización de la parametrización que me permitió producir muchos números fácilmente. Elegí no imprimir los números de Ell ya que tenía muchos. Esperemos que ahora alguien no imprima nuestros dos conjuntos de números y afirme ser el ganador :)

 /* by feersum  2014/9
   http://codegolf.stackexchange.com/questions/37627/false-positives-on-an-integer-lattice */
    #include <iostream>
    #include <cmath>
    #include <cstdlib>
    #include <cstring>
    #include <functional>
    #include <unordered_set>
    #include <thread>
using namespace std;
#define ul unsigned long long

#define h(S) unordered_##S##set
#define P 2977953206964783763LL
#define K const

#define EQ(T, F)bool operator==(K T&o)K{return!memcmp(F,o.F,sizeof(F));}

struct pparm {
    int a,b,c,d;
    int E[4];
    pparm(int A,int B,int C, int D):
        E{B*B+D*D,A*B+C*D,A*A+C*C+2*(B+D),A+C}
    {
        a=A;b=B;c=C;d=D;
    }
    EQ(pparm,E)
};

struct ans {
    int p[2][2];
    EQ(ans,p)
};
ostream&operator<<(ostream&o, ans&a)
{
    o<<'('<<a.p[0][0]<<','<<a.p[0][1]<<"),("<<a.p[1][0]<<','<<a.p[1][1]<<')'<<endl;
    return o;
}

#define HASH(N,T,F) \
struct N { \
    size_t operator() (K T&p) K { \
        size_t h = 0; \
        for(int i = 4; i--; ) \
            h=h*P+((int*)p.F)[i]; \
        return h; \
    }};

#define INS(r, a) { \
    bool new1 = r.insert(a).second; \
    n += new1; \
    if(print && new1) \
        cout<<a; }

HASH(HA,ans,p)

bool print;
int n;

void gen(h()<ans,HA>&r, K pparm&p1, K pparm&p2)
{
#ifdef DBUG
    for(int i=0;i<2;i++){
    K pparm&p=i?p2:p1;
    cout<<' '<<p.a<<' '<<p.b<<' '<<p.c<<' '<<p.d<<' ';}
    cout<<endl;
#endif

    for(ul x = 0; ; ++x) {
        ans a;
        ul s[2];
        for(int i = 0; i < 2; i++) {
            K pparm &p = i?p2:p1;
            int *pt = a.p[i];
            pt[0] = p.b+x*(p.a+x);
            pt[1] = p.d+x*(p.c+x);
            s[i] = (ul)pt[0]*pt[0] + (ul)pt[1]*pt[1];
        }
        if(*s >> 53)
            break;
if(s[1] - s[0] != 1)
exit(4);

        if(sqrt(s[0]) == sqrt(s[1])) {
             for(int i = 0; i < 2; i++)
                if(a.p[i][0] > a.p[i][1])
                    swap(a.p[i][0], a.p[i][1]);
            INS(r,a)
        }
    }
    //if(!print) cout<<n<<endl;
}

void endit()
{
    this_thread::sleep_for(chrono::seconds(599));
    exit(0);
}

int main(int ac, char**av)
{
    bool kill = false;
    for(int i = 1; i < ac; i++) {
        print |= ac>1 && !stricmp(av[1], "-P");
        kill |= !stricmp(av[i], "-K");
    }

    thread KILLER;
    if(kill)
        KILLER = thread(endit);

    h()<ans, HA> res;
    res.reserve(1<<27);

    #define JMAX 43000000
    for(ul j = 0; j < JMAX; j++) {
        pparm p1(-~j,j,~-j,0),p2(j,1,j,j);
        gen(res,p1,p2);
        if(!print && !(j%1024))
#ifdef DBUG
            cout<<j<<' ',
#endif
            cout<<n<<endl;

    }
    exit(0);
}

Escaneo circular Java, Bresenham-esque

Heurísticamente espero obtener más colisiones comenzando en el extremo más ancho del anillo. Esperaba obtener alguna mejora haciendo un escaneo para cada colisión, registrando valores para los que surplusestá entre 0e r2max - r2inclusive, pero en mis pruebas resultó ser más lento que esta versión. Del mismo modo, intenta utilizar un único int[]búfer en lugar de mucha creación de matrices y listas de dos elementos. La optimización del rendimiento es realmente una bestia extraña.

Ejecute con un argumento de línea de comandos para la salida de los pares, y sin recuentos simples.

import java.util.*;

public class CodeGolf37627 {
    public static void main(String[] args) {
        final int M = 144;
        boolean[] possible = new boolean[M];
        for (int i = 0; i <= M/2; i++) {
            for (int j = 0; j <= M/2; j++) {
                possible[(i*i+j*j)%M] = true;
            }
        }

        long count = 0;
        double sqrt = 0;
        long r2max = 0;
        List<int[]> previousPoints = null;
        for (long r2 = 1L << 53; ; r2--) {
            if (!possible[(int)(r2 % M)]) continue;

            double r = Math.sqrt(r2);
            if (r != sqrt) {
                sqrt = r;
                r2max = r2;
                previousPoints = null;
            }
            else {
                if (previousPoints == null) previousPoints = findLatticePointsBresenham(r2max, (int)r);

                if (previousPoints.size() == 0) {
                    r2max = r2;
                    previousPoints = null;
                }
                else {
                    List<int[]> points = findLatticePointsBresenham(r2, (int)r);
                    for (int[] p1 : points) {
                        for (int[] p2 : previousPoints) {
                            if (args.length > 0) System.out.format("(%d, %d) (%d, %d)\n", p1[0], p1[1], p2[0], p2[1]);
                            count++;
                        }
                    }
                    previousPoints.addAll(points);
                    System.out.println(count);
                }
            }
        }
    }

    // Surprisingly, this seems to be faster than doing one scan for all two or three r2s.
    private static List<int[]> findLatticePointsBresenham(long r2, long r) {
        List<int[]> rv = new ArrayList<int[]>();
        // Require 0 = y = x
        long x = r, y = 0, surplus = r2 - r * r;
        while (y <= x) {
            if (surplus == 0) rv.add(new int[]{(int)x, (int)y});

            // Invariant: surplus = r2 - x*x - y*y >= 0
            y++;
            surplus -= 2*y - 1;
            if (surplus < 0) {
                x--;
                surplus += 2*x + 1;
            }
        }

        return rv;
    }
}

Java - 27,817,255

La mayoría de estos son los mismos que muestra Ell , y el resto se basa en (j,0) (k,l). Para cada uno j, retrocedo algunos cuadrados y verifico si el resto da un falso positivo. Esto ocupa básicamente todo el tiempo con solo una ganancia de 25k (aproximadamente 0.1%) sobre solo (j,0) (j,1), pero una ganancia es una ganancia.

Esto terminará en menos de diez minutos en mi máquina, pero no sé lo que tienes. Por razones, si no termina antes de que acabe el tiempo, tendrá un puntaje drásticamente peor. En ese caso, puede ajustar el divisor en la línea 8 para que termine en el tiempo (esto simplemente determina qué tan atrás camina para cada uno j). Para algunos divisores diferentes, los puntajes son:

11    27817255 (best on OPs machine)
10    27818200
8     27820719
7     27822419 (best on my machine)

Para activar la salida para cada partido (y, dios, es lento si lo hace), solo descomente las líneas 10 y 19.

public class FalsePositive {
    public static void main(String[] args){
        long j = 67108864;
        long start = System.currentTimeMillis();
        long matches=0;
        while(j < 94906265 && System.currentTimeMillis()-start < 599900){
            long jSq = j*j;
            long limit = (long)Math.sqrt(j)/11; // <- tweak to fit inside 10 minutes for best results
            matches++; // count an automatic one for (j,0)(j,1)
            //System.out.println("("+j+",0) ("+j+",1)");        
            for(int i=1;i<limit;i++){
                long k = j-i;
                long kSq = k*k;
                long l = (long)Math.sqrt(jSq-kSq);
                long lSq = l*l;
                if(kSq+lSq != jSq){
                    if(Math.sqrt(kSq+lSq)==Math.sqrt(jSq)){
                        matches++;
                        //System.out.println("("+j+",0) ("+k+","+l+")");        
                    }
                }
            }
            j++;
        }
        System.out.println("\n"+matches+" Total matches, got to j="+j);
    }
}

Como referencia, las primeras 20 salidas que da (para divisor = 7, excluyendo (j,0)(j,1)tipos) son:

(67110083,0) (67109538,270462)
(67110675,0) (67109990,303218)
(67111251,0) (67110710,269470)
(67111569,0) (67110668,347756)
(67112019,0) (67111274,316222)
(67112787,0) (67111762,370918)
(67115571,0) (67115518,84346)
(67117699,0) (67117698,11586)
(67117971,0) (67117958,41774)
(67120545,0) (67120040,260368)
(67121043,0) (67120118,352382)
(67122345,0) (67122320,57932)
(67122449,0) (67122444,25908)
(67122633,0) (67122328,202348)
(67122729,0) (67121972,318784)
(67122849,0) (67122568,194224)
(67124195,0) (67123818,224970)
(67125201,0) (67125172,62396)
(67125705,0) (67124632,379540)
(67126195,0) (67125882,204990)

Julia, 530 falsos positivos

Aquí hay una búsqueda de fuerza bruta muy ingenua, que puede ver como una implementación de referencia.

num = 0
for i = 60000000:-1:0
    for j = i:-1:ifloor(0.99*i)
        s = i*i + j*j
        for x = ifloor(sqrt(s/2)):ifloor(sqrt(s))
            min_y = ifloor(sqrt(s - x*x))
            max_y = min_y+1
            for y = min_y:max_y
                r = x*x + y*y
                if r != s && sqrt(r) == sqrt(s)
                    num += 1
                    if num % 10 == 0
                        println("Found $num pairs")
                    end
                    #@printf("(i,j) = (%d,%d); (x,y) = (%d,%d); s = %d, r = %d\n", i,j,x,y,s,r)
                end
            end
        end
    end
end

Puede imprimir los pares (y sus magnitudes cuadradas exactas) descomentando la @printflínea.

Básicamente, esto inicia la búsqueda x = y = 6e7del primer par de coordenadas y escanea aproximadamente el 1% del camino hasta el eje x antes de disminuir x. Luego, para cada par de coordenadas, verifica el arco completo de la misma magnitud (redondeando hacia arriba y hacia abajo) en busca de una colisión.

El código supone que se ejecuta en un sistema de 64 bits, de modo que los tipos enteros y de coma flotante predeterminados son de 64 bits (si no, puede crearlos con int64()y float64()constructores).

Eso arroja unos escasos 530 resultados.