Тест на простоту: что можно улучшить?

[kthxbye]

Пусть

$$\operatorname{wt}(n,m)=\operatorname{wt}\left(\left\lfloor\frac{n}{m}\right\rfloor,m\right)+n\bmod m, \operatorname{wt}(0,m)=0$$

Тогда имеем последовательность

$$a(n,m)=(1+\operatorname{wt}(n,m))a\left(\left\lfloor\frac{n}{m}\right\rfloor,m\right), a(0,m)=1$$

Осуществим трансформацию

$$b(n,m)=\sum\limits_{k=0}^{n}(-1)^{\operatorname{wt}(n,m)-\operatorname{wt}(k,m)}\operatorname{if}(\binom{n}{k}\bmod m==0,0,1)a(k,m)$$

Тогда:

• если $m$ - простое, члены $b(n,m)$ всегда положительные
• если $m=2(2k+1)$ и номер $n$ первого отрицательного члена $b(n,m)$ равен $2(2k+1)$, то $2k+1$ либо простое, либо как минимум (потому что неизвестно, какие еще числа там могут вылезти) куб простого

Последнее отражено в следующей программе:

default(parisizemax,10^9)
n=7
wt(n,m)=if(n==0,0,wt(n\m,m)+n%m)
a(n,m)=if(n==0,1,(1+wt(n,m))*a(n\m,m))
b(n,m)=sum(k=0,n,(-1)^(wt(n,m)-wt(k,m))*(if(binomial(n,k)%m==0,0,1))*a(k,m))
c(n,m) = {my(k=0); while (sign(b(k,m)) != -1, k++); k; }
v=vector(2^n,i,if(c(0,2*(2*i+1))==2*(2*i+1),(2*i+1),0))
print(v)

Вот результат ее работы:

[3, 5, 7, 0, 11, 13, 0, 17, 19, 0, 23, 0, 27, 29, 31, 0, 0, 37, 0, 41, 43, 0, 47, 0, 0, 53, 0, 0, 59, 61, 0, 0, 67, 0, 71, 73, 0, 0, 79, 0, 83, 0, 0, 89, 0, 0,0, 97, 0, 101, 103, 0, 107, 109, 0, 113, 0, 0, 0, 0, 0, 125, 127, 0, 131, 0, 0,137, 139, 0, 0, 0, 0, 149, 151, 0, 0, 157, 0, 0, 163, 0, 167, 0, 0, 173, 0, 0, 179, 181, 0, 0, 0, 0, 191, 193, 0, 197, 199, 0, 0, 0, 0, 0, 211, 0, 0, 0, 0, 0, 223, 0, 227, 229, 0, 233, 0, 0, 239, 241, 0, 0, 0, 0, 251, 0, 0, 257]

Что можно улучшить в плане скорости выдачи?

[kthxbye]

За час и 42 минуты посчиталось до 1025:

[3, 5, 7, 0, 11, 13, 0, 17, 19, 0, 23, 0, 27, 29, 31, 0, 0, 37, 0, 41, 43, 0, 47, 0, 0, 53, 0, 0, 59, 61, 0, 0, 67, 0, 71, 73, 0, 0, 79, 0, 83, 0, 0, 89, 0, 0,0, 97, 0, 101, 103, 0, 107, 109, 0, 113, 0, 0, 0, 0, 0, 125, 127, 0, 131, 0, 0,137, 139, 0, 0, 0, 0, 149, 151, 0, 0, 157, 0, 0, 163, 0, 167, 0, 0, 173, 0, 0, 179, 181, 0, 0, 0, 0, 191, 193, 0, 197, 199, 0, 0, 0, 0, 0, 211, 0, 0, 0, 0, 0, 223, 0, 227, 229, 0, 233, 0, 0, 239, 241, 0, 0, 0, 0, 251, 0, 0, 257, 0, 0, 263,0, 0, 269, 271, 0, 0, 277, 0, 281, 283, 0, 0, 0, 0, 293, 0, 0, 0, 0, 0, 0, 307,0, 311, 313, 0, 317, 0, 0, 0, 0, 0, 0, 331, 0, 0, 337, 0, 0, 0, 0, 347, 349, 0,353, 0, 0, 359, 0, 0, 0, 367, 0, 0, 373, 0, 0, 379, 0, 383, 0, 0, 389, 0, 0, 0,397, 0, 401, 0, 0, 0, 409, 0, 0, 0, 0, 419, 421, 0, 0, 0, 0, 431, 433, 0, 0, 439, 0, 443, 0, 0, 449, 0, 0, 0, 457, 0, 461, 463, 0, 467, 0, 0, 0, 0, 0, 479, 0, 0, 0, 487, 0, 491, 0, 0, 0, 499, 0, 503, 0, 0, 509, 0, 0, 0, 0, 0, 521, 523, 0, 0, 0, 0, 0, 0, 0, 0, 541, 0, 0, 547, 0, 0, 0, 0, 557, 0, 0, 563, 0, 0, 569, 571,0, 0, 577, 0, 0, 0, 0, 587, 0, 0, 593, 0, 0, 599, 601, 0, 0, 607, 0, 0, 613, 0,617, 619, 0, 0, 0, 0, 0, 631, 0, 0, 0, 0, 641, 643, 0, 647, 0, 0, 653, 0, 0, 659, 661, 0, 0, 0, 0, 0, 673, 0, 677, 0, 0, 683, 0, 0, 0, 691, 0, 0, 0, 0, 701, 0,0, 0, 709, 0, 0, 0, 0, 719, 0, 0, 0, 727, 0, 0, 733, 0, 0, 739, 0, 743, 0, 0, 0, 751, 0, 0, 757, 0, 761, 0, 0, 0, 769, 0, 773, 0, 0, 0, 0, 0, 0, 787, 0, 0, 0, 0, 797, 0, 0, 0, 0, 0, 809, 811, 0, 0, 0, 0, 821, 823, 0, 827, 829, 0, 0, 0, 0, 839, 0, 0, 0, 0, 0, 0, 853, 0, 857, 859, 0, 863, 0, 0, 0, 0, 0, 0, 877, 0, 881, 883, 0, 887, 0, 0, 0, 0, 0, 0, 0, 0, 0, 907, 0, 911, 0, 0, 0, 919, 0, 0, 0, 0, 929, 0, 0, 0, 937, 0, 941, 0, 0, 947, 0, 0, 953, 0, 0, 0, 0, 0, 0, 967, 0, 971, 0, 0, 977, 0, 0, 983, 0, 0, 0, 991, 0, 0, 997, 0, 0, 0, 0, 0, 1009, 0, 1013, 0, 0, 1019, 1021, 0, 0]

Составных всего два - это 27 и 125.

[Dmitriy40]

Судя по всему простые всегда присутствуют в массиве, а составные не являющиеся степенью простого всегда отсутствуют.
Соответственно проверять только степени простых выше 2:

stop=257; v=vector(stop\2);
{forprime(p=3,stop,
   v[p\2]=p;\\Только если нужны и простые тоже
   for(k=3,oo,
      m=p^k;
      if(m>stop, break);
      if(c(0, 2*m) == 2*m, v[m\2]=m);
   );
);}
print(v);
[3, 5, 7, 0, 11, 13, 0, 17, 19, 0, 23, 0, 27, 29, 31, 0, 0, 37, 0, 41, 43, 0, 47, 0, 0, 53, 0, 0, 59, 61, 0, 0, 67, 0, 71, 73, 0, 0, 79, 0, 83, 0, 0, 89, 0, 0, 0, 97, 0, 101, 103, 0, 107, 109, 0, 113, 0, 0, 0, 0, 0, 125, 127, 0, 131, 0, 0, 137, 139, 0, 0, 0, 0, 149, 151, 0, 0, 157, 0, 0, 163, 0, 167, 0, 0, 173, 0, 0, 179, 181, 0, 0, 0, 0, 191, 193, 0, 197, 199, 0, 0, 0, 0, 0, 211, 0, 0, 0, 0, 0, 223, 0, 227, 229, 0, 233, 0, 0, 239, 241, 0, 0, 0, 0, 251, 0, 0, 257]
Time: 0.218s

\\Вывод составных чисел по 1025:
stop=1025; v=vector(stop);
{forprime(p=3,stop, for(k=3,oo,
   m=p^k;
   if(m>stop, break);
   if(c(0, 2*m) == 2*m, v[m]=m);
););}
print(select(x->(x>0),v,1));
Vecsmall([27, 125])
Time: 2.948s

Ускорение на порядки.

-- 27.11.2021, 18:07 --

Удобно выводить с показом какие да, а какие нет:

stop=1025; v=vector(stop);
{forprime(p=3,stop, for(k=3,oo,
      m=p^k; if(m>stop, break);
      v[m]=if(c(0, 2*m) == 2*m, m, -m);
););}
print(select(x->(x!=0),v));
[27, -81, 125, -243, -343, -625, -729]

[Dmitriy40]

Что ещё можно улучшить. Вот этот кусок кода (чуть оптимизировал, на скорость не повлияло):

b(n,m) = my(k,w); w=wt(n,m)%2; sum(k=0, n, if(binomial(n,k)%m==0, 0, wt(k,m)%2==w, a(k,m), -a(k,m)));
c(n,m) = my(k=0); while (b(k,m)>=0, k++); k;

Тут для последовательно увеличивающегося $n$ вызывается $b(n,m)$, которая снова начинает считать сумму с нуля, хотя достаточно сохранить сумму по $n-1$ и досчитать лишь одно значение для $n$. Этому мешает $binomial(n,k)$, но если его расписать в факториалах, то вроде бы там $m$ (которое всегда удвоенное простое в какой-то степени) может встретиться лишь один раз, при $k=0$, т.е. от биномиала похоже можно отказаться заменив на упрощённое выражение, которое уже не будет мешать не пересчитывать сумму каждый раз с нуля.

[kthxbye]

Dmitriy40, запустил вашим предыдущим кодом проверку до 4097, на что через час и 25 минут он вернул мне все значения с минусами, кроме 27 и 125.

Относительно второго вашего сообщения: $\binom{n}{k}\bmod 2=1$ например упрощается до A295989. Там сразу возвращаются индексы. Выше двойки случаи не рассматривал.

[Dmitriy40]

Если я правильно понимаю, то в $b(n,m)$ перебор и $n$ и соответственно $k$ идёт лишь до $m$ (или $m+2$), которое в свою очередь всегда равно $2p^k$ (удвоенной степени простого) (так $c(0,m)$ вызывается). Т.е. биномиал считается вовсе не для произвольных аргументов, а практически всегда строго меньше $m$ и потому ну никак не может на $m$ делиться, ведь у него в числите факториалы, куда простые входят ровно один раз (а степени простых вообще не входят), соответственно биномиал может делиться на $m$ лишь очень в исключительных случаях, которые можно попробовать выписать явно, учитывая ограничения на $n,k$. И если в них не будет зависимостей от $n-1$ и менее, либо их можно будет быстро проверить явно, это позволит не пересчитывать сумму каждый раз с нуля, а считать итерационно (прямо в цикле в $c(n,m)$). По идее это снижает сложность с квадратичной до линейной. Тут правда по моей оценке сложность примерно пятая степень (непонятно откуда), но всё равно ...

Пробовал оптимизировать вычисление $wt(n,m)$ и $a(n,m)$ кэшированием вычисленных значений в массив, они обе явно часто вызываются с повторными аргументами, заметного эффекта не получил (жалкие пара процентов). Что странно. Виноват тормозной PARI/GP. В соседней теме показывал как устранение хвостовой рекурсии вчетверо замедляет программу.

-- 27.11.2021, 22:08 --

Если посмотреть с какими параметрами вызывается $wt(n,m)$, то там полный мрак: $m=2p^k$, а вот $n$ перебирается с редкими пропусками от $0$ до $m$, причём треугольно, т.е. квадратично (перебор с 0 до увеличивающегося верхнего предела). С другой стороны, в $b(n,m)$ от функции $wt(n,m)$ нужна лишь её чётность, а так как $m$ всегда чётная, то чётность $wt(n,m)$ всегда тождественно равна просто чётности $n$, за исключением когда $n=m$, его надо обработать отдельно. Т.е. в $b(n,m)$ оба вызова $wt(n,m), wt(k,m)$ заменимы почти просто на $n\%2, k\%2$:

\\Было
b(n,m) = my(k,w); w=wt(n,m)%2; sum(k=0, n, if(binomial(n,k)%m==0, 0, wt(k,m)%2==w, a(k,m), -a(k,m)));
\\Стало
b(n,m) = my(k,w); w=if(n==m, 1, n%2); sum(k=0, n, if(binomial(n,k)%m==0, 0, w==if(k==m, 1, k%2), a(k,m), -a(k,m)));

Скорости правда почти не добавило, процента три.
Уменьшив предел в сумме можно последнее вычисление вынести за цикл суммы и отказаться от условного оператора для $k\%2$:

b(n,m) = my(k,w); w=if(n==m, 1, n%2); sum(k=0, n-1, if(binomial(n,k)%m==0, 0, k%2==w, a(k,m), -a(k,m))) + a(n,m);

Это на скорость не влияет.

Аналогичную замену можно сделать и в $a(n,m)$, только уже по модулю $m$. И нужда в $wt(n,m)$ отпадает.

-- 27.11.2021, 22:15 --

Про $binomial(n,k)\%m$ я не прав, посмотрел когда оно равно нулю, для $m=2\cdot3^3=54$ таких случаев 110 штук, причём для 18 разных значений $n$. Не представляю как это упростить чтобы считать сумму итерационно.

[Dmitriy40]

Зато вот такая замена уменьшает время счёта чисел до $2200$ с 144с до 5.3с:

b(n,m) = my(k,w,bn); w=if(n==m, 1, n%2); bn=binomial(n); sum(k=0, n-1, if(bn[k+1]%m==0, 0, w==k%2, a(k,m), -a(k,m))) + a(n,m);

До $4097$ считается за 12.9с, до $8193$ за 101с.
Предрасчёт всех возможных $a(n,m)$ и замена в коде выше обращений к функции $a(n,m)$ на выборку из массива уменьшает время с 101с до 88с.
За 808с просчиталось от $8193$ до $16385$.

PS. Положительных значений не обнаружено.

[Dmitriy40]

Теперь если внимательно посмотреть на значения $a(n,m)$, то окажется интересная вещь — все нужные значения этой функции тривиальны: $a(n<m,m)=n+1,\; a(m,m)=4$. Т.е. от неё смело можно отказаться и упростить код $b(n,m)$ до такого:

b(n,m) = my(k,w,bn); w=if(n==m, 1, n%2); bn=binomial(n)%m; sum(k=1, n, if(bn[k]==0, 0, w!=k%2, k, -k)) + if(n==m, 4, n+1);

Это уменьшает время счёта до $8193$ с 88с до 73с.
Уже ликвидировались $wt(n,m), a(n,m)$, остались $b(n,m), c(n,m)$.

[Dmitriy40]

Посмотрев на значения $binomial(n)\%m$ обнаружил что для всех $n<\lfloor m/2 \rfloor = p^s$ там не встречается ни одного нуля, т.е. в принципе можно вообще запустить в $c(n,m)$ цикл по $k$ сразу с $k=p^s$, а меньшие и не считать, они всегда проходятся.
Проверил, все числа до $8193$ вместо 75с считаются жалкие 0.42с! До $22000$ считает 2.8с. До $32769$ за 5.9с. За 15с до $65537$. Программа:

b(n,m) = my(k,w,bn); w=if(n==m, 1, n%2); bn=binomial(n)%m; sum(k=1, n, if(bn[k]==0, 0, w!=k%2, k, -k)) + if(n==m, 4, n+1);
c(n,m) = my(k=m\2); while(b(k,m)>=0, k++); k;
stop=65537; v=vector(stop);
{forprime(p=3,stop, for(k=3,oo,
   x=p^k; if(x>stop, break);
   v[x]=if(c(0, 2*x) == 2*x, x, -x);
););}
print(select(x->(x!=0),v));
[27, -81, 125, -243, -343, -625, -729, -1331, -2187, -2197, -2401, -3125, -4913, -6561, -6859, -12167, -14641, -15625, -16807, -19683, -24389, -28561, -29791, -50653, -59049]
Time: 14.695s

Но вот следующие, $41^3, 43^3$ уже снова долго считаются, больше 200с, хотя и всё равно отбрасываются, но видимо далеко не сразу. А ещё следующие снова быстро отбрасываются, за 0.1с-15с каждое.

-- 28.11.2021, 03:11 --

Честно говоря, если все предыдущие оптимизации можно было строго доказать, то вот насчёт этой сомневаюсь, чисто эмпирическое наблюдение. Но похоже ложноотрицательных срабатываний быть не может, лишь ложноположительные, но положительных так больше и не нашлось. Но если найдутся, стоит их перепроверить предыдущей версией.

Пока же вот что насчиталось (вместе с затраченным временем на каждое):

-50653. Time: 0.562s
-59049. Time: 8.237s
-68921. Time: 53.587s
-78125. Time: 1.311s
-79507. Time: 130.683s
-83521. Time: 1.466s
-103823. Time: 2.294s
-117649. Time: 14.664s
-130321. Time: 3.574s
-148877. Time: 41.981s
-161051. Time: 5.475s
-177147. Time: 152.882s
-205379. Time: 2.559s
-226981. Time: 15.819s
-279841. Time: 4.868s
-300763. Time: 5.649s

Значения от 200 тысяч и далее уже для gp64, у него скорость примерно вчетверо выше.
Дальше не хватило 5ГБ памяти для $binomial(n)$.

[Dmitriy40]

Ещё одна идея для ускорения счёта, для не слишком больших чисел: $binomial(n)$ вызывается для каждого нового $m$ снова от $m/2$ и далее, хотя от $m$ он никак не зависит и мог быть посчитан для предыдущего $m$. Т.е. можно сохранять посчитанные $binomial(n)$ в массив и если они были посчитаны ранее, то брать оттуда.
Засада в том что требуемая память быстро растёт и например для хранения всех $binomial(n=1000\ldots2000)$, нужных для расчёта любых $x=m/2<1000$, нужно уже 230МБ памяти. С другой стороны, многие и многие числа отвергаются достаточно быстро и не перебирают все $n=m/2\ldots m$, значит для последовательного увеличивающихся чисел можно хранить лишь малую часть биномиалов, только несколько после $m/2$, в надежде что следующее проверяемое число будет незначительно больше и не придётся снова вычислять эти биномиалы.
Реально это должно помогать для проверки почти всех чисел до десятков тысяч: если проверяемые числа реже, то не будет перекрытия, если числа больше, то не хватит памяти хранить несколько биномиалов.

Следующая идея: можно заранее просчитать суммы чётных и нечётных $k$ для всех $n$. Потом лишь вычитать одну из другой в зависимости от чётности $n$, предварительно вычтя из каждой суммы те $a(n,m)$, которые попадают на кратные $m$ биномиалы. Довольно часто нулей (чисел кратных $m$) в биномиале мало и это позволит не считать суммы каждый раз с нуля. Т.е. это замена суммы по всему $n$ на сумму только по нулям биномиала, а остальное посчитать заранее.

[Snegovik]

Можно использовать алгоритм - решето Эратосфена. Пример на фортране для чисел до 100. Это работает много быстрее перебора. Тест с формированием 1 млн. чисел на очень старой машине завершился за 1 секунду.

Используется синтаксис Fortran

program eratosfen
implicit none
integer :: a(2:100), j, n

forall (j=2:100) a(j) = j

        n = 2   
        do while (n*n < 100)
           forall (j = n*n : 100 : n) a(j) = 0
           n = n + 1
              do while (a(n)==0)
                 n = n + 1
              end do
        end do

  n = count (a/=0); n = n + 1
  a(2:n) = pack (a, a/=0)

print 100, a(2:n)

100 format (10i5)
end program eratosfen
 

Результат работы

    2    3    5    7   11   13   17   19   23   29
   31   37   41   43   47   53   59   61   67   71
   73   79   83   89   97

[Dmitriy40]

Snegovik
Задача не в поиске простых чисел. Их можно находить в миллиардных интервалах за секунду.