Classement

 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

Contexte

Lorsque vous travaillez sur une grille 2D de coordonnées entières, vous voulez parfois savoir si deux vecteurs (avec des composants entiers) ont la même amplitude. Bien sûr, dans la géométrie euclidienne, la grandeur d'un vecteur (x,y)est donnée par

√(x² + y²)

Une implémentation naïve pourrait donc calculer cette valeur pour les deux vecteurs et comparer les résultats. Non seulement cela entraîne un calcul inutile de la racine carrée, mais cela cause également des problèmes avec les inexactitudes en virgule flottante, qui pourraient produire des faux positifs: des vecteurs dont les amplitudes sont différentes, mais où les chiffres significatifs de la représentation en virgule flottante sont tous identiques.

Aux fins de ce défi, nous définissons un faux positif comme une paire de paires de coordonnées (a,b)et (c,d)pour lequel:

• Leur ampleur au carré est différente lorsqu'elle est représentée sous forme d'entiers non signés 64 bits.
• Leur amplitude est identique lorsqu'elle est représentée sous la forme d'un nombre à virgule flottante binaire 64 bits et calculée via une racine carrée 64 bits (selon IEEE 754 ).

Par exemple, en utilisant des représentations 16 bits (au lieu de 64), la plus petite ¹ paire de vecteurs qui donne un faux positif est

(25,20) and (32,0)

Leurs grandeurs carrées au carré sont 1025et 1024. Prendre les rendements des racines carrées

32.01562118716424 and 32.0

Mais en 16 bits, les deux flottants sont tronqués 32.0.

De même, la plus petite paire de ² donnant un faux positif pour les représentations 32 bits est

(1659,1220) and (1951,659)

_{¹ "Plus petit" mesuré par leur amplitude à virgule flottante 16 bits.
² "Plus petit" mesuré par leur amplitude à virgule flottante 32 bits.}

Enfin, voici quelques cas valides 64 bits:

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

* Le dernier cas est l'un des nombreux avec la plus petite amplitude possible pour les faux positifs 64 bits.

Le défi

En moins de 10 000 octets de code, à l'aide d'un seul thread, vous devez trouver autant de faux positifs pour les nombres à virgule flottante 64 bits (binaires) dans la plage de coordonnées 0 ≤ y ≤ x(c'est-à-dire uniquement dans le premier octant du plan euclidien) de telle sorte que dans les 10 minutes . Si deux soumissions sont à égalité pour le même nombre de paires, le bris d'égalité est le temps réel nécessaire pour trouver la dernière de ces paires.x² + y² ≤ 253

Votre programme ne doit pas utiliser plus de 4 Go de mémoire à tout moment (pour des raisons pratiques).

Il doit être possible d'exécuter votre programme en deux modes: un qui génère chaque paire comme il le trouve, et un qui ne produit que le nombre de paires trouvées à la fin. Le premier sera utilisé pour vérifier la validité de vos paires (en regardant un échantillon de sorties) et le second sera utilisé pour chronométrer réellement votre soumission. Notez que l'impression doit être la seule différence. En particulier, le programme de comptage peut ne pas coder en dur le nombre de paires qu'il a pu trouver. Il doit toujours effectuer la même boucle qui serait utilisée pour imprimer tous les nombres, et ne supprimer que l'impression elle-même!

Je vais tester toutes les soumissions sur mon ordinateur portable Windows 8, veuillez donc demander dans les commentaires si vous souhaitez utiliser un langage pas trop courant.

Notez que les paires ne doivent pas être comptées deux fois lors de la commutation des première et deuxième paires de coordonnées.

Notez également que je vais exécuter votre processus via un contrôleur Ruby, qui tuera votre processus s'il n'est pas terminé après 10 minutes. Assurez-vous de sortir le nombre de paires trouvées d'ici là. Vous pouvez soit garder une trace du temps vous-même et imprimer le résultat juste avant les 10 minutes, soit vous pouvez simplement afficher le nombre de paires trouvées sporadiquement, et je prendrai le dernier nombre comme votre score.

C ++, 275 000 000+

Nous ferons référence à des paires dont la magnitude est représentable avec précision, telles que (x, 0) , en tant que paires honnêtes et à toutes les autres paires comme des paires malhonnêtes de magnitude m , où m est la magnitude mal déclarée de la paire. Le premier programme dans le post précédent a utilisé un ensemble de couples étroitement liés de paires honnêtes et malhonnêtes:
(x, 0) et (x, 1) , respectivement, pour assez grand x. Le deuxième programme a utilisé le même ensemble de paires malhonnêtes mais a étendu l'ensemble de paires honnêtes en recherchant toutes les paires honnêtes de magnitude intégrale. Le programme ne s'arrête pas dans les dix minutes, mais il trouve la grande majorité de ses résultats très tôt, ce qui signifie que la plupart du temps d'exécution est gaspillé. Au lieu de continuer à chercher des paires honnêtes de moins en moins fréquentes, ce programme utilise le temps libre pour faire la prochaine chose logique: étendre l'ensemble des paires malhonnêtes .

De l'article précédent, nous savons que pour tous les entiers suffisamment grands r , sqrt (r ² + 1) = r , où sqrt est la fonction racine carrée à virgule flottante. Notre plan d'attaque est de trouver des paires P = (x, y) telles que x ² + y ² = r ² + 1 pour un entier r assez grand . C'est assez simple à faire, mais chercher naïvement de telles paires individuelles est trop lent pour être intéressant. Nous voulons trouver ces paires en vrac, tout comme nous l'avons fait pour les paires honnêtes dans le programme précédent.

Soit { v , w } une paire orthonormée de vecteurs. Pour tous les scalaires réels r , || r v + w || ² = r ² + 1 . Dans ℝ ² , c'est un résultat direct du théorème de Pythagore:

Nous recherchons des vecteurs v et w tels qu'il existe un entier r pour lequel x et y sont également des entiers. En remarque, notons que l'ensemble des paires malhonnêtes que nous avons utilisées dans les deux programmes précédents n'était qu'un cas particulier de ceci, où { v , w } était la base standard de ℝ ² ; cette fois, nous souhaitons trouver une solution plus générale. C'est là que les triplets de Pythagore (triplets entiers (a, b, c) satisfaisant a ² + b ² = c ², que nous utilisions dans le programme précédent) font leur retour.

Soit (a, b, c) un triplet pythagoricien. Les vecteurs v = (b / c, a / c) et w = (-a / c, b / c) (et aussi
w = (a / c, -b / c) ) sont orthonormés, comme cela est facile à vérifier . En fait, pour tout choix de triplet de Pythagore, il existe un entier r tel que x et y sont des entiers. Pour le prouver et pour trouver efficacement r et P , nous avons besoin d'une petite théorie des nombres / groupes; Je vais épargner les détails. Quoi qu'il en soit, supposons que nous ayons nos r , x et y intégraux . Nous sommes encore à court de quelques choses: nous avons besoin de rêtre assez grand et nous voulons une méthode rapide pour dériver beaucoup plus de paires similaires de celle-ci. Heureusement, il existe un moyen simple d'y parvenir.

Notez que la projection de P sur v est r v , donc r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , tout cela pour dire que xb + ya = rc . Par conséquent, pour tous les entiers n , (x + bn) ² + (y + an) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1. En d'autres termes, la magnitude au carré des paires de la forme
(x + bn, y + an) est (r + cn) ² + 1 , ce qui est exactement le type de paires que nous recherchons! Pour n assez grand , ce sont des paires de magnitude malhonnêtes r + cn .

C'est toujours agréable de regarder un exemple concret. Si nous prenons le triplet de Pythagore (3, 4, 5) , alors à r = 2 nous avons P = (1, 2) (vous pouvez vérifier que (1, 2) · (4/5, 3/5) = 2 et, clairement, 1 ² + 2 ² = 2 ² + 1. ) Ajouter 5 à r et (4, 3) à P nous amène à r '= 2 + 5 = 7 et P' = (1 + 4, 2 + 3) = (5, 5) . Et voilà, 5 ² + 5 ² = 7 ² + 1. Les coordonnées suivantes sont r '' = 12 et P '' = (9, 8) , et encore, 9 ² + 8 ² = 12 ² + 1 , et ainsi de suite, et ainsi de suite ...

Une fois que r est assez grand, nous commençons à obtenir des paires malhonnêtes avec des incréments de magnitude de 5 . Cela représente environ 27 797 402/5 paires malhonnêtes.

Alors maintenant, nous avons beaucoup de paires malhonnêtes de magnitude intégrale. Nous pouvons facilement les coupler avec les paires honnêtes du premier programme pour former des faux positifs, et avec prudence, nous pouvons également utiliser les paires honnêtes du deuxième programme. C'est essentiellement ce que fait ce programme. Comme le programme précédent, il trouve lui aussi la plupart de ses résultats très tôt --- il atteint 200 000 000 de faux positifs en quelques secondes --- puis ralentit considérablement.

Compilez avec g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Pour vérifier les résultats, ajoutez -DVERIFY(ce sera notamment plus lent.)

Courez avec flspos. Tout argument de ligne de commande pour le mode détaillé.

#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

Juste pour mettre la barre un peu plus haut ...

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

Il est facile de vérifier que pour tous les 67 108 864 <= x <= 94 906 265 = étage (sqrt (2 ⁵³ )), les paires (x, 0) et (x, 1) sont des faux positifs.

Pourquoi ça marche : 67.108.864 = 2 ²⁶ . Par conséquent, tous les nombres x dans la plage ci-dessus sont de la forme 2 ²⁶ + x ' pour quelque 0 <= x' <2 ²⁶ . Pour tout e positif , (x + e) ² = x ² + 2xe + e ² = x ² + 2 ²⁷ e + 2x'e + e ² . Si nous voulons avoir
(x + e) ² = x ² + 1 nous avons besoin d'au moins 2 ²⁷ e <= 1 , c'est-à-dire e <= 2 ^-27 Cependant, comme la mantisse des nombres à virgule flottante double précision a une largeur de 52 bits, le plus petit e tel que x + e> x est e = 2 ^{26 - 52} = 2 ^-26 . En d'autres termes, le plus petit nombre représentable supérieur à x est x + 2 ^-26 tandis que le résultat de sqrt (x ² + 1) est au plus x + 2 ^-27 . Étant donné que le mode d'arrondi IEEE-754 par défaut est arrondi au plus proche; égal à égal, il arrondira toujours à x et jamais à x + 2 ^-26 (où le bris d'égalité n'est vraiment pertinent que pour x = 67,108,864, le cas échéant. Tout nombre plus grand arrondira à x indépendamment).

C ++, 75 000 000+

Rappelons que 3 ² + 4 ² = 5 ² . Cela signifie que le point (4, 3) se trouve sur le cercle de rayon 5 centré autour de l'origine. En fait, pour tous les entiers n , (4n, 3n) se trouve sur un tel cercle de rayon 5n . Pour n assez grand (à savoir tel que 5n> = 2 ²⁶ ), nous connaissons déjà un faux positif pour tous les points de ce cercle: (5n, 1) . Génial! C'est encore 27 797 402/5 paires de faux positifs gratuites! Mais pourquoi s'arrêter ici? (3, 4, 5) n'est pas le seul de ces triplets.

Ce programme recherche tous les triplets entiers positifs (a, b, c) tels que a ² + b ² = c ² , et compte les faux positifs de cette façon. Il atteint assez rapidement 70 000 000 de faux positifs, puis ralentit considérablement à mesure que les chiffres augmentent.

Compilez avec g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Pour vérifier les résultats, ajoutez -DVERIFY(ce sera notamment plus lent.)

Courez avec flspos. Tout argument de ligne de commande pour le mode détaillé.

#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: Nouveau programme.

L'ancien utilisait trop de mémoire. Celui-ci réduit de moitié l'utilisation de la mémoire en utilisant un tableau vectoriel géant au lieu d'une table de hachage. En outre, il supprime la rupture de fil aléatoire.

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

Exécutez avec un -Pargument pour imprimer les points au lieu de leur nombre.

Pour moi, cela prend moins de 2 minutes en mode comptage et environ 5 minutes avec une impression dirigée vers un fichier (~ 4 Go), donc cela n'a pas tout à fait limité les E / S.

Mon programme d'origine était soigné, mais j'en ai abandonné la majeure partie car il ne pouvait produire que de l'ordre de 10 ^ 5 résultats. Il a cherché des paramétrisations de la forme (x ^ 2 + Ax + B, x ^ 2 + Cx + D), (x ^ 2 + ax + b, x ^ 2 + cx + d) telles que pour tout x, (x ^ 2 + Ax + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Lorsqu'il a trouvé un tel ensemble de paramètres {a, b, c, d, A, B, C, D}, il a procédé à la vérification de toutes les valeurs x sous le maximum. En regardant ma sortie de débogage de ce programme, j'ai remarqué une certaine paramétrisation de la paramétrisation de la paramétrisation qui m'a permis de produire beaucoup de nombres facilement. J'ai choisi de ne pas imprimer les numéros d'Ell car j'en avais plein. Espérons que maintenant, personne n'imprimera nos deux séries de chiffres et ne prétendra être le gagnant :)

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

Analyse du cercle Java, Bresenham-esque

Heuristiquement, je m'attends à obtenir plus de collisions en commençant à l'extrémité la plus large de l'anneau. Je m'attendais à obtenir une certaine amélioration en effectuant un balayage pour chaque collision, en enregistrant des valeurs surplusentre 0et r2max - r2inclusives, mais dans mes tests, cela s'est avéré plus lent que cette version. De même, tente d'utiliser un seul int[]tampon plutôt que de créer de nombreux tableaux et listes à deux éléments. L'optimisation des performances est une étrange bête en effet.

Exécutez avec un argument de ligne de commande pour la sortie des paires et sans pour les comptes simples.

import java.util.*;

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

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

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

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

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

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

        return rv;
    }
}

Java - 27 817 255

La plupart d'entre eux sont les mêmes que ce qu'Ell montre , et les autres sont basés sur (j,0) (k,l). Pour chacun j, je recule quelques carrés et vérifie si le reste donne un faux positif. Cela prend essentiellement tout le temps avec seulement un gain de 25k (environ 0,1%) sur juste (j,0) (j,1), mais un gain est un gain.

Cela se terminera en moins de dix minutes sur ma machine, mais je ne sais pas ce que vous avez. Pour des raisons, s'il ne se termine pas avant la fin du temps imparti, il aura un score considérablement pire. Dans ce cas, vous pouvez modifier le diviseur de la ligne 8 pour qu'il se termine dans le temps (cela détermine simplement la distance parcourue pour chacun j). Pour certains diviseurs divers, les scores sont les suivants:

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

Pour activer la sortie pour chaque match (et, mon Dieu, c'est lent si vous le faites), décommentez simplement les lignes 10 et 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);
    }
}

Pour référence, les 20 premières sorties qu'il donne (pour diviseur = 7, hors (j,0)(j,1)types) sont:

(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 faux positifs

Voici une recherche de force brute très naïve, que vous pouvez voir comme une implémentation de référence.

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

Vous pouvez imprimer les paires (et leurs amplitudes exactes au carré) en décommentant la @printfligne.

Fondamentalement, cela démarre la recherche x = y = 6e7de la première paire de coordonnées et balaye environ 1% du chemin jusqu'à l'axe des x avant de décrémenter x. Ensuite, pour chacune de ces paires de coordonnées, il vérifie l'arc entier de même amplitude (arrondi vers le haut et vers le bas) pour une collision.

Le code suppose qu'il est exécuté sur un système 64 bits, de sorte que les types entier et virgule flottante par défaut sont ceux de 64 bits (sinon, vous pouvez les créer avec int64()et les float64()constructeurs).

Cela donne un maigre 530 résultats.