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

[wrest]

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

Судя по размеру текстовых файлов (кроме одного), у Одлыжко точность гораздо хуже, чем у Платта.

Да, но какая вам нужна точность нулей и почему?

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

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

Уже единолично вашей? Вроде это была ваша с Квеном программа.

Вы ещё спросите не взять Платта в соавторы, вроде программа по его статье написана

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

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

Соответственно я и не понял чем нули Одлыжко лучше или проще в использовании,

Там у меня функция, которая получает, нули - отдельно. Можно заменить на свою.
Я не думаю что в диапазонах до 10^12 будет какая-то разница.
Платт же занимался математически строгим вычислением, контролируя погрешность при помощи интервальной арифметики, а у нас тут некая proof of concept.
Впрочем, если где-тотначнёт расходиться то можно будет заменить нули на более точные и посмотреть из-за них ли.

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

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

или как Вам выше ИИ посоветовал, или одной командой z=lfunzeros(1,9878) (только надо подобрать верхний предел для нужного количества нулей, указал для 10 тысяч), причём с любой желаемой точностью (через команду \p), правда очень долго (часы). Но можно один раз посчитать и сохранить в файл, а потом уже пользоваться.

А вот как раз про эту функцию написано что

Цитата:

Use a naive algorithm which may miss some zeros.

Вот и у меня была наивная реализация, которая нули пропускала.

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

Версия 2

(Оффтоп)

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

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

load_zeros(file, n) = {
  my(v = readvec(file));
  if(#v < n, error("Не хватает нулей в файле ", file," Надо: ",n," Есть:",#v));
  Vec(v)[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 для всех нулей; адаптивный шаг Del = ratio*|1/2+I*t|

Phihat_re_all(x, lam, zeros, K) = {
  my(n=#zeros, GatZ=vector(n), G=0, t=zeros[1], j);
  GatZ[1] = 0;
  for(j=2, n,
    G += step_incr(1/2 + I*t, zeros[j] - t, x, lam, K);
    t = zeros[j];
    GatZ[j] = G;
  );
  vector(n, j, imag(G - GatZ[j]))
}

\\ Re Phihat(1): главный член, адаптивный шаг Del = ratio*|1+I*t|
Phihat_re_at_1(x, lam, ratio, K) = {
  my(T = 8/lam, t=0, tnext, G=0, Del);
  while(t < T,
    Del = ratio*sqrt(1+t^2);
    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 раз в сегмент)
\\ Ширина окна = 2*c*lam*x, c=5 фиксировано (проверено: c=5,8,10 дают одно и то же)
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
ftarget(y) = sum(n=1, floor(log(y)/log(2)), primepi(sqrtn(y,n))/n);

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

\\ сглаженное f(x)
platt_f(x, nzeros, lam, nseg) = {
  my(res);
  my(t0=getwalltime());
  print1("Получаем ", nzeros, " нулей ...");
  my(zeros = load_zeros("zeros1", nzeros));
  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, 12));
  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, " мс");
  return(res);
}

\\ pi(x) = f(x) + поправки Мёбиуса при n>=2
platt_pi(x, nzeros, lam, nseg) = {
    my(s = platt_f(x, nzeros, lam, nseg), n);
   my(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, " мс");
  return(s);
}

Запуск

~/gp-scripts $ gp -q platt_v2.gp
? platt_pi(10^11,56000,0.0001,100)
Получаем 56000 нулей ... готово за 119 мс
Вычисляем главный член ... готово за 45 мс
Вычисляем оставшийся ряд phihat ... готово за 11172 мс
Суммируем ряд ... готово за 9 мс
Вычисляем поправку решетом ... готово за 6081 мс
Обращение Мёбиуса ... готово за 1 мс
time = 17,327 ms.
4118054813.0290594618784539930078960821
?

Перед первым запуском, скачать нули в файл zeros1
Команда для линукс (запускать в той папке которую gp будет считать домашней, в случае выше это ~/gp-scripts ).

wget -O zeros1 https://www-users.cse.umn.edu/~odlyzko/zeta_tables/zeros1

Или скачать в файл с таким именем более точные нули, если надо. Файл zeros1 - текстовый, по одному числу на строку, чтобы в pari/gp работала команда v=readvec("zeros1")

[Yadryara]

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

Да, но какая вам нужна точность нулей и почему?

Понятия не имею, это опытным путём будет установлено, когда до этого по-настоящему дойдут руки. Так же как и значения других параметров.

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

Вы ещё спросите не взять Платта в соавторы, вроде программа по его статье написана

Не понял юмора. Я один-то злобный смайлик не люблю, а тут их целых три.

Что меня интересовало о том и спросил. На всякий случай: у меня по-прежнему нет привычки задавать риторические вопросы.

[wrest]

(Оффтоп)

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

Что меня интересовало о том и спросил.

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

Функция оценки погрешности

(Оффтоп)

platt_error_estimates(x, nzeros, zero_acc, lam, nseg) = {
\\ Оценка составляющих погрешности и суммарной
  my(v = readvec("zeros1"));
  if(#v < nzeros, error("Не хватает нулей в файле ", file," Надо: ",n," Есть:",#v));
  my(zeros=Vec(v)[1..nzeros]);
  my(c=5, T1=zeros[nzeros], T=8/lam);
  my(E1, E2, E3, E4, E5, E6, Etotal);

  \\ 1. Хвост нулей (Лемма A.4, консервативная оценка)
  E1 = 2*exp(lam^2*(1-T1^2)/2) * (sqrt(x)/(T1*log(x)) + 1/(lam^2*T1^2*x)) * (lam^2*T1^2+1);

  \\ 2. Хвост главного члена (Лемма A.2)
  E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));

  \\ 3. I_{-1} (Лемма 4.5)
  E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);

  \\ 4. Хвост решета (вне окна [x*e^{-c*lam}, x*e^{c*lam}])
  E4 = x*lam/log(x) * exp(-c^2/2)/(sqrt(Pi)*c);

  \\ 5. Проверка нулей по формуле N(T)
  my(N_exp = T1/(2*Pi)*log(T1/(2*Pi)) - T1/(2*Pi) + 7/8);
  my(bound = 0.137*log(T1) + 0.443*log(log(T1)) + 1.588);
  my(diff = abs(N_exp - nzeros));
  E5 = if(diff > bound, diff, 0);

  \\ 6. Ошибка от неточности нулей
  E6 = zero_acc * sum(j=1, nzeros,
    my(g=zeros[j]);
    sqrt(x)/sqrt(1/4+g^2)*exp(lam^2*(1/4-g^2)/2)
  );

  Etotal = E1 + E2 + E3 + E4 + E5 + E6;

  print("Составляющие погрешности:");
  print("  1. Хвост нулей:        ", E1);
  print("  2. Хвост главн. члена: ", E2);
  print("  3. I_{-1}:             ", E3);
  print("  4. Хвост решета:       ", E4);
  print("  5. Пропущенные нули:   ", E5);
  print("  6. Неточность нулей:   ", E6);
  print("  Суммарная:             ", Etotal);
  if(Etotal < 0.5,
    print("  OK: погрешность < 0.5, можно округлять"),
    print("  ВНИМАНИЕ: погрешность >= 0.5")
  );
}

Для работы нужен файл zeros1 с нулями , заодно оцнивается есть ли в нём пропуски нулей (не во всём файле, а только до запрашиваемого количества, 56000 в случае запуска выше).

Запуск:

? platt_error_estimates(10^11, 56000, 4e-9, 1e-4, 1e2);
Составляющие погрешности:
  1. Хвост нулей:        0.00053543048285827827003239889931909372067
  2. Хвост главн. члена: 6.2499741402053998968848197263833762216 E-10
  3. I_{-1}:             0.00031834163459346596101071978766488784684
  4. Хвост решета:       0.16602200777641179714083432175110640678
  5. Пропущенные нули:   0
  6. Неточность нулей:   0.0055770500119285549944178360390794870799
  Суммарная:             0.17245283053078951038683526616565184806
  OK: погрешность < 0.5, можно округлять
cpu time = 283 ms, real time = 288 ms.
?

То есть при точности нулей 4e-9 и 5.6e4 нулях оценка погрешности от них получается 6e-3 (если все нули с ошибкой 4e-9 в одну сторону).
А от решета на конце интервала -- погрешность 0,17 то есть в 30 раз больше.
Параметры (числа):
1. Интервал
2. Кол-во нулей
3. Точность нулей
4. Параметр сглаживания ($\lambda$ у Платта)
5. Количество сегментов ряда Тейлора для интеграла по нулям

[Dmitriy40]

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

чтобы в pari/gp работала команда v=readvec("zeros1")

Я конечно понимаю что одна команда чтения это удобно, но если надо всего 2000 нулей, то зачем читать из файла миллион? Они ещё и в память могут не влезть ... Лучше бы функцию чтения сделать циклом строго до требуемого количества.

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

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

my(bound = 0.137*log(T1) + 0.443*log(log(T1)) + 1.588);

Очевидно точность этих констант (если они не равны в точности вот именно этому) недостаточна при количестве нулей более тысячи.

[wrest]

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

Очевидно точность этих констант (если они не равны в точности вот именно этому) недостаточна при количестве нулей более тысячи).

Ну, Платту хватило :) См. стр 13:

Цитата:

Lemma A.1. Let $t\ge 2$. Then
$ \left| N(t) - \dfrac{t}{2\pi} \log\left(\dfrac{t}{2\pi e}\right) - \dfrac{7}{8} \right| < 0.137 \log(t) + 0.443 \log\log(t) + 1.588$
Proof. See [17] Theorem 19.

[Dmitriy40]

Да, был не прав, не количество нулей больше тысячи, а $\log t$ больше тысячи, а $t$ похоже наибольший ноль, так что его логарифм намного меньше тысячи.

[wrest]

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

Я конечно понимаю что одна команда чтения это удобно, но если надо всего 2000 нулей, то зачем читать из файла миллион? Они ещё и в память могут не влезть ... Лучше бы функцию чтения сделать циклом строго до требуемого количества.

Может и так :) Два миллиона нулей Одлыжко вроде влазит, и хорошо.
Можно поправить, потом, если надо будет читать из большого файла, или например будет надо немного нулей, но с большими номерами.
Сейчас вся интрига - как же конвертировать это всё для кортежей...
Последние тайминги - 20сек для диапазона 1e11 - я считаю уже норм. Дальше можно начать ускорение внимательным разбором типов, повторных вычислений, параллелизацией и так далее, и получить ещё некоторое ускорение. Встроенная в pari /gp primepi() считает это за 4сек, то есть отставание Платта уже всего на порядок. Но конечно, с учётом того что нули предвычислены, иначе это были бы десятки минут, а не доли минуты.

Тут ещё раз хотелось бы обратить внимание на то, что Платт, для вычисления $\pi(10^{24})$ решетил диапазон в примерно $10^{15}$ простых на конце интервала, и это решето отняло половину времени вычисления. В кортежных темах исследования идут на высоте $10^{40}$ и выше, а это означает невозможность точного вычисления там, а приближённые формулы уже дадены в HL1.

Хотя, если Yadryara как-то применит всё это к CPAP-7, минимальный из которых ожидается где-то ниже $10^{23}$ -- ну будет наверное круто

[Yadryara]

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

В кортежных темах исследования идут на высоте $10^{40}$ и выше, а это означает невозможность точного вычисления там, а приближённые формулы уже дадены в HL1.

Не согласен с таким обобщением. Вот например для [0, 60] интересна зона вблизи 106-го триллиона, то бишь порядок $10^{14}$. Я, кстати, досчитываю реальное количество кортежей в 0 — 130 триллионов.

[wrest]

Улучшенная оценка погрешности. По хвостам нулей считается реальный интеграл, а не формула из Платта.
Затем хвост делится на 10, т.к. нули осциллируют. Но надёжная ("с гарантией") оценка в 10 раз хуже, то есть скорее всего если оценка проходит, то и подсчёт будет точный (при округлении до целого), но гарантий нет.
Деление на 10 -- чисто волюнтаристское решение.

(Оффтоп)

platt_error_estimates(x, nzeros, zero_acc, lam, nseg) = {
\\ Оценка составляющих погрешности и суммарной
  my(v = readvec("zeros1"));
  if(#v < nzeros, error("Не хватает нулей в файле ", file," Надо: ",n," Есть:",#v));
  my(zeros=Vec(v)[1..nzeros]);
  my(c=5, T1=zeros[nzeros], T=8/lam);
  my(E1, E2, E3, E4, E5, E6, Etotal);

  \\ 1. Хвост нулей (Лемма A.4, консервативная оценка)
  \\ E1 = 2*exp(lam^2*(1-T1^2)/2) * (sqrt(x)/(T1*log(x)) + 1/(lam^2*T1^2*x)) * (lam^2*T1^2+1);
  \\ 1. Хвост нулей (численное интегрирование, всегда корректно)
  E1 = 2*intnum(t=T1, 8/lam, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
  \\ В print_error_estimates, после вычисления E1:
  E1 = E1 / 10;  \\ эмпирический коэффициент сокращения осцилляций

  \\ 2. Хвост главного члена (Лемма A.2)
  E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));

  \\ 3. I_{-1} (Лемма 4.5)
  E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);

  \\ 4. Хвост решета (вне окна [x*e^{-c*lam}, x*e^{c*lam}])
  E4 = x*lam/log(x) * exp(-c^2/2)/(sqrt(Pi)*c);

  \\ 5. Проверка нулей по формуле N(T)
  my(N_exp = T1/(2*Pi)*log(T1/(2*Pi)) - T1/(2*Pi) + 7/8);
  my(bound = 0.137*log(T1) + 0.443*log(log(T1)) + 1.588);
  my(diff = abs(N_exp - nzeros));
  E5 = if(diff > bound, diff, 0);

  \\ 6. Ошибка от неточности нулей
  E6 = zero_acc * sum(j=1, nzeros,
    my(g=zeros[j]);
    sqrt(x)/sqrt(1/4+g^2)*exp(lam^2*(1/4-g^2)/2)
  );

  Etotal = E1 + E2 + E3 + E4 + E5 + E6;

  print("Составляющие погрешности (хвост нулей может быть больше в 10 раз!):");
  print("  1. Хвост нулей (оптимистично): ", E1);
  print("  2. Хвост главн. члена:         ", E2);
  print("  3. I_{-1}:                     ", E3);
  print("  4. Хвост решета:               ", E4);
  print("  5. Пропущенные нули:           ", E5);
  print("  6. Неточность нулей:           ", E6);
  print("  Суммарная:                     ", Etotal);
  if(Etotal < 0.5,
    print("  OK: оптимистично погрешность < 0.5, можно округлять"),
    print("  ВНИМАНИЕ: погрешность >= 0.5")
  );
}

Запуск:

? platt_error_estimates(10^11, 55000, 4e-9, 1e-4, 1e2);
Составляющие погрешности (хвост нулей может быть больше в 10 раз!):
  1. Хвост нулей (оптимистично): 0.31893411871298175845611926969408738343
  2. Хвост главн. члена:         6.2499741402053998968848197263833762216 E-10
  3. I_{-1}:                     0.00031834163459346596101071978766488784684
  4. Хвост решета:               0.16602200777641179714083432175110640678
  5. Пропущенные нули:           0
  6. Неточность нулей:           0.0055770484800786496157555152504386738177
  Суммарная:                     0.49085151722906308519425981617177932450
OK: оптимистично погрешность < 0.5, можно округлять
time = 391 ms.
?

[wrest]

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

Не согласен с таким обобщением. Вот например для [0, 60] интересна зона вблизи 106-го триллиона, то бишь порядок $10^{14}$.

Там вам были нужны интегральные характеристики (плотность на больших промежутках, т.е. сглаженная), формула (на основе HL1) для этого есть и она вполне нормальная.

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

Dmitriy40
Смотрите: погрешность нулей в 4е-9 на количестве 5e4 нулей не вносит существенной ошибки, но если бы нули были с ошибкой 1e-7 то уже бы было заметно.
Так что если считать более высокие интервалы, то и нули надо брать поточнее, но до интервалов 1e12..1e13 хватит и Одлыжко.

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

--
А вот и сенсация

? platt_error_estimates(10^12, 90000, 4e-9, 7e-5, 1e2);
Составляющие погрешности (хвост нулей может быть больше в 10 раз!):
  1. Хвост нулей (оптимистично): 0.16256374140121178472860806733846929054
  2. Хвост главн. члена:         4.0104000630719447474628860462499721325 E-9
  3. I_{-1}:                     6.4965736680873739546058682364574415895 E-5
  4. Хвост решета:               1.0653078832319756983203535645695994435
  5. Пропущенные нули:           0
  6. Неточность нулей:           0.019363884437598127502418740027510299810
  Суммарная:                     1.2473004788178665473628711780808296545
  ВНИМАНИЕ: погрешность >= 0.5
time = 554 ms.
? platt_pi(10^12,90000,0.00007,100)
Получаем 90000 нулей ... готово за 133 мс
Вычисляем главный член ... готово за 48 мс
Вычисляем оставшийся ряд phihat ... готово за 18005 мс
Суммируем ряд ... готово за 14 мс
Вычисляем поправку решетом ... готово за 43033 мс
Обращение Мёбиуса ... готово за 1 мс
time = 1min, 931 ms.
37607912018.452964917722339618161031514
? primepi(10^12)
time = 1min, 58,726 ms.
37607912018
?

На предвычисленных нулях Платт почти в джва раза обгоняет встроенную в pari/gp функцию primepi()
Хотя, конечно, и "на тоненького".

Неожиданно, однако.

Но, конечно, спецрешения на всё это смотрят свысока

~/primecount/build $ time ./primecount 1e12
37607912018

real    0m0.048s
user    0m0.088s
sys     0m0.050s
~/primecount/build $

Попутно, нашёлся ещё один баг в оценке ошибки (теперь в оценке хвоста решета).
Исправленная версия platt_error_estimates(x, nzeros, zero_acc, lam, nseg)

(Оффтоп)

platt_error_estimates(x, nzeros, zero_acc, lam, nseg) = {
\\ Оценка составляющих погрешности и суммарной
  my(v = readvec("zeros1"));
  if(#v < nzeros, error("Не хватает нулей в файле ", file," Надо: ",n," Есть:",#v));
  my(zeros=Vec(v)[1..nzeros]);
  my(c=5, T1=zeros[nzeros], T=8/lam);
  my(E1, E2, E3, E4, E5, E6, Etotal);

  \\ 1. Хвост нулей (Лемма A.4, консервативная оценка)
  \\ E1 = 2*exp(lam^2*(1-T1^2)/2) * (sqrt(x)/(T1*log(x)) + 1/(lam^2*T1^2*x)) * (lam^2*T1^2+1);
  \\ 1. Хвост нулей (численное интегрирование, всегда корректно)
  E1 = 2*intnum(t=T1, 8/lam, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
  \\ В print_error_estimates, после вычисления E1:
  E1 = E1 / 10;  \\ эмпирический коэффициент сокращения осцилляций

  \\ 2. Хвост главного члена (Лемма A.2)
  E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));

  \\ 3. I_{-1} (Лемма 4.5)
  E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);

  \\ 4. Хвост решета (вне окна [x*e^{-c*lam}, x*e^{c*lam}])

  \\ Было (завышено в c/sqrt(2) раз):
  \\E4 = x*lam/log(x) * exp(-c^2/2)/(sqrt(Pi)*c);

  \\ Правильно:
  E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);

  \\ 5. Проверка нулей по формуле N(T)
  my(N_exp = T1/(2*Pi)*log(T1/(2*Pi)) - T1/(2*Pi) + 7/8);
  my(bound = 0.137*log(T1) + 0.443*log(log(T1)) + 1.588);
  my(diff = abs(N_exp - nzeros));
  E5 = if(diff > bound, diff, 0);

  \\ 6. Ошибка от неточности нулей
  E6 = zero_acc * sum(j=1, nzeros,
    my(g=zeros[j]);
    sqrt(x)/sqrt(1/4+g^2)*exp(lam^2*(1/4-g^2)/2)
  );

  Etotal = E1 + E2 + E3 + E4 + E5 + E6;

  print("Составляющие погрешности (хвост нулей может быть больше в 10 раз!):");
  print("  1. Хвост нулей (оптимистично): ", E1);
  print("  2. Хвост главн. члена:         ", E2);
  print("  3. I_{-1}:                     ", E3);
  print("  4. Хвост решета:               ", E4);
  print("  5. Пропущенные нули:           ", E5);
  print("  6. Неточность нулей:           ", E6);
  print("  Суммарная:                     ", Etotal);
  if(Etotal < 0.5,
    print("  OK: оптимистично погрешность < 0.5, можно округлять"),
    print("  ВНИМАНИЕ: погрешность >= 0.5")
  );
}

[Dmitriy40]

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

На предвычисленных нулях Платт почти в джва раза обгоняет встроенную в pari/gp функцию primepi()

Так пишете же что погрешность 1.2473, т.е .вычисление неправильное!
Ну и округление 0.453 до 0 не слишком надёжное, на грани, для большей уверенности в результате хотелось бы погрешности поменьше. А значит вероятно надо больше нулей - и дольше.

[wrest]

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

Так пишете же что погрешность 1.2473, т.е .вычисление неправильное!

Ага, это баг в оценке погрешности был.
После исправления:

? platt_error_estimates(10^12, 90000, 4e-9, 7e-5, 1e2);
Составляющие погрешности (хвост нулей может быть больше в 10 раз!):
  1. Хвост нулей (оптимистично): 0.16256374140121178472860806733846929054
  2. Хвост главн. члена:         4.0104000630719447474628860462499721325 E-9
  3. I_{-1}:                     6.4965736680873739546058682364574415895 E-5
  4. Хвост решета:               0.30131457131392670426046446354956253364
  5. Пропущенные нули:           0
  6. Неточность нулей:           0.019363884437598127502418740027510299810
  Суммарная:                     0.48330716689981755330298207706079274466
OK: оптимистично погрешность < 0.5, можно округлять
time = 561 ms.
?

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

Ну и округление 0.453 до 0 не слишком надёжное,

Ну да, я ж и написал что "на тоненького"

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

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

А значит вероятно надо больше нулей - и дольше.

Да, дольше, но всё ещё быстрее primepi()

Но вывод тут простой: это не Платт такой хороший, это встроенная primepi() неоптимизирована/не ускорена.
Что разрабы pari/gp и подтверждают в доках:

Цитата:

Uses checkpointing and a naive O(x) algorithm; make sure to start gp with primelimit at least sqrt{x}.

Но вот и рашёлся предел по нулям zeros1 Одлыжко: интервал 1e12
Функция которая считает сколькот надо нулей при данной $\lambda$

needed_zeros(lam)=my(T=6/lam);round(T/(2*Pi)*log(T/(2*Pi*exp(1)))+7/8)

Запуск:

? needed_zeros(7e-5)
116242
?

Это оценка сверху, но натпрактике оказалось что и 90000 достаточно.
Но уже при $\lambda = 8\cdot 10^{-5}$ имеем

? needed_zeros(8e-5)
100118
?

Тютелька-в-тютельку.

Ну и при 100000 нулей получаем:

? platt_pi(10^12,100000,8e-5,100)
Получаем 100000 нулей ... готово за 138 мс
Вычисляем главный член ... готово за 50 мс
Вычисляем оставшийся ряд phihat ... готово за 20131 мс
Суммируем ряд ... готово за 17 мс
Вычисляем поправку решетом ... готово за 49721 мс
Обращение Мёбиуса ... готово за 1 мс
time = 1min, 9,689 ms.
37607912017.932542925046422344791422919?

1 минута 10 сек, ошибка 0.07
А primepi() две минуты.

[Dmitriy40]

Так primepi() надо сравнивать не с Платтом и не с primecount, а с primesieve, именно решето Эратосфена:

? primepi(10^12)
time = 2min, 40,167 ms.
%1 = 37607912018

>primesieve 1e12 -t 1 -c
Sieve size = 128 KiB
Threads = 1
100%
Seconds: 254.634
Primes: 37607912018

Выигрывает однако, хотя primesieve тоже очень хорошо оптимизирована.

[wrest]

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

Выигрывает однако, хотя primesieve тоже очень хорошо оптимизирована.

Ну что тут сказать... Наверное, что простое решето Эратосфена -- это такое себе.

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

Я думаю ради интереса выдавить из Квена скрипт для комбинаторного primepi средствами pari/gp, возможно готового к компиляции в gp2c. Такого же выигрыша как у primecount (три порядка на 10^12) конечно не ожидаю, но пара порядков я думаю будет
Вообще это тема... Алгоритмы общего назначения же не патентуемая штука, вроде, разрабы pari/gp вполне могли бы улучшить то из базового, что у них работает медленно, при помощи ИИ, чтобы не затрачивать чрезмерных усилий. Возможно, они это и делают уже

[Yadryara]

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

1 минута 10 сек, ошибка 0.07
А primepi() две минуты.

Ну это как раз не сенсация, я как раз надеялся что при правильном подборе параметров аналитическо-комбинаторный метод обгонит комбинаторные primepi и primesieve. Надеюсь, что повыше, например для $10^{13}$ выигрыш будет ещё больше.

А вот сумасшедшая скорость PrimeCount для меня радостный сюрприз. И если уже давно такая скорость, то почему например, в проекте SPT использовалась primesieve, а не PrimeCount?