Leaderboard

 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

Фон

При работе с двумерной сеткой целочисленных координат иногда требуется узнать, имеют ли два вектора (с целочисленными компонентами) одинаковую величину. Конечно, в евклидовой геометрии величина вектора (x,y)определяется как

√(x² + y²)

Таким образом, наивная реализация может вычислить это значение для обоих векторов и сравнить результаты. Это не только приводит к ненужным вычислениям квадратного корня, но также вызывает проблемы с неточностями с плавающей запятой, которые могут привести к ложным срабатываниям: векторам, значения которых различны, но где значащие цифры в представлении с плавающей запятой все идентичны.

Для целей этой задачи мы определяем ложное срабатывание как пару координатных пар (a,b)и (c,d)для которых:

• Их квадратная величина отличается, когда она представлена ​​в виде 64-разрядных целых чисел без знака.
• Их величина идентична, когда она представлена ​​в виде 64-разрядного двоичного числа с плавающей запятой и вычисляется через 64-разрядный квадратный корень (согласно IEEE 754 ).

Например, используя 16-битные представления (вместо 64), наименьшая ¹ пара векторов, которая дает ложное срабатывание,

(25,20) and (32,0)

Их квадратичные величины есть 1025и 1024. Принимая квадратный корень дает

32.01562118716424 and 32.0

Но в 16-битных числах оба они усекаются 32.0.

Аналогично, наименьшая ² пара, дающая ложное срабатывание для 32-битных представлений,

(1659,1220) and (1951,659)

_{¹ «Наименьший», измеренный по 16-битной величине с плавающей запятой.
² «Наименьший», измеренный по 32-битной величине с плавающей запятой.}

Наконец, вот несколько допустимых 64-битных случаев:

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

* Последний случай - один из нескольких с наименьшей возможной величиной для 64-битных ложных срабатываний.

Соревнование

В менее чем 10 000 байтов кода, используя один поток, вы найдете столько ложных срабатываний для 64-битных (двоичных) чисел с плавающей запятой в диапазоне координат 0 ≤ y ≤ x(то есть только в пределах первого октанта евклидовой плоскости) такой, что в течение 10 минут . Если два представления связаны для одного и того же числа пар, прерыватель связи - это фактическое время, необходимое для нахождения последней из этих пар.x² + y² ≤ 253

Ваша программа не должна использовать более 4 ГБ памяти в любое время (по практическим соображениям).

Должна быть предусмотрена возможность запуска вашей программы в двух режимах: один, который выводит каждую пару в том виде, в котором он ее нашел, и один, который выводит только количество найденных пар в конце. Первый будет использоваться для проверки правильности ваших пар (путем просмотра некоторой выборки выходных данных), а второй будет использоваться для фактической синхронизации вашего представления. Обратите внимание, что печать должна быть единственным отличием. В частности, программа подсчета не может жестко число пара он может найти. Он должен по-прежнему выполнять тот же цикл, который будет использоваться для печати всех чисел, и только пропустить саму печать!

Я протестирую все материалы на своем ноутбуке с Windows 8, поэтому, пожалуйста, спросите в комментариях, хотите ли вы использовать не слишком распространенный язык.

Обратите внимание, что пары не должны учитываться дважды при переключении первой и второй пары координат.

Также обратите внимание, что я буду запускать ваш процесс через контроллер Ruby, который убьет ваш процесс, если он не завершится через 10 минут. Обязательно выведите количество найденных пар к тому времени. Вы можете сами следить за временем и распечатывать результат незадолго до истечения 10 минут, или вы просто выводите число найденных пар время от времени, и я возьму последнее такое число как ваш результат.

C ++, 275 000 000+

Мы будем ссылаться на пары, чья величина точно представима, например (x, 0) , на честные пары, а на все остальные пары - на нечестные пары величины m , где m - ошибочно сообщенная величина пары. Первая программа в предыдущем посте использовала набор тесно связанных пар честных и нечестных пар:
(x, 0) и (x, 1) соответственно для достаточно большого x, Вторая программа использовала тот же набор нечестных пар, но расширила набор честных пар, ища все честные пары интегральной величины. Программа не завершается в течение десяти минут, но обнаруживает, что подавляющее большинство ее результатов очень рано, что означает, что большая часть времени выполнения тратится впустую. Вместо того, чтобы постоянно искать честные пары, эта программа использует свободное время для выполнения следующей логической задачи: расширения набора нечестных пар.

Из предыдущего поста мы знаем, что для всех достаточно больших целых чисел r , sqrt (r ² + 1) = r , где sqrt - функция квадратного корня с плавающей точкой. Наш план атаки состоит в том, чтобы найти пары P = (x, y) , для которых x ² + y ² = r ² + 1 для некоторого достаточно большого целого числа r . Это достаточно просто сделать, но наивно искать отдельные такие пары слишком медленно, чтобы быть интересным. Мы хотим найти эти пары навалом, как мы это делали для честных пар в предыдущей программе.

Пусть { v , w } - ортонормированная пара векторов. Для всех вещественных скаляров r , || r v + w || ² = r ² + 1 . В § ² это прямой результат теоремы Пифагора:

Мы ищем векторы v и w такие, что существует целое число r, для которого x и y также являются целыми числами. В качестве примечания, заметим , что множество недобросовестных пар мы использовали в предыдущих двух программ было просто частным случаем этого, где { v , ш } был стандартный базис ℝ ² ; на этот раз мы хотим найти более общее решение. Это где пифагорейские триплеты (целые триплеты (a, b, c), удовлетворяющие a ² + b ² = c ², который мы использовали в предыдущей программе) сделать их возвращение.

Пусть (a, b, c) - пифагорейский триплет. Векторы v = (b / c, a / c) и w = (-a / c, b / c) (а также
w = (a / c, -b / c) ) ортонормированы, что легко проверить , Как выясняется, для любого выбора трифлета Пифагора существует целое число r, такое что x и y являются целыми числами. Чтобы доказать это и эффективно найти r и P , нам нужна небольшая теория чисел / групп; Я собираюсь сэкономить детали. В любом случае, предположим, что у нас есть наши интегралы r , x и y . Мы по- прежнему не хватает нескольких вещей: мы должны гбыть достаточно большим, и мы хотим, чтобы быстрый метод вывел еще много подобных пар из этой. К счастью, есть простой способ сделать это.

Обратите внимание, что проекция P на v есть r v , следовательно, r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , все это говорит о том, что xb + ya = rc . В результате для всех целых чисел n , (x + bn) ² + (y + an) ² = (x ² + y ² ) + 2 (xb + ya) n + (a ² + b ² ) n ² = ( r ² + 1) + 2 (rc) n + (c ² ) n ² = (r + cn) ² + 1, Другими словами, квадратичная величина пар вида
(x + bn, y + an) равна (r + cn) ² + 1 , и это именно тот тип пар, который мы ищем! Для достаточно больших n это нечестные пары величин r + cn .

Всегда приятно взглянуть на конкретный пример. Если мы возьмем трифлет Пифагора (3, 4, 5) , то при r = 2 имеем P = (1, 2) (вы можете проверить, что (1, 2) · (4/5, 3/5) = 2 и, очевидно, 1 ² + 2 ² = 2 ² + 1. ) Добавление 5 к r и (4, 3) к P приводит нас к r '= 2 + 5 = 7 и P' = (1 + 4, 2 + 3) = (5, 5) . И вот, 5 ² + 5 ² = 7 ² + 1, Следующие координаты r '' = 12 и P '' = (9, 8) , и снова, 9 ² + 8 ² = 12 ² + 1 , и так далее, и так далее ...

Как только r станет достаточно большим, мы начнем получать нечестные пары с приращениями величины 5 . Это примерно 27 797 402/5 нечестных пар.

Так что теперь у нас есть множество нечестных пар целой величины. Мы можем легко связать их с честными парами первой программы, чтобы сформировать ложные срабатывания, и с должной осторожностью мы также можем использовать честные пары второй программы. Это в основном то, что делает эта программа. Как и предыдущая программа, она также находит большинство своих результатов очень рано - она ​​достигает 200 000 000 ложных срабатываний в течение нескольких секунд - и затем значительно замедляется.

Компилировать с g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Чтобы проверить результаты, добавьте -DVERIFY(это будет заметно медленнее.)

Беги с flspos. Любой аргумент командной строки для подробного режима.

#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

Просто чтобы установить планку чуть выше ...

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

Легко проверить, что для всех 67 108 864 <= x <= 94 906 265 = floor (sqrt (2 ⁵³ )) пары (x, 0) и (x, 1) являются ложными срабатываниями.

Почему это работает : 67 108 864 = 2 ²⁶ . Поэтому все числа x в указанном выше диапазоне имеют вид 2 ²⁶ + x ' для некоторого 0 <= x' <2 ²⁶ . Для всех положительных е , (х + е) ² = х ² + 2xe + е ² = х ² + 2 ²⁷ е + + е 2x'e ² . Если мы хотим иметь
(x + e) ² = x ² + 1, нам нужно по крайней мере 2 ²⁷ e <= 1 , то есть e <= 2 ^-27 Однако, поскольку мантисса чисел с плавающей запятой двойной точности имеет ширину 52 бита, наименьшее e, такое что x + e> x, равно e = 2 ^{26 - 52} = 2 ^-26 . Другими словами, наименьшее представимое число, большее x, равно x + 2 ^{-26, в} то время как результат sqrt (x ² + 1) не превышает x + 2 ^-27 . Так как по умолчанию режим округления IEEE-754 является округлением до ближайшего; привязка к четному , она всегда будет округляться до x, а не до x + 2 ^-26 (где разрыв связей действительно имеет значение только для x = 67,108,864, если вообще. Любое большее число будет округлено до x независимо).

C ++, 75 000 000+

Напомним, что 3 ² + 4 ² = 5 ² . Это означает, что точка (4, 3) лежит на окружности радиуса 5 с центром вокруг начала координат. На самом деле, для всех целых п , (4n, 3n) лежит на такой окружности радиуса 5n . Для достаточно большого n (а именно, такого, что 5n> = 2 ²⁶ ), мы уже знаем ложноположительное значение для всех точек этого круга: (5n, 1) . Большой! Это еще 27,797,402 / 5 свободных ложноположительных пар прямо здесь! Но зачем останавливаться на достигнутом? (3, 4, 5) - не единственный такой триплет.

Эта программа ищет все положительные целочисленные триплеты (a, b, c), такие что a ² + b ² = c ² , и подсчитывает ложные срабатывания таким образом. Он достигает 70 000 000 ложных срабатываний довольно быстро, но затем значительно замедляется по мере роста числа.

Компилировать с g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Чтобы проверить результаты, добавьте -DVERIFY(это будет заметно медленнее.)

Беги с flspos. Любой аргумент командной строки для подробного режима.

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

С ++ 11 - 100 993 667

РЕДАКТИРОВАТЬ: Новая программа.

Старый использовал слишком много памяти. Это вдвое уменьшает использование памяти, используя гигантский векторный массив вместо хеш-таблицы. Также он удаляет случайную нить.

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

Запустите с -Pаргументом, чтобы распечатать точки вместо их числа.

Для меня это занимает менее 2 минут в режиме подсчета и около 5 минут с печатью, направленной в файл (~ 4 ГБ), поэтому он не совсем ограничен вводом / выводом.

Моя оригинальная программа была аккуратной, но я отбросил большую ее часть, поскольку она могла дать только порядка 10 ^ 5 результатов. То, что он сделал, - это поиск параметризаций вида (x ^ 2 + Ax + B, x ^ 2 + Cx + D), (x ^ 2 + ax + b, x ^ 2 + cx + d), таких, что для любого x, (x ^ 2 + Ax + B) ^ 2 + (x ^ 2 + Cx + D) ^ 2 = (x ^ 2 + ax + b) ^ 2 + (x ^ 2 + cx + d) ^ 2 + 1. Когда он нашел такой набор параметров {a, b, c, d, A, B, C, D}, он приступил к проверке всех значений x ниже максимального. Глядя на мой отладочный вывод из этой программы, я заметил определенную параметризацию параметризации параметризации, которая позволила мне легко получить много чисел. Я решил не распечатывать номера Элл, так как у меня было много моих собственных. Надеюсь, теперь кто-то не распечатает оба наших набора номеров и не объявит себя победителем :)

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

Ява, брезенхамский круг

Эвристически я ожидаю получить больше столкновений, начиная с более широкого конца кольца. Я ожидал получить некоторое улучшение, выполняя одно сканирование для каждого столкновения, записывая значения для которого surplusмежду 0и r2max - r2включительно, но в моем тестировании это оказалось медленнее, чем в этой версии. Точно так же пытается использовать один int[]буфер, а не много создания двухэлементных массивов и списков. Оптимизация производительности - действительно странный зверь.

Запускайте с аргументом командной строки для вывода пар и без простого подсчета.

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

Ява - 27 817 255

Большинство из них совпадают с тем, что показывает Элл , а остальные основаны на (j,0) (k,l). Для каждого jя возвращаюсь на несколько квадратов и проверяю, дает ли остаток ложное срабатывание. Это занимает в основном все время с только 25k (около 0,1%) выигрыша по сравнению с просто (j,0) (j,1), но выигрыш - это выигрыш.

Это закончится через десять минут на моей машине, но я не знаю, что у вас есть. Потому что причины, если он не закончится до истечения времени, будет иметь значительно худший результат. В этом случае вы можете настроить делитель в строке 8, чтобы он заканчивался вовремя (это просто определяет, как далеко он идет назад для каждого j). Для некоторых различных делителей, оценки:

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

Чтобы включить вывод для каждого совпадения (и, боже, медленно, если вы делаете), просто раскомментируйте строки 10 и 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);
    }
}

Для справки, первые 20 выходов, которые он дает (для делителя = 7, исключая (j,0)(j,1)типы):

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

Юлия, 530 ложных срабатываний

Вот очень наивный поиск грубой силы, который вы можете рассматривать как справочную реализацию.

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

Вы можете распечатать пары (и их точные квадратные величины), раскомментировав @printfстроку.

В основном это начинает поиск в x = y = 6e7первой паре координат и сканирует около 1% пути к оси x перед уменьшением x. Затем для каждой такой пары координат он проверяет всю дугу одинаковой величины (округление вверх и вниз) на предмет столкновения.

Код предполагает , что он запускается на 64-битной системе, так что целые по умолчанию и с плавающей точкой типа на это 64-битные (если нет, то вы можете создать их int64()и float64()конструкторами).

Это дает скудные 530 результатов.