Entre os melhores

 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

fundo

Ao trabalhar em uma grade bidimensional de coordenadas inteiras, às vezes você deseja saber se dois vetores (com componentes inteiros) têm a mesma magnitude. Obviamente, na geometria euclidiana, a magnitude de um vetor (x,y)é dada por

√(x² + y²)

Portanto, uma implementação ingênua pode calcular esse valor para os dois vetores e comparar os resultados. Isso não apenas gera um cálculo desnecessário de raiz quadrada, mas também causa problemas com imprecisões de ponto flutuante, que podem gerar falsos positivos: vetores cujas magnitudes são diferentes, mas onde os dígitos significativos na representação de ponto flutuante são todos idênticos.

Para os propósitos deste desafio, definimos um falso positivo como um par de pares de coordenadas (a,b)e (c,d)para o qual:

• Sua magnitude ao quadrado é diferente quando representada como números inteiros não assinados de 64 bits.
• Sua magnitude é idêntica quando representada como número de ponto flutuante binário de 64 bits e calculada através de uma raiz quadrada de 64 bits (conforme IEEE 754 ).

Como um exemplo, usando representações de 16 bits em vez de (64), o mais pequeno ^um par de vectores que produz um falso positivo é

(25,20) and (32,0)

Suas magnitudes ao quadrado ao quadrado são 1025e 1024. Tomando os rendimentos da raiz quadrada

32.01562118716424 and 32.0

Mas em carros alegóricos de 16 bits, ambos são truncados 32.0.

Da mesma forma, o menor par ^{2 que} produz um falso positivo para representações de 32 bits é

(1659,1220) and (1951,659)

_{¹ "Menor", medido pela magnitude do ponto flutuante de 16 bits.
² "Menor", medido pela magnitude do ponto flutuante de 32 bits.}

Por fim, aqui estão alguns 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) *

* O último caso é um dos vários com a menor magnitude possível para falsos positivos de 64 bits.

O desafio

Em menos de 10.000 bytes de código, usando um único encadeamento, você encontrará o número de falsos positivos para números de ponto flutuante (binário) de 64 bits no intervalo de coordenadas 0 ≤ y ≤ x(ou seja, apenas dentro do primeiro octante do plano euclidiano) tal que dentro de 10 minutos . Se dois envios empatam para o mesmo número de pares, o desempate é o tempo real necessário para encontrar o último desses pares.x² + y² ≤ 253

Seu programa não deve usar mais de 4 GB de memória a qualquer momento (por motivos práticos).

Deve ser possível executar o seu programa de dois modos: um que gera todos os pares como o encontra e outro que gera apenas o número de pares encontrados no final. O primeiro será usado para verificar a validade de seus pares (observando algumas amostras de saídas) e o segundo será usado para realmente cronometrar o envio. Observe que a impressão deve ser a única diferença. Em particular, o programa de contagem pode não codificar o número de pares que poderia encontrar. Ele ainda deve executar o mesmo loop que seria usado para imprimir todos os números e omitir apenas a impressão em si!

Testarei todos os envios no meu laptop com Windows 8; portanto, pergunte nos comentários se você deseja usar um idioma não muito comum.

Observe que os pares não devem ser contados duas vezes sob a comutação do primeiro e do segundo par de coordenadas.

Observe também que executarei seu processo por meio de um controlador Ruby, que matará seu processo se não terminar após 10 minutos. Certifique-se de gerar o número de pares encontrados até então. Você pode acompanhar o tempo e imprimir o resultado imediatamente antes dos 10 minutos decorridos, ou você apenas gera o número de pares encontrados esporadicamente, e tomarei o último número como sua pontuação.

C ++, 275.000.000+

Vamos nos referir a pares cuja magnitude é exatamente representável, como (x, 0) , como pares honestos e a todos os outros pares como pares desonestos de magnitude m , onde m é a magnitude relatada incorretamente do par. O primeiro programa do post anterior usou um conjunto de pares estreitamente relacionados de pares honestos e desonestos:
(x, 0) e (x, 1) , respectivamente, para x suficientemente grandes. O segundo programa usou o mesmo conjunto de pares desonestos, mas estendeu o conjunto de pares honestos, procurando todos os pares honestos de magnitude integral. O programa não termina em dez minutos, mas encontra a grande maioria de seus resultados muito cedo, o que significa que a maior parte do tempo de execução é desperdiçada. Em vez de continuar procurando pares honestos cada vez menos frequentes, esse programa usa o tempo livre para fazer a próxima coisa lógica: estender o conjunto de pares desonestos .

No post anterior, sabemos que para todos os números inteiros grandes o suficiente r , sqrt (r ² + 1) = r , em que sqrt é a função de raiz quadrada de ponto flutuante. Nosso plano de ataque é encontrar pares P = (x, y) tais que x ² + y ² = r ² + 1 para algum número inteiro grande o suficiente r . Isso é simples o suficiente, mas procurar ingenuamente esses pares individuais é muito lento para ser interessante. Queremos encontrar esses pares em massa, assim como fizemos para pares honestos no programa anterior.

Seja { v , w } um par de vetores ortonormal. Para todos os escalares reais r , || r v + w || ² = r ² + 1 . No ℝ ² , este é um resultado direto do teorema de Pitágoras:

Estamos à procura de vetores V e w tal que existe um número inteiro r para o qual x e y também são inteiros. Como observação lateral, observe que o conjunto de pares desonestos que usamos nos dois programas anteriores era simplesmente um caso especial disso, onde { v , w } era a base padrão de ℝ ² ; desta vez, queremos encontrar uma solução mais geral. É aqui que os trigêmeos pitagóricos (trigêmeos inteiros (a, b, c) satisfazem a ² + b ² = c ², que usamos no programa anterior) voltam.

Seja (a, b, c) um trigêmeo pitagórico. Os vectores de v = (b / c, a / c) e w = (-a / C, b / c), (e também
w = (a / c, b / c) ) são ortonormais, como é fácil de verificar . Como se vê, para qualquer escolha do trigêmeo pitagórico, existe um número inteiro r tal que x e y são números inteiros. Para provar isso e encontrar efetivamente r e P , precisamos de um pouco de teoria dos números / grupos; Vou poupar os detalhes. De qualquer maneira, suponha que tenhamos nossa integral r , x e y . Ainda faltam algumas coisas: precisamos de rser grande o suficiente e queremos um método rápido para derivar muitos outros pares semelhantes deste. Felizmente, existe uma maneira simples de fazer isso.

Observe que a projeção de P em v é r v , portanto r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , tudo isso para dizer que xb + ya = rc . Como resultado, para todos os números inteiros n , (x + bn) ² + (y + an) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1. Em outras palavras, a magnitude ao quadrado dos pares da forma
(x + bn, y + an) é (r + cn) ² + 1 , que é exatamente o tipo de pares que estamos procurando! Para n grande o suficiente , esses são pares desonestos de magnitude r + cn .

É sempre bom olhar para um exemplo concreto. Se pegarmos o trigêmeo pitagórico (3, 4, 5) , então em r = 2 , temos P = (1, 2) (você pode verificar que (1, 2) · (4/5, 3/5) = 2 e, claramente, um ² + 2 ² = 2 ² + 1 .) a adição de 5 a r e (4, 3) de P nos leva a r '= 2 + 5 = 7 e P' = (1 + 4, 2 + 3) = (5, 5) . Lo e eis que 5 ² + 5 ² = 7 ² + 1. As próximas coordenadas são r '' = 12 e P '' = (9, 8) e, novamente, 9 ² + 8 ² = 12 ² + 1 , e assim por diante ...

Uma vez que r é grande o suficiente, nós começar a receber desonestas pares com incrementos de magnitude de 5 . São aproximadamente 27.797.402 / 5 pares desonestos.

Então agora temos muitos pares desonestos de magnitude integral. Podemos facilmente associá-los aos pares honestos do primeiro programa para formar falsos positivos e, com o devido cuidado, também podemos usar os pares honestos do segundo programa. Isso é basicamente o que esse programa faz. Como o programa anterior, ele também encontra a maioria dos resultados muito cedo - chega a 200.000.000 de falsos positivos em alguns segundos - e depois diminui consideravelmente.

Compile com g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Para verificar os resultados, adicione -DVERIFY(isso será notavelmente mais lento).

Corra com flspos. Qualquer argumento da linha de comandos para o modo detalhado.

#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);
}

Python, 27.797.402

Só para deixar a barra um pouco mais alta ...

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

É fácil verificar se, para todos os 67.108.864 <= x <= 94.906.265 = floor (sqrt (2 ⁵³ )), os pares (x, 0) e (x, 1) são falsos positivos.

Por que funciona : 67.108.864 = 2 ²⁶ . Portanto, todos os números x no intervalo acima têm a forma 2 ²⁶ + x ' para alguns 0 <= x' <2 ²⁶ . Para todo positivo e , (x + e) ² = x ² + 2xe + e ² = x ² + 2 ²⁷ e + 2x'e + e ² . Se queremos ter
(x + e) ² = x ² + 1 , precisamos de pelo menos 2 ²⁷ e <= 1 , isto é, e <= ^2-27 No entanto, como a mantissa de números de ponto flutuante de precisão dupla tem 52 bits de largura, o menor e é tal que x + e> x é e = 2 ^{26 - 52} = 2 ^-26 . Em outras palavras, o menor número representável maior que x é x + ^2-26, enquanto o resultado de sqrt (x ² + 1) é no máximo x + ^2-27 . Como o modo de arredondamento padrão IEEE-754 é o arredondamento para o mais próximo; empates para pares, sempre arredondará para x e nunca para x + ^2-26 (onde o desempate é realmente relevante apenas para x = 67.108.864, se houver. Qualquer número maior será arredondado para x independentemente).

C ++, 75.000.000+

Lembre-se de que 3 ² + 4 ² = 5 ² . O que isto significa é que o ponto (4, 3) se encontra no círculo do raio 5, centralizado em torno da origem. Na verdade, para todos os números inteiros n , (4n, 3n) se encontra esse círculo de raio 5n . Para n grande o suficiente (ou seja, tal que 5n> = 2 ²⁶ ), já sabemos um falso positivo para todos os pontos deste círculo: (5n, 1) . Ótimo! São outros 27.797.402 / 5 pares de falsos positivos gratuitos aqui! Mas por que parar aqui? (3, 4, 5) não é o único desses trigêmeos.

Este programa procura todos os trigêmeos inteiros positivos (a, b, c) de modo que a ² + b ² = c ² e conte falso-positivos dessa maneira. Chega a 70.000.000 de falsos positivos bastante rápido, mas diminui consideravelmente à medida que os números aumentam.

Compile com g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Para verificar os resultados, adicione -DVERIFY(isso será notavelmente mais lento).

Corra com flspos. Qualquer argumento da linha de comandos para o modo detalhado.

#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

EDIT: Novo programa.

O antigo usava muita memória. Este reduz pela metade o uso de memória usando uma matriz de vetores gigante em vez de uma tabela de hash. Também remove a crosta aleatória do encadeamento.

   /* 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;
}

Execute um -Pargumento para imprimir os pontos em vez do número.

Para mim, leva menos de 2 minutos no modo de contagem e cerca de 5 minutos com a impressão direcionada para um arquivo (~ 4 GB), para que não se torne limitado à E / S.

Meu programa original foi arrumado, mas perdi a maior parte dele, pois só podia produzir na ordem de 10 ^ 5 resultados. O que ele fez foi procurar parametrizações da forma (x ^ 2 + Ax + B, x ^ 2 + Cx + D), (x ^ 2 + ax + b, x ^ 2 + cx + d) para que, para qualquer A soma de dois algarismos distintos é igual a: (x ^ 2 + Ax + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Quando encontrou esse conjunto de parâmetros {a, b, c, d, A, B, C, D}, passou a verificar todos os valores de x abaixo do máximo. Enquanto olhava para minha saída de depuração deste programa, notei uma certa parametrização da parametrização da parametrização que me permitiu produzir muitos números facilmente. Optei por não imprimir os números de Ell, já que eu tinha muitos. Espero que agora alguém não imprima ambos os nossos conjuntos de números e afirme ser o vencedor :)

 /* 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);
}

Varredura de círculo em Java, estilo Bresenham

Heuristicamente, espero obter mais colisões começando no final mais amplo do espaço anular. Eu esperava obter alguma melhoria fazendo uma varredura para cada colisão, registrando valores surplusentre 0e r2max - r2inclusivos, mas em meus testes isso se mostrou mais lento que esta versão. Da mesma forma, tenta usar um único int[]buffer em vez de criar muitas listas e matrizes de dois elementos. A otimização de desempenho é realmente uma besta estranha.

Execute com um argumento de linha de comando para saída dos pares e sem contagens 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

A maioria é igual ao que Ell mostra , e o restante é baseado (j,0) (k,l). Para cada uma j, eu recuo alguns quadrados e verifico se o restante dá um falso positivo. Isso ocupa basicamente o tempo todo com apenas um ganho de 25k (cerca de 0,1%) (j,0) (j,1), mas um ganho é um ganho.

Isso terminará em menos de dez minutos na minha máquina, mas não sei o que você tem. Por razões, se não terminar antes que o tempo se esgote, terá uma pontuação drasticamente pior. Nesse caso, você pode ajustar o divisor na linha 8 para que ele termine no tempo (isso simplesmente determina a que distância ele caminha para cada um j). Para alguns vários divisores, as pontuações são:

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

Para ativar a saída para cada correspondência (e, Deus, é lento se você o fizer), apenas descomente as linhas 10 e 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);
    }
}

Para referência, as primeiras 20 saídas fornecidas (para o divisor = 7, excluindo (j,0)(j,1)tipos) são:

(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

Aqui está uma pesquisa de força bruta muito ingênua, que você pode ver como uma implementação de referência.

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

Você pode imprimir os pares (e suas magnitudes quadradas exatas) descomentando a @printflinha.

Basicamente, isso inicia a busca no x = y = 6e7primeiro par de coordenadas e escaneia cerca de 1% do caminho até o eixo x antes de decrementar x. Em seguida, para cada par de coordenadas, ele verifica o arco inteiro da mesma magnitude (arredondando para cima e para baixo) em busca de uma colisão.

O código assume que ele é executado em um sistema de 64 bits, para que os tipos inteiros padrão e de ponto flutuante são aquelas de 64 bits (se não, você pode criá-los com int64()e float64()construtores).

Isso produz apenas 530 resultados.