C ++, 275 000 000+
Będziemy odnosić się do par, których wielkość jest dokładnie reprezentowalna, takich jak (x, 0) , jako uczciwe pary, a do wszystkich innych par jako nieuczciwe pary wielkości m , gdzie m jest błędnie zgłoszoną wielkością pary. Pierwszy program w poprzednim poście wykorzystywał zestaw ściśle ze sobą powiązanych par uczciwych i nieuczciwych par: odpowiednio
(x, 0) i (x, 1) dla wystarczająco dużego x. Drugi program wykorzystywał ten sam zestaw nieuczciwych par, ale rozszerzył zestaw uczciwych par, szukając wszystkich uczciwych par integralnej wielkości. Program nie kończy się w ciągu dziesięciu minut, ale większość swoich wyników znajduje bardzo wcześnie, co oznacza, że większość czasu działania marnuje się. Zamiast szukać coraz rzadziej uczciwych par, ten program wykorzystuje wolny czas, aby zrobić następną logiczną rzecz: rozszerzenie zestawu nieuczciwych par.
Z poprzedniego postu wiemy, że dla wszystkich wystarczająco dużych liczb całkowitych r , sqrt (r 2 + 1) = r , gdzie sqrt jest zmiennoprzecinkową funkcją pierwiastka zmiennoprzecinkowego. Nasz plan ataku polega na znalezieniu par P = (x, y) takich, że x 2 + y 2 = r 2 + 1 dla pewnej dużej liczby całkowitej r . To dość proste, ale naiwne szukanie pojedynczych takich par jest zbyt wolne, aby było interesujące. Chcemy znaleźć te pary masowo, tak jak robiliśmy to dla uczciwych par w poprzednim programie.
Niech { v , w } będzie parą wektorów ortonormalnych. Dla wszystkich rzeczywistych skalarów r , || r v + w || 2 = r 2 + 1 . W ℝ 2 jest to bezpośredni wynik twierdzenia Pitagorasa:
Szukamy wektorów v i w takich, że istnieje liczba całkowita r, dla której x i y są również liczbami całkowitymi. Na marginesie, zauważmy, że zestaw nieuczciwych par, których użyliśmy w poprzednich dwóch programach, był po prostu szczególnym przypadkiem tego, gdzie { v , w } było standardową podstawą ℝ 2 ; tym razem chcemy znaleźć bardziej ogólne rozwiązanie. To jest, gdy pitagorejską tryplety (liczby całkowite tryplety (a, b, c) , spełniających w 2 + b 2 = C 2, które wykorzystaliśmy w poprzednim programie) powracają.
Niech (a, b, c) będzie triolą pitagorejską. Wektory v = (b / c, a / c) i w = (-a / c, b / c) (a także
w = (a / c, -b / c) ) są ortonormalne, co łatwo zweryfikować . Jak się okazało, dla każdego wyboru Pitagorasa triplet, istnieje całkowita R w taki sposób, x i y są liczbami całkowitymi. Aby to udowodnić i skutecznie znaleźć r i P , potrzebujemy małej teorii liczb / grup; Oszczędzę szczegóły. Tak czy inaczej, załóżmy, że mamy naszą całkę r , x i y . Wciąż brakuje nam kilku rzeczy: potrzebujemy rbyć wystarczająco duże i chcemy szybkiej metody, aby uzyskać z niej wiele podobnych par. Na szczęście istnieje prosty sposób na osiągnięcie tego.
Zauważ, że rzut P na v jest r v , stąd r = P · v = (x, y) · (b / c, a / c) = xb / c + ya / c , wszystko to oznacza, że xb + ya = rc . W rezultacie dla wszystkich liczb całkowitych n , (x + bn) 2 + (y + an) 2 = (x 2 + y 2 ) + 2 (xb + ya) n + (a 2 + b 2 ) n 2 = ( r 2 + 1) + 2 (rc) n + (c 2 ) n 2 = (r + cn) 2 + 1. Innymi słowy, kwadratowa wielkość par postaci
(x + bn, y + an) wynosi (r + cn) 2 + 1 , a dokładnie takich par szukamy! Dla wystarczająco dużego n są to nieuczciwe pary wielkości r + cn .
Zawsze miło jest spojrzeć na konkretny przykład. Jeśli weźmiemy triplet Pitagorasa (3, 4, 5) , to przy r = 2 mamy P = (1, 2) (możesz sprawdzić, że (1, 2) · (4/5, 3/5) = 2 i oczywiście 1 2 + 2 2 = 2 2 + 1 ). Dodanie 5 do r oraz (4, 3) do P prowadzi nas do r '= 2 + 5 = 7 i P' = (1 + 4, 2 + 3) = (5, 5) . Lo i oto 5 2 + 5 2 = 7 2 + 1. Kolejne współrzędne to r '' = 12 i P '' = (9, 8) , i znowu 9 2 + 8 2 = 12 2 + 1 itd. I tak dalej ...
Gdy r jest wystarczająco duże, zaczynamy otrzymywać nieuczciwe pary o przyrostach wielkości 5 . To mniej więcej 27 797 402/5 nieuczciwych par.
Więc teraz mamy wiele nieuczciwych par integralnej wielkości. Możemy z łatwością połączyć je z uczciwymi parami pierwszego programu, aby utworzyć fałszywie dodatnie, a przy należytej staranności możemy również użyć uczciwych par drugiego programu. Zasadniczo to robi ten program. Podobnie jak poprzedni program, większość wyników bardzo wcześnie wykrywa - osiąga 200 000 000 fałszywych trafień w ciągu kilku sekund --- a następnie znacznie zwalnia.
Kompiluj z g++ flspos.cpp -oflspos -std=c++11 -msse2 -mfpmath=sse -O3. Aby zweryfikować wyniki, dodaj -DVERIFY(będzie to znacznie wolniejsze).
Uruchom z flspos. Dowolny argument wiersza polecenia dla trybu pełnego.
#include <cstdio>
#define _USE_MATH_DEFINES
#undef __STRICT_ANSI__
#include <cmath>
#include <cfloat>
#include <vector>
#include <iterator>
#include <algorithm>
using namespace std;
/* Make sure we actually work with 64-bit precision */
#if defined(VERIFY) && FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1
# error "invalid FLT_EVAL_METHOD (did you forget `-msse2 -mfpmath=sse'?)"
#endif
template <typename T> struct widen;
template <> struct widen<int> { typedef long long type; };
template <typename T>
inline typename widen<T>::type mul(T x, T y) {
return typename widen<T>::type(x) * typename widen<T>::type(y);
}
template <typename T> inline T div_ceil(T a, T b) { return (a + b - 1) / b; }
template <typename T> inline typename widen<T>::type sq(T x) { return mul(x, x); }
template <typename T>
T gcd(T a, T b) { while (b) { T t = a; a = b; b = t % b; } return a; }
template <typename T>
inline typename widen<T>::type lcm(T a, T b) { return mul(a, b) / gcd(a, b); }
template <typename T>
T div_mod_n(T a, T b, T n) {
if (b == 0) return a == 0 ? 0 : -1;
const T n_over_b = n / b, n_mod_b = n % b;
for (T m = 0; m < n; m += n_over_b + 1) {
if (a % b == 0) return m + a / b;
a -= b - n_mod_b;
if (a < 0) a += n;
}
return -1;
}
template <typename T> struct pythagorean_triplet { T a, b, c; };
template <typename T>
struct pythagorean_triplet_generator {
typedef pythagorean_triplet<T> result_type;
private:
typedef typename widen<T>::type WT;
result_type p_triplet;
WT p_c2b2;
public:
pythagorean_triplet_generator(const result_type& triplet = {3, 4, 5}) :
p_triplet(triplet), p_c2b2(sq(triplet.c) - sq(triplet.b))
{}
const result_type& operator*() const { return p_triplet; }
const result_type* operator->() const { return &p_triplet; }
pythagorean_triplet_generator& operator++() {
do {
if (++p_triplet.b == p_triplet.c) {
++p_triplet.c;
p_triplet.b = ceil(p_triplet.c * M_SQRT1_2);
p_c2b2 = sq(p_triplet.c) - sq(p_triplet.b);
} else
p_c2b2 -= 2 * p_triplet.b - 1;
p_triplet.a = sqrt(p_c2b2);
} while (sq(p_triplet.a) != p_c2b2 || gcd(p_triplet.b, p_triplet.a) != 1);
return *this;
}
result_type operator()() { result_type t = **this; ++*this; return t; }
};
int main(int argc, const char* argv[]) {
const bool verbose = argc > 1;
const int min = 1 << 26;
const int max = sqrt(1ll << 53);
const size_t small_triplet_count = 1000;
vector<pythagorean_triplet<int>> small_triplets;
small_triplets.reserve(small_triplet_count);
generate_n(
back_inserter(small_triplets),
small_triplet_count,
pythagorean_triplet_generator<int>()
);
int found = 0;
auto add = [&] (int x1, int y1, int x2, int y2) {
#ifdef VERIFY
auto n1 = sq(x1) + sq(y1), n2 = sq(x2) + sq(y2);
if (x1 < y1 || x2 < y2 || x1 > max || x2 > max ||
n1 == n2 || sqrt(n1) != sqrt(n2)
) {
fprintf(stderr, "Wrong false-positive: (%d, %d) (%d, %d)\n",
x1, y1, x2, y2);
return;
}
#endif
if (verbose) printf("(%d, %d) (%d, %d)\n", x1, y1, x2, y2);
++found;
};
int output_counter = 0;
for (int x = min; x <= max; ++x) add(x, 0, x, 1);
for (pythagorean_triplet_generator<int> i; i->c <= max; ++i) {
const auto& t1 = *i;
for (int n = div_ceil(min, t1.c); n <= max / t1.c; ++n)
add(n * t1.b, n * t1.a, n * t1.c, 1);
auto find_false_positives = [&] (int r, int x, int y) {
{
int n = div_ceil(min - r, t1.c);
int min_r = r + n * t1.c;
int max_n = n + (max - min_r) / t1.c;
for (; n <= max_n; ++n)
add(r + n * t1.c, 0, x + n * t1.b, y + n * t1.a);
}
for (const auto t2 : small_triplets) {
int m = div_mod_n((t2.c - r % t2.c) % t2.c, t1.c % t2.c, t2.c);
if (m < 0) continue;
int sr = r + m * t1.c;
int c = lcm(t1.c, t2.c);
int min_n = div_ceil(min - sr, c);
int min_r = sr + min_n * c;
if (min_r > max) continue;
int x1 = x + m * t1.b, y1 = y + m * t1.a;
int x2 = t2.b * (sr / t2.c), y2 = t2.a * (sr / t2.c);
int a1 = t1.a * (c / t1.c), b1 = t1.b * (c / t1.c);
int a2 = t2.a * (c / t2.c), b2 = t2.b * (c / t2.c);
int max_n = min_n + (max - min_r) / c;
int max_r = sr + max_n * c;
for (int n = min_n; n <= max_n; ++n) {
add(
x2 + n * b2, y2 + n * a2,
x1 + n * b1, y1 + n * a1
);
}
}
};
{
int m = div_mod_n((t1.a - t1.c % t1.a) % t1.a, t1.b % t1.a, t1.a);
find_false_positives(
/* r = */ (mul(m, t1.c) + t1.b) / t1.a,
/* x = */ (mul(m, t1.b) + t1.c) / t1.a,
/* y = */ m
);
} {
int m = div_mod_n((t1.b - t1.c % t1.b) % t1.b, t1.a, t1.b);
find_false_positives(
/* r = */ (mul(m, t1.c) + t1.a) / t1.b,
/* x = */ m,
/* y = */ (mul(m, t1.a) + t1.c) / t1.b
);
}
if (output_counter++ % 50 == 0)
printf("%d\n", found), fflush(stdout);
}
printf("%d\n", found);
}