Tabela liderów

 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

tło

Pracując na siatce 2D współrzędnych całkowitych, czasami chcesz wiedzieć, czy dwa wektory (ze składowymi liczb całkowitych) mają tę samą wielkość. Oczywiście w geometrii euklidesowej wielkość wektora (x,y)podaje

√(x² + y²)

Zatem naiwna implementacja może obliczyć tę wartość dla obu wektorów i porównać wyniki. Nie tylko powoduje to niepotrzebne obliczanie pierwiastka kwadratowego, ale także powoduje problemy z niedokładnościami zmiennoprzecinkowymi, które mogą dawać fałszywie dodatnie: wektory o różnej wielkości, ale w których wszystkie cyfry znaczące w reprezentacji zmiennoprzecinkowej są identyczne.

Na potrzeby tego wyzwania definiujemy wartość fałszywie dodatnią jako parę par współrzędnych (a,b)i (c,d)dla których:

• Ich kwadratowa wielkość jest różna, gdy jest reprezentowana jako 64-bitowe liczby całkowite bez znaku.
• Ich wielkość jest identyczna, gdy jest reprezentowana jako 64-bitowa binarna liczba zmiennoprzecinkowa i obliczana na podstawie 64-bitowego pierwiastka kwadratowego (zgodnie z IEEE 754 ).

Na przykład, używając reprezentacji 16-bitowych (zamiast 64), najmniejsza ¹ para wektorów, która daje wynik fałszywie dodatni, to

(25,20) and (32,0)

Ich kwadratowe wielkości kwadratowe są 1025i 1024. Biorąc pierwiastek kwadratowy daje

32.01562118716424 and 32.0

Ale w pływakach 16-bitowych oba te zostają obcięte 32.0.

Podobnie, najmniejsza ² para dająca fałszywie dodatni dla reprezentacji 32-bitowych to

(1659,1220) and (1951,659)

_{¹ „Najmniejszy” mierzony za pomocą 16-bitowej wielkości zmiennoprzecinkowej.
² „Najmniejszy” mierzony ich 32-bitową wielkością zmiennoprzecinkową.}

Na koniec kilka garści prawidłowych 64-bitowych przypadków:

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

* Ostatni przypadek jest jednym z kilku o najmniejszej możliwej wielkości dla 64-bitowych wyników fałszywie dodatnich.

Wyzwanie

W mniej niż 10 000 bajtów kodu, przy użyciu jednego wątku, znajdziesz tyle fałszywych wyników dodatnich dla liczb zmiennoprzecinkowych 64-bitowych (binarnych) w zakresie współrzędnych 0 ≤ y ≤ x(to znaczy tylko w obrębie pierwszej oktany płaszczyzny euklidesowej) tak, że w ciągu 10 minut . Jeśli dwa zgłoszenia wiążą się dla tej samej liczby par, przerywnikiem remisu jest faktyczny czas potrzebny na znalezienie ostatniej z tych par.x² + y² ≤ 253

Twój program nie może zużywać więcej niż 4 GB pamięci (ze względów praktycznych).

Musi istnieć możliwość uruchomienia programu w dwóch trybach: jeden, który wypisuje każdą znalezioną parę, i jeden, który wypisuje tylko liczbę znalezionych par na końcu. Pierwszy zostanie użyty do zweryfikowania poprawności twoich par (poprzez sprawdzenie próbki wyników), a drugi zostanie wykorzystany do faktycznego pomiaru czasu przesłania. Pamiętaj, że druk musi być jedyną różnicą. W szczególności program zliczający może nie zakodować na stałe liczby par, które może znaleźć. Musi nadal wykonywać tę samą pętlę, która byłaby użyta do wydrukowania wszystkich liczb, i pomijać sam wydruk!

Przetestuję wszystkie zgłoszenia na moim laptopie z systemem Windows 8, więc w komentarzach zapytaj, czy chcesz użyć niezbyt popularnego języka.

Należy pamiętać, że par nie można liczyć dwukrotnie podczas przełączania pierwszej i drugiej pary współrzędnych.

Zauważ też, że uruchomię twój proces za pośrednictwem kontrolera Ruby, który zabije twój proces, jeśli nie zakończy się po 10 minutach. Pamiętaj, aby podać liczbę znalezionych par do tego czasu. Możesz albo samodzielnie śledzić czas i wydrukować wynik tuż przed upływem 10 minut, albo po prostu wypisywać sporadycznie liczbę znalezionych par, a ja wezmę ostatnią taką liczbę jako twój wynik.

C ++, 275 000 000+

Będziemy odnosić się do par, których wielkość jest dokładnie reprezentowalna, takich jak (x, 0) , jako uczciwe pary, a do wszystkich innych par jako nieuczciwe pary wielkości m , gdzie m jest błędnie zgłoszoną wielkością pary. Pierwszy program w poprzednim poście wykorzystywał zestaw ściśle ze sobą powiązanych par uczciwych i nieuczciwych par: odpowiednio
(x, 0) i (x, 1) dla wystarczająco dużego x. Drugi program wykorzystywał ten sam zestaw nieuczciwych par, ale rozszerzył zestaw uczciwych par, szukając wszystkich uczciwych par integralnej wielkości. Program nie kończy się w ciągu dziesięciu minut, ale większość swoich wyników znajduje bardzo wcześnie, co oznacza, że ​​większość czasu działania marnuje się. Zamiast szukać coraz rzadziej uczciwych par, ten program wykorzystuje wolny czas, aby zrobić następną logiczną rzecz: rozszerzenie zestawu nieuczciwych par.

Z poprzedniego postu wiemy, że dla wszystkich wystarczająco dużych liczb całkowitych r , sqrt (r ² + 1) = r , gdzie sqrt jest zmiennoprzecinkową funkcją pierwiastka zmiennoprzecinkowego. Nasz plan ataku polega na znalezieniu par P = (x, y) takich, że x ² + y ² = r ² + 1 dla pewnej dużej liczby całkowitej r . To dość proste, ale naiwne szukanie pojedynczych takich par jest zbyt wolne, aby było interesujące. Chcemy znaleźć te pary masowo, tak jak robiliśmy to dla uczciwych par w poprzednim programie.

Niech { v , w } będzie parą wektorów ortonormalnych. Dla wszystkich rzeczywistych skalarów r , || r v + w || ² = r ² + 1 . W ℝ ² jest to bezpośredni wynik twierdzenia Pitagorasa:

Szukamy wektorów v i w takich, że istnieje liczba całkowita r, dla której x i y są również liczbami całkowitymi. Na marginesie, zauważmy, że zestaw nieuczciwych par, których użyliśmy w poprzednich dwóch programach, był po prostu szczególnym przypadkiem tego, gdzie { v , w } było standardową podstawą ℝ ² ; tym razem chcemy znaleźć bardziej ogólne rozwiązanie. To jest, gdy pitagorejską tryplety (liczby całkowite tryplety (a, b, c) , spełniających w ² + b ² = C ², które wykorzystaliśmy w poprzednim programie) powracają.

Niech (a, b, c) będzie triolą pitagorejską. Wektory v = (b / c, a / c) i w = (-a / c, b / c) (a także
w = (a / c, -b / c) ) są ortonormalne, co łatwo zweryfikować . Jak się okazało, dla każdego wyboru Pitagorasa triplet, istnieje całkowita R w taki sposób, x i y są liczbami całkowitymi. Aby to udowodnić i skutecznie znaleźć r i P , potrzebujemy małej teorii liczb / grup; Oszczędzę szczegóły. Tak czy inaczej, załóżmy, że mamy naszą całkę r , x i y . Wciąż brakuje nam kilku rzeczy: potrzebujemy rbyć wystarczająco duże i chcemy szybkiej metody, aby uzyskać z niej wiele podobnych par. Na szczęście istnieje prosty sposób na osiągnięcie tego.

Zauważ, że rzut P na v jest r v , stąd r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , wszystko to oznacza, że xb + ya = rc . W rezultacie dla wszystkich liczb całkowitych n , (x + bn) ² + (y + an) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1. Innymi słowy, kwadratowa wielkość par postaci
(x + bn, y + an) wynosi (r + cn) ² + 1 , a dokładnie takich par szukamy! Dla wystarczająco dużego n są to nieuczciwe pary wielkości r + cn .

Zawsze miło jest spojrzeć na konkretny przykład. Jeśli weźmiemy triplet Pitagorasa (3, 4, 5) , to przy r = 2 mamy P = (1, 2) (możesz sprawdzić, że (1, 2) · (4/5, 3/5) = 2 i oczywiście 1 ² + 2 ² = 2 ² + 1 ). Dodanie 5 do r oraz (4, 3) do P prowadzi nas do r '= 2 + 5 = 7 i P' = (1 + 4, 2 + 3) = (5, 5) . Lo i oto 5 ² + 5 ² = 7 ² + 1. Kolejne współrzędne to r '' = 12 i P '' = (9, 8) , i znowu 9 ² + 8 ² = 12 ² + 1 itd. I tak dalej ...

Gdy r jest wystarczająco duże, zaczynamy otrzymywać nieuczciwe pary o przyrostach wielkości 5 . To mniej więcej 27 797 402/5 nieuczciwych par.

Więc teraz mamy wiele nieuczciwych par integralnej wielkości. Możemy z łatwością połączyć je z uczciwymi parami pierwszego programu, aby utworzyć fałszywie dodatnie, a przy należytej staranności możemy również użyć uczciwych par drugiego programu. Zasadniczo to robi ten program. Podobnie jak poprzedni program, większość wyników bardzo wcześnie wykrywa - osiąga 200 000 000 fałszywych trafień w ciągu kilku sekund --- a następnie znacznie zwalnia.

Kompiluj z g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Aby zweryfikować wyniki, dodaj -DVERIFY(będzie to znacznie wolniejsze).

Uruchom z flspos. Dowolny argument wiersza polecenia dla trybu pełnego.

#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

Wystarczy ustawić poprzeczkę nieco wyżej ...

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

Łatwo jest zweryfikować, że dla wszystkich 67 108 864 <= x <= 94 906 265 = piętro (sqrt (2 ⁵³ )) pary (x, 0) i (x, 1) są fałszywie dodatnie.

Dlaczego to działa : 67 108 864 = 2 ²⁶ . W związku z tym wszystkie liczby x w powyższym zakresie mają postać 2 ²⁶ + x ' dla niektórych 0 <= x' <2 ²⁶ . Dla wszystkich dodatnich e , (x + e) ² = x ² + 2xe + e ² = x ² + 2 ²⁷ e + 2x'e + e ² . Jeśli chcemy mieć
(x + e) ² = x ² + 1 , potrzebujemy co najmniej 2 ²⁷ e <= 1 , to znaczy e <= ^2-27 Ponieważ jednak mantysa liczb zmiennoprzecinkowych o podwójnej precyzji ma szerokość 52-bitową, najmniejsze e takie, że x + e> x wynosi e = 2 ^26–52 = ^2–26 . Innymi słowy, najmniejszą reprezentowalną liczbą większą niż x jest x + ^2-26, podczas gdy wynik sqrt (x ² + 1) wynosi najwyżej x + ^2-27 . Ponieważ domyślny tryb zaokrąglania IEEE-754 jest zaokrąglany do najbliższego; remis -do-parzystego, zawsze będzie zaokrąglać do x, a nigdy do x + ^2-26 (gdzie remis jest tak naprawdę istotny tylko dla x = 67 108 864, Jeżeli w ogóle. Każda większa liczba zaokrągli się do x niezależnie).

C ++, 75 000 000+

Przypomnij sobie, że 3 ² + 4 ² = 5 ² . Oznacza to, że punkt (4, 3) leży na okręgu o promieniu 5 wyśrodkowanym wokół początku. W rzeczywistości dla wszystkich liczb całkowitych n , (4n, 3n) leży na takim kole o promieniu 5n . Dla wystarczająco dużego n (mianowicie takiego, że 5n> = 2 ²⁶ ), znamy już fałszywie dodatni dla wszystkich punktów na tym okręgu: (5n, 1) . Wspaniały! To kolejne 27 797 402/5 bezpłatnych par fałszywie dodatnich właśnie tutaj! Ale po co się tu zatrzymywać? (3, 4, 5) nie jest jedyną taką trojaczką.

Ten program wyszukuje wszystkie dodatnie trojaczki całkowite (a, b, c) tak, że a ² + b ² = c ² , i w ten sposób liczy fałszywie dodatnie. Szybko osiąga 70 000 000 fałszywie dodatnich wyników, ale następnie znacznie zwalnia, gdy liczba rośnie.

Kompiluj z g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Aby zweryfikować wyniki, dodaj -DVERIFY(będzie to znacznie wolniejsze).

Uruchom z flspos. Dowolny argument wiersza polecenia dla trybu pełnego.

#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

EDYCJA: Nowy program.

Stary zużywał za dużo pamięci. Ten zmniejsza zużycie pamięci o połowę, używając gigantycznej tablicy wektorowej zamiast tablicy mieszającej. Usuwa również losowy wątek 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;
}

Uruchom z -Pargumentem, aby wydrukować punkty zamiast ich liczby.

Dla mnie w trybie zliczania zajmuje to mniej niż 2 minuty i około 5 minut z drukowaniem skierowanym do pliku (~ 4 GB), więc nie do końca było ograniczone we / wy.

Mój oryginalny program był schludny, ale większość z niego zrezygnowałem, ponieważ mógł wygenerować tylko wyniki rzędu 10 ^ 5. Poszukiwano parametryzacji postaci (x ^ 2 + Axe + B, x ^ 2 + Cx + D), (x ^ 2 + ax + b, x ^ 2 + cx + d) tak, aby dla dowolnego x, (x ^ 2 + Axe + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Po znalezieniu takiego zestawu parametrów {a, b, c, d, A, B, C, D} przystąpił do sprawdzenia wszystkich wartości x poniżej maksimum. Patrząc na moje wyniki debugowania z tego programu, zauważyłem pewną parametryzację parametryzacji, która pozwoliła mi z łatwością wygenerować wiele liczb. Zdecydowałem się nie drukować numerów Ella, ponieważ miałem wiele własnych. Mam nadzieję, że teraz ktoś nie wydrukuje obu naszych zestawów liczb i nie ogłosi się zwycięzcą :)

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

Skanowanie koła Java, Bresenham-esque

Heurystycznie spodziewam się więcej kolizji, zaczynając od szerszego końca pierścienia. Spodziewałem się ulepszenia, wykonując jeden skan dla każdego kolizji, rejestrując wartości, które surplussą pomiędzy 0i r2max - r2włącznie, ale w moich testach okazało się wolniejsze niż ta wersja. Podobnie próbuje użyć pojedynczego int[]bufora, a nie wielu tworzenia tablic i list dwuelementowych. Optymalizacja wydajności jest naprawdę dziwną bestią.

Uruchom z argumentem wiersza polecenia dla danych wyjściowych par i bez prostych obliczeń.

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

Większość z nich jest taka sama, jak pokazuje Ell , a pozostałe oparte są na (j,0) (k,l). Dla każdego jwracam kilka kwadratów i sprawdzam, czy reszta daje fałszywie dodatni wynik. Zajmuje to w zasadzie cały czas z jedynie 25k (około 0,1%) zyskiem w stosunku do zwykłego (j,0) (j,1), ale zysk jest zyskiem.

To skończy się za dziesięć minut na moim komputerze, ale nie wiem, co masz. Ponieważ powody, jeśli nie skończy się, zanim skończy się czas, będzie miał znacznie gorszy wynik. W takim przypadku możesz dostosować dzielnik na linii 8, aby skończył się w czasie (to po prostu określa, jak daleko wraca dla każdego j). W przypadku niektórych różnych dzielników wyniki są następujące:

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

Aby włączyć wyjście dla każdego meczu (i, na Boga, jest to powolne, jeśli to zrobisz), po prostu odkomentuj linie 10 i 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);
    }
}

Dla porównania pierwsze 20 wyników, jakie daje (dla dzielnika = 7, wyłączając (j,0)(j,1)typy) to:

(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 fałszywych trafień

Oto bardzo naiwne wyszukiwanie siłowe, które można traktować jako implementację referencyjną.

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

Możesz wydrukować pary (i ich dokładną kwadratową wielkość), usuwając komentarz z @printflinii.

Zasadniczo rozpoczyna się wyszukiwanie x = y = 6e7pierwszej pary współrzędnych i skanowanie około 1% drogi do osi x przed zmniejszeniem x. Następnie dla każdej takiej pary współrzędnych sprawdza cały łuk tej samej wielkości (zaokrąglając w górę i w dół) pod kątem kolizji.

Kod zakłada, że jest prowadzony w systemie 64-bitowym, tak że typy domyślne całkowite i zmiennoprzecinkowe są te, 64-bitowe (jeśli nie, możesz utworzyć je int64()i float64()konstruktorów).

To daje skromne 530 wyników.