Точное количество простых чисел в интервале

[Yadryara]

Вот здесь Saouter, Trudigian и Demichel уточняют границу первого Литлвудова нарушения. Вроде более понятно написано. Как понял, 525 миллионов дзета-нулей им хватило.

[vicvolf]

При условии справедливости HL1 расстояние между простыми кортежами имеет экспоненциальное распределение с параметром $C/\ln^k(N)$
https://ru.wikipedia.org/wiki/%D0%AD%D0 ... 0%B8%D0%B5
Теорема Галлахера (1976):
При условии гипотезы Харди-Литтлвуда, количества простых в непересекающихся интервалах асимптотически независимы и распределены по Пуассону. Это влечёт экспоненциальность расстояний.

[Yadryara]

vicvolf в сообщении #1689604 писал(а):

При условии справедливости HL1 расстояние между простыми кортежами имеет

Занятно. По Вашей просьбе переехали из кортежной темы сюда, чтобы вновь говорить о кортежах?

Если без кортежей никак, то может зря переезжали...

[vicvolf]

Yadryara в сообщении #1689630 писал(а):

vicvolf в сообщении #1689604 писал(а):

При условии справедливости HL1 расстояние между простыми кортежами имеет

Занятно. По Вашей просьбе переехали из кортежной темы сюда, чтобы вновь говорить о кортежах?
Если без кортежей никак, то может зря переезжали...

Хорошо, давайте здесь продолжим обсуждать вопросы по точному количеству простых чисел на интервале, а в кортежной теме относительно кортежей.

[Yadryara]

Вот добился-таки от ИИ ссылки: https://www.math.uni-bonn.de/people/jbuethe/topics/AnalyticPiX.html

Как понимаю, это как раз пишет J. Büthe.

J. Büthe (в переводе) писал(а):

Аналитическое вычисление функции подсчета простых чисел

История

Аналитическое вычисление $\pi(x)$ означает использование информации о дзета-функции Римана для определения количества простых чисел, не больших x. Теория, лежащая в основе этого, восходит к Б. Риману, который сформулировал явную формулу, включающую сумму по нетривиальным нулям дзета-функции Римана. Эта формула, после корректировки постоянного члена, была впоследствии доказана как верная Г. фон Мангольдтом.

В то время, вероятно, никому бы не пришло в голову использовать эту формулу для вычисления $\pi(x)$, поскольку вычисление нулей дзета-функции Римана вручную гораздо сложнее, чем нахождение простых чисел (известно, что Риман сам вычислил первые 3 нуля, хотя и не опубликовал их). Лишь более чем через сто лет Г. Рисель и Г. Гёль исследовали применение явной формулы Римана для вычисления функции подсчета простых чисел. Но оказалось, что сумма по нулям сходится слишком медленно, чтобы эффективно вычислить $\pi(x)$.

Позже Лагариас и Одлыжко показали, что аналитическое вычисление функции подсчета простых чисел все-таки может быть выполнено эффективно. Их подход не основывался на явной формуле, а скорее на формуле Перрона, которая выражает $\pi(x)$ как криволинейный интеграл по вертикальной линии в комплексной плоскости. Здесь аналогично возникает проблема медленной сходимости, но ее можно решить, модифицировав ядро ​​функции за счет необходимости вычисления суммы по степеням простых чисел в окрестности x. Это привело к алгоритму со сложностью O(x^(1/2+ε)) (как по памяти, так и по времени) для вычисления $\pi(x)$, который асимптотически является самым быстрым из имеющихся на сегодняшний день методов. Однако они так и не реализовали этот метод, поскольку ожидается, что подразумеваемая константа будет очень большой.

Метод Лагариаса-Одлыжко был дополнительно проанализирован Галвеем, который ввел другую функцию ядра для ускорения сходимости интеграла. Он также улучшил требования к объему памяти для вычисления простых чисел в окрестности x с помощью решета Аткина-Бернштейна. Однако он также не реализовал этот метод для вычисления $\pi(x)$ для больших значений x.

Недавно появились два подхода к аналитическому вычислению $\pi(x)$ на больших высотах. Один из них предложен Дж. Франке, Т. Кляйнюнгом, А. Йостом и автором (FKBJ), а другой — Д. Платтом. Оба метода схожи тем, что в них используется часть нулей дзета-функции Римана, а не численное вычисление кривых интегралов. Однако метод Франке основан на явной формуле Вейля-Барнера, тогда как Платт использует модификацию метода Галвея. Мы (FKBJ) первыми вычислили $\pi(10^{24})$, предполагая гипотезу Римана для наших вычислений. Этот результат позже был подтвержден Д. Платтом без каких-либо подобных предположений. Между тем, мы также вычислили значение $\pi(10^{25})$ безусловно, используя новый вариант метода Франке, разработанный автором.

Хотя методы, используемые Д. Платтом и нами, схожи, их реализации основаны на разных концепциях. В реализации Платта делается сильный акцент на строгость, используется интервальная арифметика на основе библиотек с множественной точностью, что исключает необходимость анализа ошибок. Наша реализация, напротив, ориентирована на скорость, используя 128-битный тип данных с фиксированной запятой, разработанный Дж. Франке. Именно поэтому мы приходим к разным выводам по вопросу о том, произошел ли уже переход от аналитического метода к комбинаторному.

[Yadryara]

Вот по этой таблице работали с ИИ:

Yadryara в сообщении #1692155 писал(а):

Давайте посчитаем по точной формуле Римана, сколько имеется простых чисел не превышающих 100:

$\tikz[scale=.09]{ \fill[green!30!blue!20] ( 0,50) rectangle ( 20,60); \fill[green!30!grey!40] (20,50) rectangle ( 40,60); \fill[green!90!blue!50] (80, 0) rectangle (100,10); \draw[step=20cm] (0,0) grid +(100,60); \draw (0,10) -- (100,10); \draw (0,30) -- (100,30); \draw (0,50) -- (100,50); \node at (10,55){\text{30.126}}; \node at (10,45){\text{-3.083}}; \node at (10,35){\text{-1.136}}; \node at (10,25){\text{-0.336}}; \node at (10,15){\text{0.209}}; \node at (30,55){\text{-0.900}}; \node at (30,45){\text{0.070}}; \node at (30,35){\text{0.075}}; \node at (30,25){\text{0.011}}; \node at (30,15){\text{0.055}}; \node at (50,55){\text{-0.693}}; \node at (50,45){\text{0.347}}; \node at (50,35){\text{0.231}}; \node at (50,25){\text{0.139}}; \node at (50,15){\text{-0.116}}; \node at (70,55){\text{0.000}}; \node at (70,45){\text{-0.001}}; \node at (70,35){\text{-0.004}}; \node at (70,25){\text{-0.013}}; \node at (70,15){\text{0.018}}; \node at (90,55){\textbf{28.533}}; \node at (90,45){\textbf{-2.667}}; \node at (90,35){\textbf{-0.834}}; \node at (90,25){\textbf{-0.199}}; \node at (90,15){\textbf{0.166}}; \node at (10,5){\textbf{25.780}}; \node at (30,5){\textbf{-0.689}}; \node at (50,5){\textbf{-0.092}}; \node at (70,5){\textbf{0.000}}; \node at (90,5){\textbf{24.999}}; }$

И результат поначалу вдохновил:

Простые до 100         Риман    Деконволюция без контура

Mиллион нулей:     24.999951
 
5 нулей:           25.246                         25.005

Но когда перешли к простым до тысячи, деконволюция без контурного интегрирования уже проиграла:

Простые до 1000        Риман    Деконволюция без контура

Mиллион нулей:    167.999876
 
10 тысяч нулей:   168.008                        168.023

100 нулей:        167.595                        167.442

А разобраться в этом самом контурном интегрировании пока не получается, хотя в статье Дэвида Платта написано вроде подробно.

[wrest]

Yadryara
А с какой точностью вы хотите вычислять primepi ?

[Yadryara]

$< 0.5$

Ну вот например как справляется классическая формула Римана из статьи 1859-го года:

Простые до              Риман
       10^      Mиллион нулей

         1           4.000004
         2          24.999951
         3         167.999876
         4        1229.001850
         5        9592.024036
         6       78497.993882
         7      664578.769876
         8     5761455.755962*

* — последнее значение при округлении уже не даст точного $\pi(10^8)= 5761455$.

[wrest]

Yadryara
Держите :)

(Оффтоп)

\\ ============================================================
\\  platt.gp — аналитическое pi(x) методом Платта
\\  (явная формула Римана + гауссово сглаживание + обращение Мёбиуса)
\\
\\  Нулевая часть — рядами (лемма 6.1), решётная поправка —
\\  сегментами с разложением Тейлора (§5). Глобальных переменных нет.
\\
\\  БЫСТРЫЙ СТАРТ
\\  ---------------
\\    \r platt.gp
\\    \p 50
\\    platt_pi(10^8, 60, 0.04, 500)     \\ -> 5761455.00...
\\    primepi(10^8)                      \\ сверка: 5761455
\\
\\  ИНТЕРФЕЙС  (единственная функция)
\\  ---------------------------------
\\    platt_pi(x, nzeros, lam, nseg)
\\      x      — аргумент pi(x);
\\      nzeros — сколько нулей дзета использовать;
\\      lam    — параметр сглаживания;
\\      nseg   — число сегментов Тейлора в решётной поправке.
\\
\\    Нули кешируются в zeta_zeros.gp автоматически: первый запуск
\\    их находит (медленно), повторные — читают из файла (быстро).
\\
\\  ПАРАМЕТРЫ (проверенные ориентиры)
\\  ---------------------------------
\\    x ~ 10^6 : nzeros=12, lam=0.10, nseg=300
\\    x ~ 10^8 : nzeros=60, lam=0.04, nseg=500
\\    Для больших x уменьшайте lam и увеличивайте nzeros.
\\    Правило: нулей нужно до высоты ~6/lam; lam балансирует
\\    длину интеграла главного члена (~1/lam) и ширину решета (~lam).
\\
\\  Точность: погрешность << 0.5, округление platt_pi даёт точное целое.
\\ ============================================================

\\ ================= НУЛИ ДЗЕТА-ФУНКЦИИ =================

RStheta(t) = imag(lngamma(1/4 + I*t/2)) - t/2*log(Pi);
Zhardy(t)  = real(exp(I*RStheta(t)) * zeta(1/2 + I*t));

nextzero(t0, h) = {
  my(a=t0, va=Zhardy(a), b, vb);
  while(1,
    b = a + h; vb = Zhardy(b);
    if(va*vb < 0, return(solve(u=a, b, Zhardy(u))));
    a = b; va = vb;
  );
}

savezeros(v, file) = { extern(concat("rm -f ", file)); write(file, v); }
loadzeros(file)    = read(file);

\\ Первые N нулей: из кэша (List по ссылке), иначе из файла, иначе
\\ досчитать. Возвращает вектор первых N нулей, пополняет кэш и файл.
getzeros(~cache, N, file, h) = {
  my(v, t, z);
  if(#cache >= N, return(Vec(cache)[1..N]));
  if(#cache == 0,
    iferr(v = loadzeros(file), e, v = []);
    for(i=1, #v, listput(~cache, v[i]));
  );
  if(#cache >= N, return(Vec(cache)[1..N]));
  t = if(#cache > 0, cache[#cache] + h, 3.0);
  while(#cache < N,
    z = nextzero(t, h);
    listput(~cache, z);
    t = z + h;
  );
  savezeros(Vec(cache), file);
  Vec(cache)[1..N]
}


\\ ============ СГЛАЖЕННАЯ НУЛЕВАЯ ЧАСТЬ (Phihat) ============

phihat(s, x, lam) = (x^s/s)*exp(lam^2*s^2/2);

\\ int_0^Del phihat(s0+I*h) dh рядом Тейлора (лемма 6.1), без числ. интегр.
step_incr(s0, Del, x, lam, K) = {
  my(A = log(x) + s0*lam^2, ph0 = phihat(s0,x,lam));
  my(eD = exp(I*Del*A), invIA = 1/(I*A), Ik = vector(K+1), c, s=0, n, m, k);
  Ik[1] = (eD - 1)*invIA;
  for(k=1, K, Ik[k+1] = Del^k*eD*invIA - k*invIA*Ik[k]);
  for(k=0, K,
    c = 0;
    for(m=0, k\2,
      n = k - 2*m;
      c += ((-I)^n/s0^n) * ((-lam^2/2)^m/m!);
    );
    s += c*Ik[k+1];
  );
  ph0*s
}

\\ Re Phihat для всех нулей одним проходом вдоль критической прямой
Phihat_re_all(x, lam, zeros, Del, K) = {
  my(n=#zeros, T, GatZ=vector(n), G=0, t, tnext, j);
  T = zeros[n] + 15/lam;
  t = zeros[1];
  for(j=1, n,
    while(t < zeros[j] - 1e-30,
      tnext = min(t + Del, zeros[j]);
      G += step_incr(1/2 + I*t, tnext - t, x, lam, K);
      t = tnext;
    );
    GatZ[j] = G;
  );
  while(t < T,
    tnext = min(t + Del, T);
    G += step_incr(1/2 + I*t, tnext - t, x, lam, K);
    t = tnext;
  );
  vector(n, j, imag(G - GatZ[j]))
}

\\ Re Phihat(1) — главный член, проход вдоль Re(s)=1 от t=0
Phihat_re_at_1(x, lam, Del, K) = {
  my(T = 15/lam, t=0, tnext, G=0);
  while(t < T,
    tnext = min(t + Del, T);
    G += step_incr(1 + I*t, tnext - t, x, lam, K);
    t = tnext;
  );
  imag(G)
}


\\ ============ РЕШЁТНАЯ ПОПРАВКА (§5 Платта) ============

\\ sum_{p^m} (chi-phi)/m; сегменты + Тейлор для phi (erfc раз в сегмент)
sieve_corr_taylor(x, lam, nseg) = {
  my(c=5, lo=x*exp(-c*lam), hi=x*exp(c*lam));
  my(sq2lam=sqrt(2)*lam, sqrtpi=sqrt(Pi), invx=1/x);
  my(B=vector(nseg+1, j, round(lo + (j-1)*(hi-lo)/nseg)));
  my(s=0, j, a, b, x0, N, S1, S2, C, d, g, eg, ph0, ph1, ph2, p);
  for(j=1, nseg,
    a=B[j]; b=B[j+1]-1;
    if(b >= a,
      x0=(a+b)/2; N=0; S1=0; S2=0; C=0;
      forprime(p=a, b,
        if(p<x, C++);
        d=x0-p; S1+=d; S2+=d*d; N++;
      );
      g=log(x0*invx)/sq2lam; eg=exp(-g*g);
      ph0=0.5*erfc(g);
      ph1=-eg/(sqrtpi*sq2lam*x0);
      ph2=-eg/(sqrtpi*sq2lam*x0*x0)*(2*g/sq2lam+1);
      s += C - (ph0*N - ph1*S1 + ph2*S2/2);
    );
  );
  my(m=2, q);
  while(2^m<=hi,
    forprime(p=ceil(lo^(1/m)), floor(hi^(1/m)),
      q=p^m;
      if(q>=lo&&q<=hi, s+=((q<x)-0.5*erfc(log(q*invx)/sq2lam))/m);
    );
    m++;
  );
  s
}


\\ ============ ОБРАЩЕНИЕ МЁБИУСА ============

\\ точное f(y)=sum (1/n) pi(y^{1/n}) встроенным primepi
\\ (нужно только для поправок при малых аргументах y=x^{1/n}, n>=2)
ftarget(y) = sum(n=1, floor(log(y)/log(2)), primepi(sqrtn(y,n))/n);


\\ ============ СБОРКА ============

\\ сглаженное f(x); нули добывает сам через getzeros (кэш в файле)
platt_f(x, nzeros, lam, nseg) = {
  my(cache = List());
  my(zeros = getzeros(~cache, nzeros, "zeta_zeros.gp", 0.25));
  my(re1 = Phihat_re_at_1(x, lam, 0.5, 40));
  my(rez = Phihat_re_all(x, lam, zeros, 2, 24));
  my(zsum = sum(j=1, #rez, rez[j]));
  re1 - 2*zsum - log(2) + sieve_corr_taylor(x, lam, nseg)
}

\\ pi(x) = f(x) + поправки Мёбиуса при n>=2
platt_pi(x, nzeros, lam, nseg) = {
  my(s = platt_f(x, nzeros, lam, nseg), n);
  for(n=2, floor(log(x)/log(2)),
    if(moebius(n)!=0, s += (moebius(n)/n)*ftarget(sqrtn(x,n)));
  );
  s
}

Запуск

? \p 50
realprecision = 57 significant digits (50 digits displayed)
? platt_pi(10^8, 60, 0.04, 500)
time = 4,905 ms.
5761454.9920227407395004742199404525311697855001998
? primepi(10^8)
time = 31 ms.
5761455
?

Запускать в линуксе, есть побочный эффект: в домашней папке (откуда запущен gp) будет сохраняться файл zeta_zeros.gp с нулями дзета функции.

Скрипт выдавлен из Qwen-а

[Yadryara]

Молодцы

Разбираюсь. Для маленьких чисел пока работает очень здорово.

[wrest]

Yadryara в сообщении #1729979 писал(а):

Для маленьких чисел пока работает очень здорово.

Баланс по времени счёта, как говорит нам Платт, регулируется параметром $\lambda$ и этот параметр балансирует подсчет по нулям и подсчет поправки (см. первые два абзаца 7. The Computation).
С увеличением интервала надо уменьшать эту лямбду. Но конкретный её выбор зависит от среды где это вычисляется. Для интервала 10^9 у меня баланс достигается при
platt_pi(10^9, 300, 0.006, 400)
То есть $\lambda=6\cdot 10^{-3}$ и 300 нулей.
Я добавил в главную функцию печать таймингов, выглядит так:

? platt_pi(10^9, 300, 0.006, 400)
Получаем 300 нулей ... готово за 17 мс
Вычисляем главный член ... готово за 4309 мс
Вычисляем ряд phihat ... готово за 693 мс
Суммируем ряд ... готово за 0 мс
Вычисляем поправку решетом ... готово за 4407 мс
Обращение Мёбиуса ... готово за 0 мс
time = 9,376 ms.
50847534.007594498779931789572117124789251893668638
?

Функция с печатью:

(Оффтоп)

\\ сглаженное f(x); нули добывает сам через getzeros (кэш в файле)
platt_f(x, nzeros, lam, nseg) = {
  my(cache = List());
  my(res);
  my(t0=getwalltime());
  print1("Получаем ", nzeros, " нулей ...");
  my(zeros = getzeros(~cache, nzeros, "zeta_zeros.gp", 0.25));
  print(" готово за ", getwalltime()-t0, " мс");
  t0=getwalltime();
  print1("Вычисляем главный член ...");
  my(re1 = Phihat_re_at_1(x, lam, 0.5, 40));
  print(" готово за ", getwalltime()-t0, " мс");
  t0=getwalltime();
  print1("Вычисляем ряд phihat ...");
  my(rez = Phihat_re_all(x, lam, zeros, 2, 24));
  print(" готово за ", getwalltime()-t0, " мс");
  t0=getwalltime();
  print1("Суммируем ряд ...");
  my(zsum = sum(j=1, #rez, rez[j]));
  print(" готово за ", getwalltime()-t0, " мс");
  t0=getwalltime();
  print1("Вычисляем поправку решетом ...");
  res = re1 - 2*zsum - log(2) + sieve_corr_taylor(x, lam, nseg);
  print(" готово за ", getwalltime()-t0, " мс");
  res
}

\\ pi(x) = f(x) + поправки Мёбиуса при n>=2
platt_pi(x, nzeros, lam, nseg) = {
  my(s = platt_f(x, nzeros, lam, nseg), n);
  t0=getwalltime();
  print1("Обращение Мёбиуса ...");
  for(n=2, floor(log(x)/log(2)),
    if(moebius(n)!=0, s += (moebius(n)/n)*ftarget(sqrtn(x,n)));
  );
  print(" готово за ", getwalltime()-t0, " мс");
  s
}

По подсчёту интегралов тоже можно как-то ускорить наверное, но вы же не гонитесь за скоростью тут (комбинаторные методы считают $\pi(x)$ намного быстрее аналитических).

[Yadryara]

wrest в сообщении #1730000 писал(а):

Я добавил в главную функцию печать таймингов

Здорово.

Ну я прощупывал подходы, как подняться выше 10^12. Всё-таки это довольно долго. Например:

? platt_pi(99999998881.1, 1000, 0.013, 10000)
time = 9min, 52,814 ms.
%163 = 4118054760.35315458615681893438695651

99999998881 — это второе число кортежа [0, 60].

Я правильно понимаю, что 12-сегментный контурный интеграл Платта в этой программе не считается?

[wrest]

Yadryara в сообщении #1730016 писал(а):

Я правильно понимаю, что 12-сегментный контурный интеграл Платта в этой программе не считается?

Всё должно соответствать Платту.
Кроме контроля ошибок интервальной арифметикой - это не делалось, но 50 значащих цифр должно хватать.

-- добавлено через 29 минут --

Yadryara в сообщении #1730016 писал(а):

Ну я прощупывал подходы, как подняться выше 10^12.

Ну, тут вас ждёт выщипывание блох от погрешностей, если захотите ускорять.
Программа посчитала вам точно, но долго :)
Вы же хотели чего-то другого, а не вычисления $\pi(x)$, верно?
Методом Платта вы не посчитаете быстрее чем комбинаторным.
У меня на планшете $\pi(10^{18})$ считается за 26c, а $\pi(99999998881)$ считается за 50мс
Это отдельная программа (не pari), вы тоже можете её установить.
Но аналитический метод вам был нужен для чего-то другого.
Ещё раз подчеркну, кстати, что у Платта предполагается, что аналитическая часть и комбинаторная (поправки наиконце интервала) считаются примерно равное время. Но если сократить комбинаторную, то аналитическая потребует больше (больше нулей надо) и в сумме время увеличится.

~/primecount/build $ time ./primecount -g 10^18
24739954287740860

real    0m26.360s
user    2m49.736s
sys     0m1.144s
~/primecount/build $ time ./primecount -g 99999998881
4118054760

real    0m0.049s
user    0m0.078s
sys     0m0.029s
~/primecount/build $

[Dmitriy40]

wrest в сообщении #1730025 писал(а):

Но аналитический метод вам был нужен для чего-то другого.

Он надеется, разобравшись с $\pi(x)$, расширить формулы для любых кортежей и получить точную замену очень примерной HL1. Т.е. заменить перебор кандидатов в кортежи на вычисление (хотя бы очень примерно) где функция их количества делает скачок - там и будет искомый кортеж, гарантированно. HL1 такой гарантии к сожалению не даёт. Но по моему эта задача по сложности даже труднее исходной формулы Римана и тянет чуть ли не на нобелевку ...

[Yadryara]

wrest в сообщении #1730025 писал(а):

Вы же хотели чего-то другого, а не вычисления $\pi(x)$, верно?

В смысле? Имеете в виду хотел ли я чего-то кроме вычисления $\pi(x)$? Да, конечно. Например, научиться точнее считать некоторые кортежи.

wrest в сообщении #1730025 писал(а):

У меня на планшете $\pi(10^{18})$ считается за 26c, а $\pi(99999998881)$ считается за 50мс
Это отдельная программа (не pari), вы тоже можете её установить.

primesieve так быстро не работает, вроде. Другая программа?

Dmitriy40 в сообщении #1730029 писал(а):

Он надеется, разобравшись с $\pi(x)$, расширить формулы для любых кортежей и получить точную замену очень примерной HL1.

Кто он?