|
Последний раз редактировалось wrest 16.08.2026, 01:22, всего редактировалось 6 раз(а).
Новая версия ( platt_v5.gp) набора скриптов подсчёта  методом Платта Немного ускоренная, немного дополненная по математике вычислений. Главное обновление: добавлена функция platt_opt(x) которая измеряет скорость вычислений на конкретном компе и исходя из результатов рекомендует параметры для запуска главной функции platt_pi() исходя из минимизации времени вычислений путём поиска таких параметров что время на решето и на интегрирование одинаковое и погрешность не превысит 0.5. При расчёте учитывается размер применимости встроенной таблицы простых и замедление forprime за его границами. В главной функции добавлено два параметра - широта решета c и отключение печати статистики (толькотвозврао результата). Количество сегментов для разложения интеграла вдоль нетривиальных нулей в ряд Тейлора оставлено, но лучше его не менять. Почти не влияет на результаттначиная со 100 сегментов. (Оффтоп)
Код: \\ ============================================================
\\ platt.gp — аналитическое pi(x) методом Платта
\\ (явная формула Римана + гауссово сглаживание + обращение Мёбиуса)
\\
\\ Нули берутся из таблицы Одлыжко (zeros1 = 10^5 нулей,
\\ zeros6 = 10^6 нулей). Нулевая часть — рядами (лемма 6.1),
\\ решётная поправка — сегментами с разложением Тейлора (§5).
\\
\\ БЫСТРЫЙ СТАРТ
\\ ---------------
\\ \r platt_v5.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 просто вернуть результат
\\
\\ ПАРАМЕТРЫ (проверенные ориентиры)
\\ ---------------------------------
\\ x ~ 10^6 : nzeros=12, lam=0.10, nseg=100
\\ x ~ 10^8 : nzeros=100, lam=0.04, nseg=100
\\ x ~ 10^10: nzeros=8000, lam=0.001, nseg=100
\\ Правило: нулей нужно до высоты ~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, 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, 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)));
my(s=0, j, a, b, x0, N, S1, S2, C, d, g, eg, ph0, ph1, ph2, ph3, p);
\\ Индекс сегмента, содержащего границу x
my(jx=1);
while(jx<nseg && B[jx+1]<=x, jx++);
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;
if(j < jx,
\\ сегмент целиком ниже x: все простые < x
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; N++);
C=N;
, if(j > jx,
\\ сегмент целиком выше x: ни один простой не < x
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; N++);
C=0;
, \\ сегмент содержит x: единственный случай с if
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;
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);
\\ ============ СБОРКА ============
platt_pi(x, nzeros, lam, nseg=100, c=5, verbose=1) = {
if(verbose,
\\ оценка и печать ошибки при выбранных параметрах
platt_error_estimates(x, nzeros, 4*10^(-9), lam, nseg, c);
);
\\ вычисление сглаженной f(x)
my(s,res,t0);
if(verbose,
t0=getwalltime();
print1("Получаем ", nzeros, " нулей ...");
);
\\ загружаем нули
my(zeros = load_zeros(nzeros));
if(verbose,
print(" готово за ", getwalltime()-t0, " мс");
t0=getwalltime();
print1("Вычисляем главный член ... ");
);
\\ главный член, вклад полюса при s=1
my(re1 = Phihat_re_at_1(x, lam, 0.5, 40));
if(verbose,
print(strprintf("%.3f",re1), " готово за ", getwalltime()-t0, " мс");
t0=getwalltime();
print1("Вычисляем остаток ряда phihat ... ");
);
\\ вклад нулей дзета функции
my(rez = Phihat_re_all(x, lam, zeros, 12));
my(zsum = sum(j=1, #rez, rez[j]));
if(verbose,
print(strprintf("%.3f", -2*zsum)," (готово за ", getwalltime()-t0, " мс)");
t0=getwalltime();
print1("Вычисляем поправку решетом ... ");
);
\\ вклад простых на конце интервала
res=sieve_corr_taylor(x, lam, nseg, c);
if(verbose,
print(strprintf("%.3f",res)," (старт: ",strprintf("%.4g",x*exp(-c*lam))," длина: ",strprintf("%.4g",x*(exp(c*lam)-exp(-c*lam)))," ) (готово за ", getwalltime()-t0, " мс)");
);
\\ общий итог f(x)
res=re1 - 2*zsum - log(2) + res;
if(verbose,
print("Сумма ... ",strprintf("%.3f",res));
t0=getwalltime();
);
\\ pi(x) = f(x) + поправки Мёбиуса при n>=2
s=0;
if(verbose,
print1("Поправки Мёбиуса ... ");
);
for(n=2, floor(log(x)/log(2)),
if(moebius(n)!=0, s += (moebius(n)/n)*ftarget(sqrtn(x,n)));
);
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/5);
\\ 1. Хвост нулей (Лемма A.4, консервативная оценка)
\\ Оценка, может быть с ошибкой
\\ E1 = kE1*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 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
\\ В print_error_estimates, после вычисления E1:
\\ 2. Хвост главного члена (Лемма A.2)
E2 = exp(lam^2*(1-T^2)/2) * (x/(T*log(x)) + 1/(lam^2*T^2*x));
\\ 3. Интеграл вокруг -1 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*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("Составляющие погрешности (хвост нулей может быть больше!):");
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, вероятна ошибка")
);
}
\\ ===== Измерение стоимостей с использованием реальных функций =====
measure_costs(x, nseg) = {
my(primelimit = default(primelimit));
my(threshold = primelimit^2);
my(width = 10^6);
my(cost_fast, cost_mid, cost_slow, seg_factor);
my(cost_zero, cost_load);
my(t, cnt, a, b, d, S1, S2, x0);
my(n_test = 2000, test_lam = 1e-4, zeros, sum_test, t0);
\\ Ограничение на число нулей (zeros6 содержит 10^6)
n_test = min(n_test, 10^6);
\\ === Стоимость простых: быстрый режим (до порога таблицы) ===
x0 = primelimit; a = x0 - width\2; b = x0 + width\2;
t = gettime(); cnt = 0; S1 = 0; S2 = 0;
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; cnt++);
t = gettime() - t;
cost_fast = if(cnt>0 && t>=0, t/cnt*1e-3, 1.3e-6);
\\ === Стоимость простых: средний режим (за порогом таблицы, до 2^63) ===
if(threshold * 10 < 2^63,
x0 = threshold * 10; a = x0 - width\2; b = x0 + width\2;
t = gettime(); cnt = 0; S1 = 0; S2 = 0;
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; cnt++);
t = gettime() - t;
cost_mid = if(cnt>0 && t>=0, t/cnt*1e-3, cost_fast*10);
,
cost_mid = cost_fast * 10;
);
\\ === Стоимость простых: медленный режим (за порогом 2^63) ===
if(2^63 * 2 < 1e30,
x0 = 2^63 * 2; a = x0 - width\2; b = x0 + width\2;
t = gettime(); cnt = 0; S1 = 0; S2 = 0;
forprime(p=a, b, d=x0-p; S1+=d; S2+=d*d; cnt++);
t = gettime() - t;
cost_slow = if(cnt>0 && t>=0, t/cnt*1e-3, cost_mid*5);
,
cost_slow = cost_mid * 5;
);
\\ Поправка на сегментацию
seg_factor = 1 + nseg/1000;
cost_fast *= seg_factor;
cost_mid *= seg_factor;
cost_slow *= seg_factor;
\\ === Стоимость нулей: используем реальные функции ===
\\ Загрузка нулей
t0 = gettime();
zeros = load_zeros(n_test);
cost_load = (gettime()-t0)/n_test*1e-3;
\\ Вычисление phihat реальной функцией
t0 = gettime();
sum_test = Phihat_re_all(x, test_lam, zeros, 12);
cost_zero = (gettime()-t0)/n_test*1e-3;
printf("Измеренные стоимости:\n");
printf(" primelimit = %d, порог таблицы = %.3g, порог 64 бит = %.3g\n",
primelimit, threshold, 2^63);
printf(" Простые (быстрый, < %.3g): %.2f мкс\n", threshold, cost_fast*1e6);
printf(" Простые (средний, < 2^63): %.2f мкс\n", cost_mid*1e6);
printf(" Простые (медленный, > 2^63): %.2f мкс\n", cost_slow*1e6);
printf(" Загрузка нуля: %.2f мкс\n", cost_load*1e6);
printf(" Вычисление нуля (phihat): %.2f мкс\n", cost_zero*1e6);
return([threshold, cost_fast, cost_mid, cost_slow, cost_zero, cost_load]);
}
\\ ===== Подбор оптимальных параметров для метода Платта =====
platt_opt(x) = {
my(target=0.2, kE1=1/5);
my(lam, T1, nzeros, c, T);
my(E1, E2, E3, E4, Etotal);
my(n_primes, time_total, i, best_time, best_lam, best_nzeros, best_c, best_T1);
my(lam_min, lam_max, lam_step, T1_new, E1_new, c_new, E4_new);
my(cost_prime, upper);
my(costs = measure_costs(x, 500));
my(threshold=costs[1], cost_fast=costs[2], cost_mid=costs[3], cost_slow=costs[4]);
my(cost_zero=costs[5], cost_load=costs[6]);
\\ Начальное приближение lam из баланса времени
cost_prime = cost_fast;
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;
for(k=0, 20,
lam = lam_min + k*lam_step;
T = 8/lam;
\\ Подбор T1 (высота нулей)
T1 = 6/lam;
E1 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
while(E1 > target && T1 < T,
T1 = min(T1*1.15, T);
E1 = kE1*2*intnum(t=T1, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
);
while(E1 < target/3 && T1 > 2/lam,
T1_new = T1/1.05;
E1_new = kE1*2*intnum(t=T1_new, T, sqrt(x)/t * exp(-lam^2*t^2/2) * log(t)/(2*Pi));
if(E1_new < target, T1 = T1_new; E1 = E1_new, break);
);
nzeros = round(T1/(2*Pi)*log(T1/(2*Pi*exp(1))) + 7/8);
\\ Подбор c (ширина окна)
c = 5;
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
while(E4 > target,
c += 0.5;
E4 = x*lam*sqrt(2)/log(x) * exp(-c^2/2)/(sqrt(Pi)*c^2);
);
while(E4 < target/3 && c > 3,
c_new = c - 0.5;
E4_new = x*lam*sqrt(2)/log(x) * exp(-c_new^2/2)/(sqrt(Pi)*c_new^2);
if(E4_new < target, c = c_new; E4 = E4_new, break);
);
\\ Остальные погрешности
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);
Etotal = E1 + E2 + E3 + E4;
\\ Стоимость простого в зависимости от режима
upper = x*exp(c*lam);
if(upper <= threshold,
cost_prime = cost_fast;
,
if(upper <= 2^63,
cost_prime = cost_mid;
,
cost_prime = cost_slow;
)
);
\\ Время
n_primes = 2*c*lam*x/log(x);
time_total = (cost_load + cost_zero)*nzeros + cost_prime*n_primes + 0.02;
\\ Обновление минимума
if(Etotal < 0.5 && 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("ВНИМАНИЕ: не найдены параметры с погрешностью < 0.5");
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;
upper = x*exp(c*lam);
if(upper <= threshold,
cost_prime = cost_fast;
,
if(upper <= 2^63,
cost_prime = cost_mid;
,
cost_prime = cost_slow;
)
);
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 = %d\n", nzeros);
printf(" c = %.1f\n", c);
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", cost_load*nzeros);
printf(" Ряд phihat: %.1f с\n", cost_zero*nzeros);
printf(" Решето: %.1f с\n", cost_prime*n_primes);
printf(" Сумма: %.1f с\n", time_total);
print("Вызов: ");
print("platt_pi(",x,",",nzeros,",",strprintf("%.4g",lam),",",500,",",strprintf("%.1f",c),")");
print();
\\ return([lam, nzeros, c, Etotal, time_total]);
} Кажется, что уже весь Платт теперь как на ладони. Так что, Yadryara , буду с нетерпением ждать применения Платта к кортежам 
|
