Инструкция с офсайта
"In general, if p|n, n=p*(hx+k) where h is the order of 2 mod p and k is the least solution to 2^k mod p = c. If k doesn't exist then p cannot divide n. Sometimes a k exist, but GCD(h,k) contains another prime that cannot divide n, thus p cannot divide n. By iteratively choosing p and constructing n=p*(hx+k), you can drastically reduce the number of calculations needed to find the next solution. I call this the "prime and line" method because you are combining a given prime with numbers chosen from the line hx+k. This style of iterating also has the benefit of eliminating all prime n by virtue of constructing n in factored form."
у меня перевелась в такой алгоритм для таблицы подходящих простых и их линий:
Код:
Clear[Y]
Y = Table[{1, 1}, {200000}];
p = 1; i = 1;
While[p < 1000000,
s = 0;
While[s != 1 ,
p = NextPrime[p]; k = 1; h = MultiplicativeOrder[2, p];
While[k < p,
k = k + 2;
If[PowerMod[2, k, p] == 3 ,
g = GCD[h, k];
If[g == 1,
s = 1; Break[],
If[PrimeQ[g],
kk = 1; hh = MultiplicativeOrder[2, g];
While[kk < g,
kk = kk + 2;
If[PowerMod[2, kk, g] == 3 ,
gg = GCD[hh, kk];
If[gg == 1,
s = 1; Break[]]]
]; If[gg == 1 && s == 1, Break[]]
]]]]];
Y[[i]] = {p, k + h r};
i++
]
Y = Take[Y, i - 1];
Length[Y]
Y
В миллионе натуральных получилось 23279 подходящих простых.
Похоже на правду? Или можно как-то проще или оптимальнее?
Добавлено спустя 25 минут 8 секунд:Далее пытаюсь посчитать количество шагов в конкретном диапазоне для отсчета от двух простых:
Код:
Clear[m, n, r, i, j]; A = 0; B = 0;
For[m = 2, m <= 24355,
p1 = First[Y[[m]]]; w1 = Last[Y[[m]]];
w1 = w1 /. r -> i;
t = Solve[w1 == 10^9/m, i]; t = i /. t;
jn = IntegerPart[t[[1]]];
t = Solve[w1 == 10^10/m, i]; t = i /. t;
jk = IntegerPart[t[[1]]];
A = A + jk - jn;
For[n = 1, n <= m - 1,
p2 = First[Y[[n]]]; w2 = Last[Y[[n]]];
w2 = w2 /. r -> j;
r = Reduce[p1 w1 == p2 w2 && i >= 0 && j >= 0 , Integers];
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p2 w2 /. j -> r];
t = Solve[r == 10^9, j]; t = j /. t;
jn = IntegerPart[t[[1]]];
t = Solve[r == 10^10, j]; t = j /. t;
jk = IntegerPart[t[[1]]];
B = B + jk - jn;
]; Clear[r, j]; n++
]; Print[" A= ", A, " B= ", B, " m= ", m]; m++
]
или от трех:
Код:
Clear[m, n, k, r, i, j]; A = 0; B = 0; F = 0;
For[m = 10, m <= 24355,
p1 = First[Y[[m]]]; w1 = Last[Y[[m]]];
w1 = w1 /. r -> i;
t = Solve[w1 == 10^17/m, i]; t = i /. t;
jn = IntegerPart[t[[1]]];
t = Solve[w1 == 10^18/m, i]; t = i /. t;
jk = IntegerPart[t[[1]]];
A = A + jk - jn;
For[n = 9, n <= m - 1,
p2 = First[Y[[n]]]; w2 = Last[Y[[n]]];
w2 = w2 /. r -> j;
r = Reduce[p1 w1 == p2 w2 && i >= 0 && j >= 0 , Integers];
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p2 w2 /. j -> r];
t = Solve[r == 10^17, j]; t = j /. t;
jn = IntegerPart[t[[1]]];
t = Solve[r == 10^18, j]; t = j /. t;
jk = IntegerPart[t[[1]]];
B = B + jk - jn;
p12 = 1; w12 = r /. j -> i;
For[k = 8, k <= n - 1,
p3 = First[Y[[k]]]; w3 = Last[Y[[k]]];
w3 = w3 /. r -> j;
r = Reduce[p12 w12 == p3 w3 && i >= 0 && j >= 0 , Integers];
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p3 w3 /. j -> r];
t = Solve[r == 10^17, j]; t = j /. t;
jn = IntegerPart[t[[1]]];
t = Solve[r == 10^18, j]; t = j /. t;
jk = IntegerPart[t[[1]]];
F = F + jk - jn;
]; Clear[r, j]; k++
]
]; Clear[r, j]; n++
]; Print[" A= ", A, " B= ", B, " F= ", F, " m= ", m]; m++
]
в А кол-во шагов для отсчета от одного простого, в В - от двух, в F - от трех.
для последнего варианта ответ получается такой
Код:
A= 652173913043478 B= 9290888426 F= 0 m= 10
A= 1204999465869031 B= 13546192847 F= 0 m= 11
A= 1667962428831994 B= 34047480072 F= 3077981 m= 12
A= 2502068082214756 B= 49264541255 F= 8063903 m= 13
A= 2875822235038676 B= 62377364818 F= 17144590 m= 14
A= 3507401182407097 B= 98890768490 F= 29668423 m= 15
A= 3794390978325465 B= 115283878496 F= 35718623 m= 16
A= 4046491818661599 B= 119956585136 F= 46352349 m= 17
A= 4466659885888490 B= 138200044170 F= 56875466 m= 18
A= 4828250886290257 B= 150708830690 F= 64329582 m= 19
A= 4996161334051451 B= 162146936337 F= 69468740 m= 20
A= 5142932371233447 B= 172368687059 F= 74333909 m= 21
A= 5406861990001775 B= 179029208703 F= 77750331 m= 22
A= 5530692479820157 B= 188365950215 F= 82411119 m= 23
A= 5740189686524067 B= 196267259699 F= 87337943 m= 24
A= 5835427781762162 B= 202124737919 F= 89673630 m= 25
A= 6016660161947422 B= 208664909728 F= 93485568 m= 26
A= 6102570814868384 B= 213413422489 F= 95143710 m= 27
A= 6260133840078468 B= 225607491838 F= 100765194 m= 28
A= 6981865997255693 B= 254034636441 F= 117558950 m= 29
A= 7047083388560041 B= 257003285105 F= 120391424 m= 30
A= 7168557271675171 B= 262054323452 F= 123660393 m= 31
A= 7337984982518544 B= 272881730896 F= 132081550 m= 32
A= 7446641266473633 B= 277860773362 F= 135465237 m= 33
A= 7498748723629724 B= 282097014554 F= 137946903 m= 34
A= 7544997438943188 B= 284654787661 F= 140917228 m= 35
A= 7628609479076967 B= 287954441444 F= 143044154 m= 36
A= 7742274546013062 B= 295283297370 F= 149277820 m= 37
A= 7815600275196702 B= 297970950563 F= 150969641 m= 38
A= 7850994328995664 B= 299971699543 F= 153322705 m= 39
A= 7884278352664303 B= 301989450274 F= 155703296 m= 40
A= 7915637237681725 B= 303572775631 F= 156301191 m= 41
A= 7975326851688888 B= 306090648447 F= 158072131 m= 42
A= 8031742572869857 B= 308209351835 F= 159416874 m= 43
Добавлено спустя 10 минут 14 секунд:ну и сами вычисления от двух простых:
Код:
Off[Last::"normal"]; Off[PowerMod::"ninv"];
Clear[m, n, r, i, j];
For[m = 2, m <= 23279,
p1 = First[Y[[m]]]; w1 = Last[Y[[m]]];
w1 = w1 /. r -> i;
For[n = 1, n <= m-1,
p2 = First[Y[[n]]]; w2 = Last[Y[[n]]];
w2 = w2 /. r -> j;
(*Print[" p1= ",p1," w1= ",w1," p2= ",p2," w2= ",w2];*)
r = Reduce[p1 w1 == p2 w2 && i >= 0 && j >= 0 , Integers];(*Print[
r];Print[Length[r]];*)
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p2 w2 /. j -> r];(*Print[" p1= ",p1," p2= ",p2," r= ",
r];*)
t = Solve[r == 10^17, j]; t = j /. t;
jk = IntegerPart[t[[1]]];
(*If[jk==0,Break[]];*)(*Print[jk];*)
t = Solve[r == 10^16, j]; t = j /. t;
jn = IntegerPart[t[[1]]];(*Print[jn];Print[jk-jn];Print[r];*)
For[j = jn, j <= jk,
f = PowerMod[2, r, r];
If[f == 3,
Print["f= ", f, " p1= ", p1, " p2= ", p2, " r= ",
FactorInteger[r], " j= ", j]; Break[]];
j++]
]; Clear[r, j]; n++
]; m++
]
и от трех простых:
Код:
Clear[m, n, k, r, i, j]; A = 0; B = 0; F = 0;
For[m = 3, m <= 23279,
p1 = First[Y[[m]]]; w1 = Last[Y[[m]]];
w1 = w1 /. r -> i;
For[n = 2, n <= m - 1,
p2 = First[Y[[n]]]; w2 = Last[Y[[n]]];
w2 = w2 /. r -> j;
r = Reduce[p1 w1 == p2 w2 && i >= 0 && j >= 0 , Integers];
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p2 w2 /. j -> r];
p12 = 1; w12 = r /. j -> i;
For[k = 1, k <= n - 1,
p3 = First[Y[[k]]]; w3 = Last[Y[[k]]];
w3 = w3 /. r -> j;
r = Reduce[p12 w12 == p3 w3 && i >= 0 && j >= 0 , Integers];
If[Length[r] != 0,
r = Last[Last[r]];
r = r /. C[1] -> j;
r = Expand[p3 w3 /. j -> r];
t = Solve[r == 10^17, j]; t = j /. t;
jn = IntegerPart[t[[1]]];
t = Solve[r == 10^18, j]; t = j /. t;
jk = IntegerPart[t[[1]]];
For[j = jn, j <= jk,
f = PowerMod[2, r, r];
If[f == 3,
Print["f= ", f, " p1= ", p1, " p2= ", p2, " r= ",
FactorInteger[r], " j= ", j]; Break[]];
j++]
]; Clear[r, j]; k++
]
]; Clear[r, j]; n++
];m++
]
Добавлено спустя 2 минуты 34 секунды:номера простых из известных решений, для проверок
Код:
Y[[1]]
{5, 3 + 4 r}
Y[[8]]
{97, 19 + 48 r}
Y[[2]]
{19, 13 + 18 r}
Y[[5]]
{47, 19 + 23 r}
Y[[255]]
{6793, 1027 + 1698 r}
Y[[515]]
{15287, 7219 + 7643 r}
Y[[743]]
{23333, 14367 + 23332 r}
Y[[5420]]
{205319, 57357 + 102659 r}