|
Последний раз редактировалось wrest 17.08.2026, 02:33, всего редактировалось 15 раз(а).
platt_v7_par.gp ВЕРСИЯ 7 ПараллельнаяCodebase (набор скриптов), скачать в файл например platt_v7_par.gpГлобальных переменных нет, всё в функциях. (Оффтоп)
Код: \\ ============================================================
\\ platt_v7_par.gp — аналитическое pi(x) методом Платта
\\ (явная формула Римана + гауссово сглаживание + обращение Мёбиуса)
\\ Статья на основе которой работает алгоритм
\\ David J. Platt Computing π(x) Analytically
\\ https://arxiv.org/abs/1203.5712
\\ https://doi.org/10.48550/arXiv.1203.5712
\\
\\ ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ
\\
\\ Нули берутся из таблицы Одлыжко (zeros1 = 10^5 нулей,
\\ zeros6 = 10^6 нулей). Нулевая часть — рядами (лемма 6.1),
\\ решётная поправка — сегментами с разложением Тейлора (§5).
\\
\\ БЫСТРЫЙ СТАРТ
\\ ---------------
\\ \r platt_v7_par.gp
\\ \p 50
\\ platt_pi(10^8, 100, 0.04, 100,5,0) \\ -> 5761455.00...
\\ primepi(10^8) \\ сверка: 5761455
\\
\\ ИНТЕРФЕЙС
\\ ---------
\\ platt_pi(x, nzeros, lam, nseg,prime_sieve,verbose)
\\ x — аргумент pi(x);
\\ nzeros — сколько нулей дзета использовать;
\\ lam — параметр сглаживания;
\\ nseg — число сегментов Тейлора в решётной поправке, рекомендуется 100.
\\ prime_sieve — ширина решета в конце интервала для поправки, от 4 до 7
\\ verbose — 1 печатать статистику, 0 просто вернуть результат
\\
\\ ПАРАМЕТРЫ (проверенные ориентиры)
\\ ---------------------------------
\\ Правило: нулей нужно до высоты ~6/lam; lam балансирует
\\ длину интеграла главного члена (~1/lam) и ширину решета (~lam).
\\ Функция platt_opt(x) проводит измерение производительности
\\ и на его основе выдаёт оптимальные по времен вычислений параметры
\\
\\ Точность: погрешность << 0.5, округление platt_pi даёт точное целое.
\\ ============================================================
\\ ================= НУЛИ ДЗЕТА-ФУНКЦИИ =================
\\ должны быть в файле zeros1 (0^5 шт.) и zeros6 (2*10^6 шт.)
\\ например нули из таблиц Одлыжко
\\ отсюда: https://www-users.cse.umn.edu/~odlyzko/zeta_tables/index.html
\\ ================= НУЛИ ДЗЕТА-ФУНКЦИИ =================
load_zeros(n) = {
my(filename=if(n>10^5,"zeros6", "zeros1"));
my(fd, vals, val, i);
fd = fileopen(filename, "r");
vals = vector(n);
for(i = 1, n,
val = eval(fileread(fd));
if(val == 0, fileclose(fd);error("Не хватает нулей в файле ", filename," Надо: ",n," Есть:",i));
vals[i] = val;
);
fileclose(fd);
return(vals);
}
\\ ============ СГЛАЖЕННАЯ НУЛЕВАЯ ЧАСТЬ (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);
\\ Предвычисление степеней один раз на вызов
my(pow_s0 = vector(K+1), pow_lam = vector(K\2+1));
my(inv_s0 = -I/s0, half_lam2 = -lam^2/2);
pow_s0[1] = 1;
for(n=1, K, pow_s0[n+1] = pow_s0[n]*inv_s0);
pow_lam[1] = 1;
for(m=1, K\2, pow_lam[m+1] = pow_lam[m]*half_lam2/m);
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 += pow_s0[n+1] * pow_lam[m+1];
);
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);
if(n < 2, return(vector(n)));
\\ Экспорт функций для параллельного использования
export(step_incr);
export(phihat);
\\ Параллельное вычисление всех шагов
my(steps = parvector(n-1, j, step_incr(1/2 + I*zeros[j], zeros[j+1] - zeros[j], x, lam, K)));
\\ Последовательное накопление
my(G = 0, GatZ = vector(n));
GatZ[1] = 0;
for(j=2, n, G += steps[j-1]; 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 Платта) ============
\\ Вспомогательная функция для обработки одного сегмента
\\ Все параметры передаются явно (не через нужно для параллельности)
process_segment(j, a, b, x0, jx, x, lo, sq2lam, sqrtpi, invx, w) = {
my(N=0, S1=0, S2=0, C=0, d, g, eg, ph0, ph1, ph2, ph3);
if(b >= a,
if(j < jx,
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; N++);
C=N;
, if(j > jx,
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; N++);
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);
ph3=eg/(sqrtpi*sq2lam*x0^3)*(2/sq2lam^2 - 2 - (2*g/sq2lam)^2 - 3*(2*g/sq2lam));
ph1 += ph3*w^2/8;
C - (ph0*N - ph1*S1 + ph2*S2/2);
, 0)
}
\\ sum_{p^m} (chi-phi)/m; сегменты + Тейлор для phi (erfc один раз в сегмент)
\\ Ширина окна = 2*c*lam*x, c=5 рекомендуется (проверено: c=5,8,10 дают одно и то же)
\\ ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ
sieve_corr_taylor(x, lam, nseg, c=5) = {
my(lo=x*exp(-c*lam), hi=x*exp(c*lam));
my(sq2lam=sqrt(2)*lam, sqrtpi=sqrt(Pi), invx=1/x);
my(w=(hi-lo)/nseg);
my(B=vector(nseg+1, j, round(lo + (j-1)*(hi-lo)/nseg)));
\\ Индекс сегмента, содержащего границу x
my(jx=1);
while(jx<nseg && B[jx+1]<=x, jx++);
\\ Экспорт функции для параллельного использования
export(process_segment);
\\ Параллельная обработка всех сегментов
my(s = vecsum(parvector(nseg, j,
my(a=B[j], b=B[j+1]-1, x0);
if(b >= a, x0=(a+b)/2, x0=0);
process_segment(j, a, b, x0, jx, x, lo, sq2lam, sqrtpi, invx, w)
)));
\\ Степени простых (последовательно, их мало)
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);
\\ ============ СБОРКА ============
platt_pi(x, nzeros, lam, nseg=100, c=5, verbose=2) = {
if(verbose>1,
platt_error_estimates(x, nzeros, 4*10^(-9), lam, nseg, c);
);
my(s,res,t0);
if(verbose,
t0=getwalltime();
print1("Получаем ", nzeros, " нулей ...");
);
my(zeros = load_zeros(nzeros));
if(verbose,
print(" готово за ", getwalltime()-t0, " мс");
t0=getwalltime();
print1("Вычисляем главный член ... ");
);
my(re1 = Phihat_re_at_1(x, lam, 0.5, 40));
if(verbose,
print(strprintf("%.3f",re1), " готово за ", getwalltime()-t0, " мс");
t0=getwalltime();
print1("Вычисляем остаток ряда phihat ... ");
);
\\ ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ: Phihat_re_all теперь использует parvector
my(rez = Phihat_re_all(x, lam, zeros, 20));
my(zsum = sum(j=1, #rez, rez[j]));
if(verbose,
print(strprintf("%.3f", -2*zsum)," готово за ", getwalltime()-t0, " мс");
t0=getwalltime();
print1("Вычисляем поправку решетом ... ");
);
\\ ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ: sieve_corr_taylor теперь использует parvector
res=sieve_corr_taylor(x, lam, nseg, c);
if(verbose,
print(strprintf("%.3f",res)," длина ",strprintf("%.4g",x*(exp(c*lam)-exp(-c*lam)))," готово за ", getwalltime()-t0, " мс");
);
res=re1 - 2*zsum - log(2) + res;
if(verbose,
print("Сумма Phi ... ",strprintf("%.3f",res));
t0=getwalltime();
print1("Поправка Мёбиуса ... ");
);
\\ ПАРАЛЛЕЛЬНАЯ ВЕРСИЯ: поправки Мёбиуса
export(ftarget);
my(nmax = floor(log(x)/log(2)));
s = vecsum(parvector(nmax-1, n,
my(nn = n+1);
if(moebius(nn)!=0, (moebius(nn)/nn)*ftarget(sqrtn(x,nn)), 0)
));
if(verbose,
print(strprintf("%.3f",s)," готово за ", getwalltime()-t0, " мс");
);
return(res+s);
}
\\ Оценка погрешностей
platt_error_estimates(x, nzeros, zero_acc, lam, nseg, c) = {
my(zeros = load_zeros(nzeros));
my( T1=zeros[nzeros], T=8/lam);
my(E1, E2, E3, E4, E5, E6, Etotal);
\\ Эмпирический коэффициент, в надежде что нули осциллируют
my(kE1 = 1/4);
E1 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));
E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
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);
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("Составляющие погрешности (хвост нулей может быть больше!):");
printf(" 1. Хвост нулей (оптимистично): %.3g\n", E1);
printf(" 2. Хвост главн. члена: %.3g\n", E2);
printf(" 3. I_{-1}: %.3g\n", E3);
print(" 4. Размер решета: ",strprintf("%.4g",x*(exp(c*lam)-exp(-c*lam))));
printf(" 5. Хвост решета: %.3g\n", E4);
print(" 6. Пропущенные нули: ", E5);
printf(" 7. Неточность нулей: %.3g\n", E6);
printf(" Суммарная: %.3g\n", Etotal);
if(Etotal < 0.5,
print(" OK: оптимистично погрешность < 0.5, можно вычислять"),
print(" ВНИМАНИЕ: погрешность > 0.5, вероятна ошибка")
);
}
\\ ===== Подбор оптимальных параметров с измерением реальными функциями =====
platt_opt(x) = {
my(E_max=0.15, kE1=1/4);
my(lam, T1, nzeros, c, T);
my(E1, E2, E3, E4, Etotal, budget, budget_E4);
my(n_primes, time_total, i, best_time, best_lam, best_nzeros, best_c, best_T1);
my(lam_min, lam_max, lam_step);
my(cost_prime, cost_zero, cost_load);
my(N_max = 10^6); \\ максимальное доступное число нулей (zeros6)
\\ === Измерение стоимостей реальными функциями ===
\\ 1. Стоимость простого: запуск sieve_corr_taylor на окне шириной 10^8
my(lam_test = 1e7/x);
my(nseg_test = 500);
my(t0 = getwalltime());
sieve_corr_taylor(x, lam_test, nseg_test, 5);
my(time_sieve = (getwalltime()-t0)*1e-3);
my(width_test = x*(exp(5*lam_test)-exp(-5*lam_test)));
my(n_primes_test = width_test/log(x));
cost_prime = time_sieve / n_primes_test;
\\ 2. Стоимость загрузки нулей
my(n_test = 5000);
n_test = min(n_test, 10^6);
t0 = getwalltime();
my(zeros = load_zeros(n_test));
cost_load = (getwalltime()-t0)/n_test*1e-3;
\\ 3. Стоимость вычисления нуля: запуск Phihat_re_all
my(test_lam = 1e-4);
t0 = getwalltime();
Phihat_re_all(x, test_lam, zeros, 20);
cost_zero = (getwalltime()-t0)/n_test*1e-3;
printf("Измеренные стоимости:\n");
printf(" Простое (решето): %.2f мкс\n", cost_prime*1e6);
printf(" Загрузка нуля: %.2f мкс\n", cost_load*1e6);
printf(" Вычисление нуля: %.2f мкс\n", cost_zero*1e6);
\\ === Начальное приближение lam ===
lam = 1e-4;
for(i=1, 10,
lam = sqrt(3*cost_zero*log(x)/(2*Pi*cost_prime*5*x) * log(6/(2*Pi*lam*exp(1))));
);
lam_min = lam/3;
lam_max = lam*3;
lam_step = (lam_max - lam_min)/20;
best_time = 1e30;
best_lam = 0; best_nzeros = 0; best_c = 0; best_T1 = 0;
\\ === Перебор lam ===
for(k=0, 20,
lam = lam_min + k*lam_step;
T = 8/lam;
\\ Фиксированные погрешности для данного lam
E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));
E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);
\\ Бюджет погрешности для E1 + E4
budget = E_max - E2 - E3;
if(budget <= 0, next);
\\ Перебор T1 (высота нулей)
my(T1_min = 4/lam, T1_max = T);
my(n_T1 = 10);
for(j=0, n_T1,
T1 = T1_min + j*(T1_max - T1_min)/n_T1;
\\ Ограничение на число нулей
nzeros = round(T1/(2*Pi)*log(T1/(2*Pi*exp(1))) + 7/8);
\\ if(nzeros > N_max, next);
\\ Погрешность хвоста нулей
E1 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
if(E1 >= budget, next);
\\ Оставшийся бюджет для E4
budget_E4 = budget - E1;
\\ Подбор c: минимальное c такое, что E4(c) <= budget_E4
c = 3;
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
while(E4 > budget_E4 && c < 15,
c += 0.25;
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
);
if(E4 > budget_E4, next);
\\ Время
n_primes = 2*c*lam*x/log(x);
time_total = (cost_load + cost_zero)*nzeros + cost_prime*n_primes + 0.02;
if(time_total < best_time,
best_time = time_total;
best_lam = lam;
best_nzeros = nzeros;
best_c = c;
best_T1 = T1;
);
);
);
if(best_time >= 1e30,
print("ВНИМАНИЕ: не найдены параметры с погрешностью < ", E_max);
return([0, 0, 0, 0, 0]);
);
\\ === Вывод результата ===
lam = best_lam; nzeros = best_nzeros; c = best_c; T1 = best_T1;
T = 8/lam;
E1 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));
E3 = exp(lam^2/2)/(2*Pi*x*lam) * (5/sqrt(2*Pi) + 2/lam);
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
Etotal = E1 + E2 + E3 + E4;
n_primes = 2*c*lam*x/log(x);
time_total = (cost_load + cost_zero)*nzeros + cost_prime*n_primes + 0.02;
printf("\nОптимальные параметры для x=%.3g (по времени счёта):\n", x);
printf(" lam = %.4g\n", lam);
printf(" nzeros = %.4g = (%d)\n", nzeros, nzeros);
printf(" c = %.2f\n", c);
printf(" sieve = %.4g\n", x*(exp(c*lam)-exp(-c*lam)));
printf("\nОценка погрешности:\n");
printf(" Хвост нулей: %.3g\n", E1);
printf(" Хвост главн. члена: %.3g\n", E2);
printf(" I_{-1}: %.3g\n", E3);
printf(" Хвост решета: %.3g\n", E4);
printf(" Суммарная: %.3g\n", Etotal);
if(Etotal < 0.5,
print(" OK: погрешность < 0.5"),
print(" ВНИМАНИЕ: погрешность >= 0.5")
);
printf("\nПрогноз времени:\n");
printf(" Загрузка нулей: %.1f с\n", 1.2*cost_load*nzeros);
printf(" Ряд phihat: %.1f с\n", 1.2*cost_zero*nzeros);
printf(" Решето: %.1f с\n", 1.2*cost_prime*n_primes);
printf(" Итого: %.1f с\n", 1.2*time_total);
print("Вызов: ");
print("platt_pi(",x,",",nzeros,",",strprintf("%.4g",lam),",",500,",",strprintf("%.2f",c),",1)");
print();
} Запуск на ноуте: (Оффтоп)
Код: ? platt_opt(10^14)
Измеренные стоимости:
Простое (решето): 0.34 мкс
Загрузка нуля: 3.40 мкс
Вычисление нуля: 91.00 мкс
Оптимальные параметры для x=1.00 e14 (по времени счёта):
lam = 8.085 e-6
nzeros = 1.134 e6 = (1133595)
c = 5.5
sieve = 8.894 e9
Оценка погрешности:
Хвост нулей: 0.101
Хвост главн. члена: 3.97 e-8
I_{-1}: 4.87 e-5
Хвост решета: 0.179
Суммарная: 0.279
OK: погрешность < 0.5
Прогноз времени:
Загрузка нулей: 3.9 с
Ряд phihat: 103.2 с
Решето: 93.6 с
Итого: 200.6 с
Вызов:
platt_pi(100000000000000,1133595,8.085 e-6,500,5.5,1)
cpu time = 10,676 ms, real time = 1,701 ms.
? platt_pi(100000000000000,1133595,8.085 e-6,500,5.5,1)
Получаем 1133595 нулей ... готово за 3167 мс
Вычисляем главный член ... 3204942065789.476 готово за 14 мс
Вычисляем остаток ряда phihat ... 18188.361 готово за 121628 мс
Вычисляем поправку решетом ... 867.002 длина 8.894 e9 готово за 126806 мс
Сумма Phi ... 3204942084844.146
Поправка Мёбиуса ... -334042.484 готово за 3 мс
cpu time = 27min, 46,689 ms, real time = 4min, 11,643 ms.
3204941750801.6614053754253971336487334
? Коэффициент параллельности очень хороший cpu time = 27min, 46,689 ms, real time = 4min, 11,643 ms. (у меня 8 потоков) 10:14 считается за 4 минуты, но это ноутбук с тротлингом и вот этим всем. Некоторые функции разделены на параллельную и последовательную части выделением новых функций наружу. Поправлены некоторые коэффициенты, влияющие на точность. Переделан бенчмарк производительности и подбор параметров: две функции сведены в одну platt_opt(), и теперь должно попадать в погрешность 0.5 всегда. Ну или по крайней мере до 10^14  Бенчмарк запускает реальные функции и запускает параллельно. Ошибка бенчмарка по времени выполнения у меня получатся 10-30%, думаю норм. При verbose=0 в основной функции platt_pi() возвращается только результат, verbose=1 печатается статистика и тайминги, verbose=2 дополнительно печатается оценка погрешностей до вычислений. Поскольку ноутбук беру в руки редко, а на планшете нет параллельного pari/gp, вряд ли буду что-то обновлять в ближайшее время.
|
![$$\pi^*(x) = \frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} \log \zeta(s) \hat{\varphi}(s) ds + \sum_{p^m} \frac{1}{m} [\chi_x(p^m) - \varphi(p^m)] \qquad\qquad\qquad \text{(3.2)}$$](https://dxdy.ru/math/fd854c27adc4361b802e5ade3d94cbda82.png "$$\pi^*(x) = \frac{1}{2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} \log \zeta(s) \hat{\varphi}(s) ds + \sum_{p^m} \frac{1}{m} [\chi_x(p^m) - \varphi(p^m)] \qquad\qquad\qquad \text{(3.2)}$$")
