Bestenliste

 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

Hintergrund

Wenn Sie an einem 2D-Gitter mit Ganzzahlkoordinaten arbeiten, möchten Sie manchmal wissen, ob zwei Vektoren (mit Ganzzahlkomponenten) die gleiche Größe haben. Natürlich ist in der euklidischen Geometrie die Größe eines Vektors (x,y)gegeben durch

√(x² + y²)

Daher könnte eine naive Implementierung diesen Wert für beide Vektoren berechnen und die Ergebnisse vergleichen. Dies verursacht nicht nur eine unnötige Quadratwurzelberechnung, sondern auch Probleme mit Gleitkommaungenauigkeiten, die zu falsch positiven Ergebnissen führen können: Vektoren, deren Größen unterschiedlich sind, bei denen jedoch alle signifikanten Ziffern in der Gleitkommadarstellung identisch sind.

Für die Zwecke dieser Herausforderung, definieren wir eine falsche positive als ein Paar von Koordinatenpaaren (a,b)und (c,d)für die gilt :

• Ihre quadratische Größe unterscheidet sich, wenn sie als 64-Bit-Ganzzahlen ohne Vorzeichen dargestellt werden.
• Ihre Größe ist identisch, wenn sie als 64-Bit-Gleitkommazahl dargestellt und über eine 64-Bit-Quadratwurzel (gemäß IEEE 754 ) berechnet wird .

Beispielsweise ist bei Verwendung von 16-Bit-Darstellungen (anstelle von 64) das kleinste ^1- Paar von Vektoren, das ein falsches Positiv ergibt, gleich

(25,20) and (32,0)

Ihre quadratischen Beträge sind 1025und 1024. Die Quadratwurzelerträge nehmen

32.01562118716424 and 32.0

Aber in 16-Bit-Floats werden beide abgeschnitten 32.0.

Ebenso ist das kleinste ^2- Paar , das für 32-Bit-Darstellungen ein falsches Positiv ergibt,

(1659,1220) and (1951,659)

_{¹ "Kleinste" gemessen an der 16-Bit-Fließkomma-Größe.
² "Kleinste" gemessen an der 32-Bit-Gleitkomma-Größe.}

Zum Schluss noch eine Handvoll gültiger 64-Bit-Fälle:

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

* Der letzte Fall ist einer von mehreren mit der kleinstmöglichen Größe für 64-Bit-Fehlalarme.

Die Herausforderung

In weniger als 10.000 Byte Code finden Sie mit einem einzigen Thread so viele falsch-positive Ergebnisse für 64-Bit-Gleitkommazahlen (binär) im Koordinatenbereich 0 ≤ y ≤ x(dh nur innerhalb des ersten Oktanten der Euklidischen Ebene). so dass innerhalb von 10 Minuten . Wenn zwei Einsendungen die gleiche Anzahl von Paaren ergeben, ist der Gleichstand die tatsächliche Zeit, die benötigt wird, um das letzte dieser Paare zu finden.x² + y² ≤ 253

Ihr Programm darf zu keinem Zeitpunkt mehr als 4 GB Speicher belegen (aus praktischen Gründen).

Es muss möglich sein, Ihr Programm in zwei Modi auszuführen: einen, der jedes gefundene Paar ausgibt, und einen, der nur die Anzahl der gefundenen Paare am Ende ausgibt. Der erste wird verwendet, um die Gültigkeit Ihrer Paare zu überprüfen (indem Sie sich einige Beispiele von Ausgaben ansehen), und der zweite wird verwendet, um Ihre Einreichung tatsächlich zeitlich zu steuern. Beachten Sie, dass der Druck der einzige Unterschied sein muss. Insbesondere kann das Zählen Programm nicht die Anzahl der Paare codiert es könnte finden. Es muss immer noch dieselbe Schleife ausgeführt werden, die zum Drucken aller Zahlen verwendet wird, und der Druck selbst muss nur weggelassen werden!

Ich werde alle Einsendungen auf meinem Windows 8-Laptop testen. Bitte fragen Sie in den Kommentaren, ob Sie eine nicht allzu verbreitete Sprache verwenden möchten.

Beachten Sie, dass Paare beim Umschalten des ersten und zweiten Koordinatenpaars nicht zweimal gezählt werden dürfen .

Beachten Sie auch, dass ich Ihren Prozess über einen Ruby-Controller ausführen werde, der Ihren Prozess beendet, wenn er nicht nach 10 Minuten beendet ist. Stellen Sie sicher, dass Sie die Anzahl der bis dahin gefundenen Paare ausgeben. Sie können entweder die Zeit selbst verfolgen und das Ergebnis kurz vor Ablauf der 10 Minuten ausdrucken, oder Sie geben nur die Anzahl der sporadisch gefundenen Paare aus, und ich nehme die letzte Zahl als Ihre Punktzahl.

C ++, 275.000.000+

Wir bezeichnen Paare, deren Größe genau darstellbar ist, wie z. B. (x, 0) , als ehrliche Paare und alle anderen Paare als unehrliche Paare der Größe m , wobei m die falsch angegebene Größe des Paares ist. Das erste Programm im vorherigen Beitrag verwendete eine Reihe eng verwandter Paare von ehrlichen und unehrlichen Paaren:
(x, 0) bzw. (x, 1) für x , das groß genug ist. Das zweite Programm verwendete dieselbe Menge unehrlicher Paare, erweiterte jedoch die Menge ehrlicher Paare, indem es nach allen ehrlichen Paaren mit ganzzahliger Größe suchte. Das Programm wird nicht innerhalb von zehn Minuten beendet, findet jedoch den größten Teil seiner Ergebnisse sehr früh. Dies bedeutet, dass ein Großteil der Laufzeit verschwendet wird. Anstatt immer seltener nach ehrlichen Paaren zu suchen, nutzt dieses Programm die freie Zeit, um die nächste logische Sache zu tun: die Menge der unehrlichen Paare zu erweitern.

Von der früheren Post wissen wir , dass für alle groß genug , um ganze Zahlen r , sqrt (r ² + 1) = r , wobei sqrt die Fließkommaquadratwurzelfunktion ist. Unser Angriffsplan besteht darin, Paare P = (x, y) zu finden, so dass x ² + y ² = r ² + 1 für eine ausreichend große ganze Zahl r gilt . Das ist einfach genug, aber naiv nach solchen Paaren zu suchen, ist zu langsam, um interessant zu sein. Wir möchten diese Paare in großen Mengen finden, so wie wir es im vorherigen Programm für ehrliche Paare getan haben.

Sei { v , w } ein orthonormales Vektorpaar. Für alle reellen Skalare gilt r , || r v + w || ² = r ² + 1 . In ℝ ² ist dies ein direktes Ergebnis des Satzes von Pythagoras:

Wir suchen nach Vektoren v und w, so dass es eine ganze Zahl r gibt, für die x und y auch ganze Zahlen sind. Als Randnotiz sei angemerkt, dass die Menge der unehrlichen Paare, die wir in den beiden vorhergehenden Programmen verwendet haben, lediglich ein Sonderfall davon war, wobei { v , w } die Standardbasis von ℝ ^{2 war} ; Dieses Mal möchten wir eine allgemeinere Lösung finden. Hier sind pythagoreische Tripletts (ganzzahlige Tripletts (a, b, c), die a ² + b ² = c ² erfüllen, die wir im vorherigen Programm verwendet haben) machen ihr Comeback.

Sei (a, b, c) ein pythagoreisches Triplett. Die Vektoren v = (b / c, a / c) und w = (-a / c, b / c) (und auch
w = (a / c, -b / c) ) sind orthonormal, wie leicht zu überprüfen ist . Wie sich herausstellt, existiert für jede Wahl des pythagoreischen Triplets eine ganze Zahl r, so dass x und y ganze Zahlen sind. Um dies zu beweisen und r und P effektiv zu finden , brauchen wir eine kleine Zahl / Gruppentheorie; Ich werde die Details schonen. Angenommen, wir haben unser Integral r , x und y . Wir haben noch ein paar Kleinigkeiten: Wir brauchen rum groß genug zu sein und wir wollen eine schnelle Methode, um noch viele ähnliche Paare von diesem abzuleiten. Glücklicherweise gibt es einen einfachen Weg, dies zu erreichen.

Man beachte , dass die Projektion von P auf v ist r v , also R = P · V = (x, y) · (b / c, a / c) = XB / c + YA / c , alles dies , dass sagen xb + ya = rc . Als Ergebnis für alle ganzen Zahlen n , (x + bn) ² + (y + AN) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1. Mit anderen Worten, die quadratische Größe von Paaren der Form
(x + bn, y + an) ist (r + cn) ² + 1 , was genau die Art von Paaren ist, nach denen wir suchen! Für groß genug n sind dies unehrliche Größenpaare r + cn .

Es ist immer schön, sich ein konkretes Beispiel anzuschauen. Wenn wir das pythagoreische Triplett (3, 4, 5) nehmen , dann haben wir bei r = 2 P = (1, 2) (Sie können überprüfen, dass (1, 2) · (4/5, 3/5) = 2 und 1 ² + 2 ² = 2 ² + 1. ) Addieren von 5 zu r und (4, 3) zu P führt uns zu r '= 2 + 5 = 7 und P' = (1 + 4, 2 + 3) = (5, 5) . Siehe da, 5 ² + 5 ² = 7 ² + 1. Die nächsten Koordinaten sind r '' = 12 und P '' = (9, 8) und wieder 9 ² + 8 ² = 12 ² + 1 und so weiter und so fort ...

Sobald r groß genug ist, erhalten wir unehrliche Paare mit Größeninkrementen von 5 . Das sind ungefähr 27.797.402 / 5 unehrliche Paare.

Jetzt haben wir also viele unehrliche Paare mit ganzzahliger Größe. Wir können sie leicht mit den ehrlichen Paaren des ersten Programms koppeln, um False-Positives zu bilden, und mit der gebotenen Sorgfalt können wir auch die ehrlichen Paare des zweiten Programms verwenden. Dies ist im Grunde, was dieses Programm tut. Wie das vorherige Programm findet auch es die meisten seiner Ergebnisse sehr früh - es erreicht innerhalb weniger Sekunden 200.000.000 falsch positive Ergebnisse - und verlangsamt sich dann beträchtlich.

Kompilieren mit g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Fügen Sie hinzu, um die Ergebnisse zu überprüfen -DVERIFY(dies ist erheblich langsamer).

Laufen Sie mit flspos. Ein beliebiges Befehlszeilenargument für den ausführlichen Modus.

#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

Nur um die Messlatte etwas höher zu legen ...

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 ist leicht zu überprüfen, ob für alle 67.108.864 <= x <= 94.906.265 = floor (sqrt (2 ⁵³ )) die Paare (x, 0) und (x, 1) falsch positiv sind.

Warum es funktioniert : 67.108.864 = 2 ²⁶ . Daher haben alle Zahlen x im obigen Bereich die Form 2 ²⁶ + x ' für einige 0 <= x' <2 ²⁶ . Für alle positiven e ist (x + e) ² = x ² + 2xe + e ² = x ² + 2 ²⁷ e + 2x'e + e ² . Wenn wir
(x + e) ² = x ² + 1 haben wollen, brauchen wir mindestens 2 ²⁷ e <= 1 , ^dh e <= 2 ^-27 Da jedoch die Mantisse von Gleitkommazahlen mit doppelter Genauigkeit 52 Bit breit ist, ist das kleinste e, so dass x + e> x ist, e = 2 ^{26 - 52} = 2 ^-26 . Mit anderen Worten, die kleinste darstellbare Zahl, die größer als x ist, ist x + 2 ^-26, während das Ergebnis von sqrt (x ² + 1) höchstens x + 2 ^{-27 ist} . Da der standardmäßige IEEE-754-Rundungsmodus auf die nächste Runde eingestellt ist; Unentschieden, es wird immer auf x und nie auf x + 2 ^-26 gerundet (wobei der Unentschieden wirklich nur für x = 67,108,864 relevant ist, wenn überhaupt. Jede größere Zahl wird unabhängig davon auf x gerundet .

C ++, 75.000.000+

Denken Sie daran, dass 3 ² + 4 ² = 5 ² . Dies bedeutet, dass der Punkt (4, 3) auf dem Kreis mit dem Radius 5 liegt, der um den Ursprung zentriert ist. Eigentlich für alle ganzen Zahlen n , (4n, 3n) liegt auf einem solchen Kreis mit dem Radius 5N . Für ausreichend großes n (nämlich so, dass 5n> = 2 ²⁶ ) kennen wir bereits ein falsch positives Ergebnis für alle Punkte auf diesem Kreis: (5n, 1) . Groß! Das sind weitere 27.797.402 / 5 freie falsch-positive Paare genau dort! Aber warum hier aufhören? (3, 4, 5) ist nicht das einzige derartige Triplett.

Dieses Programm sucht nach allen positiven ganzzahligen Triplets (a, b, c), so dass a ² + b ² = c ² , und zählt auf diese Weise falsch-positive. Es kommt ziemlich schnell zu 70.000.000 Fehlalarmen, verlangsamt sich dann aber beträchtlich, wenn die Zahlen wachsen.

Kompilieren mit g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Fügen Sie hinzu, um die Ergebnisse zu überprüfen -DVERIFY(dies ist erheblich langsamer).

Laufen Sie mit flspos. Ein beliebiges Befehlszeilenargument für den ausführlichen Modus.

#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: Neues Programm.

Der alte benutzte zu viel Speicher. Dieser halbiert die Speichernutzung, indem ein riesiges Vektorarray anstelle einer Hash-Tabelle verwendet wird. Außerdem wird die zufällige Fadenkruft entfernt.

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

Führen Sie ein -PArgument aus, um die Punkte anstelle der Anzahl auszudrucken.

Bei mir dauert es im Zählmodus weniger als 2 Minuten und beim Drucken auf eine Datei (~ 4 GB) ungefähr 5 Minuten, sodass es nicht ganz zu einer I / O-Beschränkung kam.

Mein ursprüngliches Programm war ordentlich, aber ich ließ das meiste davon fallen, da es nur in der Größenordnung von 10 ^ 5 Ergebnissen produzieren konnte. Es wurde nach Parametrisierungen der Form (x ^ 2 + Axe + B, x ^ 2 + Cx + D), (x ^ 2 + Axe + b, x ^ 2 + cx + d) gesucht, so dass für jede x, (x ^ 2 + Ax + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + Ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Als ein solcher Satz von Parametern {a, b, c, d, A, B, C, D} gefunden wurde, wurden alle x-Werte unter dem Maximum geprüft. Bei der Betrachtung meiner Debug-Ausgabe dieses Programms fiel mir eine bestimmte Parametrisierung der Parametrisierung auf, die es mir ermöglichte, auf einfache Weise viele Zahlen zu erzeugen. Ich habe mich dafür entschieden, Ell's Nummern nicht auszudrucken, da ich genug eigene hatte. Hoffentlich druckt jetzt jemand nicht beide Zahlen aus und behauptet, der Gewinner zu sein :)

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

Java, Bresenham-ähnlicher Kreisscan

Aus heuristischer Sicht erwarte ich mehr Kollisionen, wenn ich am breiteren Ende des Rings beginne. Ich erwartete eine gewisse Verbesserung, indem ich für jede Kollision einen Scan durchführte, bei dem Werte surpluszwischen 0und r2max - r2einschließlich aufgezeichnet wurden, aber in meinen Tests erwies sich dies als langsamer als diese Version. In ähnlicher Weise wird versucht, einen einzelnen int[]Puffer zu verwenden, anstatt viele Arrays und Listen mit zwei Elementen zu erstellen. Leistungsoptimierung ist in der Tat ein seltsames Biest.

Führen Sie mit einem Befehlszeilenargument für die Ausgabe der Paare und ohne für einfache Zählungen aus.

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

Die meisten davon sind die gleichen wie die, die Ell zeigt , und der Rest basiert auf (j,0) (k,l). Für jedes jgehe ich ein paar Felder zurück und überprüfe, ob der Rest falsch positiv ist. Dies nimmt im Grunde genommen die gesamte Zeit in Anspruch, mit nur 25.000 (ca. 0,1%) Gewinn (j,0) (j,1), aber ein Gewinn ist ein Gewinn.

Dies wird in weniger als zehn Minuten auf meinem Computer abgeschlossen sein, aber ich weiß nicht, was Sie haben. Wenn es nicht vor Ablauf der Zeit beendet wird, hat es eine drastisch schlechtere Punktzahl. In diesem Fall können Sie den Divisor in Zeile 8 so einstellen, dass er pünktlich endet (dies bestimmt einfach, wie weit er für jede Zeile zurückgeht j). Für einige verschiedene Teiler lauten die Punkte:

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

Um die Ausgabe für jedes Match einzuschalten (und, Gott, es ist langsam, wenn Sie das tun), entfernen Sie einfach die Kommentare in den Zeilen 10 und 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);
    }
}

Als Referenz sind die ersten 20 Ausgaben (für Divisor = 7, ausgenommen (j,0)(j,1)Typen):

(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 Fehlalarme

Hier ist eine sehr naive Brute-Force-Suche, die Sie als Referenzimplementierung ansehen können.

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

Sie können die Paare (und ihre exakten quadratischen Größen) ausdrucken, indem Sie die @printfZeile auskommentieren .

Grundsätzlich startet dies die Suche x = y = 6e7nach dem ersten Koordinatenpaar und tastet ungefähr 1% des Weges zur x-Achse ab, bevor x dekrementiert wird. Dann prüft es für jedes derartige Koordinatenpaar den gesamten Bogen der gleichen Größe (Auf- und Abrunden) auf eine Kollision.

Der Code geht davon aus, dass er auf einem 64-Bit-System ausgeführt wird, sodass die standardmäßigen Ganzzahl- und Gleitkommatypen 64-Bit-Typen sind (andernfalls können Sie sie mit int64()und float64()Konstruktoren erstellen ).

Das ergibt magere 530 Ergebnisse.