Эвристический метод нахождения некоторых простых близнецов

[Dmitriy40]

wrest
ИИ такой ИИ ...

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

/* Приводим к стандартному виду a*x^2 + b*x + c */
[всё пропущенное - лишнее]
a = polcoef(A, 2); b = polcoef(A, 1); c = polcoef(A, 0)

Ну и ИИ не знал что в PARI допустима конструкция if-elseif, так что вложенный if со своей закрывающей скобкой не нужен, можно сразу писать следующее (произвольное) условие в ветке else и его ветку true, а в ветку else можно снова произвольное условие и так далее до посинения.

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

/* Находим НОД коэффициентов */

Зачем? Если далее не используется ...

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

N = 500

Стоило бы наверное проверять не до фиксированной границы, а до НОК коэффициентов или хоть до максимального из них.

[wrest]

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

Зачем? Если далее не используется ...

Это адаптациия другого текста где используется

Впрочем, адаптация не очень-то, согласен.
Запощу попозже получшее.

[wrest]

Ну вот, вроде готово.
Итак, даётся пара функций
PrimeRateSinge(A) и PrimeRateTwin(A)
в качестве A принимается квадратный трёхчлен.
Функция PrimeRateSingle считает поправочную константу для простых которую все мы любим, но названия у неё нет.
Функция PrimeRateTwin возвращает эту константу для простых близнецов.
Связь с плотностью тут такая. Примем плотность простых близнецов, которые возвращает функция$ f(n)=n$ за единицу. Имеется в виду что $f(n)$ и $f(n)+2$ оба простые. Тогда, подав на вход функции PrimeRateTwin интересующий нас полином, например $g(n)=33n^2+2475n-1$ мы получим результат около 41. Делим его на учетверённую "константу близнецов 0.66" равную $4\cdot 0.66=2.64$ и получаем примерно 15
Таким образом, $g(n)=33n^2+2475n-1 $возвращает левого простого близнеца в 15 раз чаще чем $f(n)=n$ для чисел того же размера.
Константа близнецов для справки $0.66016181584686957392..$.

Функции

(Оффтоп)

PrimeRateSingle(A, N = 50) = {
  my(D = poldisc(A), S, P, lim, S1, S2, L, v, Ebad, B, n, d, fan, X, Y, a, b, c, D0, f, p, s, t, n_temp);
  
  if(poldegree(A) != 2, error("не квадратный многочлен!"));
  [a, b, c] = Vec(A);
  if(issquare(D) || gcd([2*a, a+b, c]) > 1, return(0));
  N = max(N, 3);
  
  S = if((a + b) % 2, 1., 2.);
  
  n_temp = a >> valuation(a, 2);
  P = factor(n_temp)[,1];
  foreach(P, p, S *= if(b % p, (p - 1) / (p - 2), p / (p - 1)));
  
  [D0, f] = coredisc(D, 1);
  n_temp = f >> valuation(f, 2);
  P = factor(n_temp)[,1];
  S /= vecprod([1 - kronecker(D0, p) / (p - 1) | p <- P]);
  
  S *= prodeuler(p = 3, N, 1 - kronecker(D0, p) / (p - 1));
  
  B = getlocalbitprec();
  lim = ceil(B * log(2) / log(N / 2));
  localbitprec(B + lim + exponent(lim));
  
  L = lfuninit(D0, [1/2, lim, 0]);
  v = vector(lim);
  
  forfactored(X = 1, lim,
    [n, fan] = X;
    S1 = 0;
    P = fan[,1];
    if(n % 2 == 0, P = P[^1]);
    X = matconcat([P, vectorv(#P, i, 1)]);
    fordivfactored(X, Y,
      [d] = Y;
      S1 += moebius(Y) << (n/d)
    );
    v[n] = S1 / (2*n)
  );
  
  Ebad = [];
  P = setunion(factor(abs(D0))[,1]~, primes([2, N]));
  for(t = 1, #P,
    p = P[t];
    s = kronecker(D0, p);
    if(s, 
      Ebad = concat(Ebad, [1 - s * 'x]);
      Ebad = concat(Ebad, p)
    )
  );
  
  my(Pbad = [], Epoly = []);
  forstep(t = 1, #Ebad, 2,
    Pbad = concat(Pbad, Ebad[t+1]);
    Epoly = concat(Epoly, Ebad[t])
  );
  Ebad = [Pbad, Epoly];
  
  S1 = sum(n = 1, lim, 
    v[n] * log(
      lfun(L, n) * prod(j = 1, #Pbad, 
        subst(Epoly[j], 'x, Pbad[j]^(-n))
      )
    )
  );
  
  S2 = sum(n = 2, lim, 
    (v[n] - if(n % 2 == 0, v[n/2])) * 
    log(zeta(n) * prod(j = 1, #P, 1 - P[j]^(-n)))
  );
  
  return(S * exp(-(S1 + S2)));
};


PrimeRateTwin(A, N = 500) = {
  my(A2, C1, C2, S_corr, a, b, c, D, Dp, D0, D0p, f, fp, bad, p,
     wp, wp1, wp2, x, num, den);
  
  if(poldegree(A) != 2, error("не квадратный многочлен"));
  [a, b, c] = Vec(A);
  
  D = b^2 - 4*a*c;
  Dp = b^2 - 4*a*(c+2);
  
  if(issquare(D) || issquare(Dp) || gcd([2*a, a+b, c, c+2]) > 1, return(0));
  
  [D0, f] = coredisc(D, 1);
  [D0p, fp] = coredisc(Dp, 1);
  
  A2 = A + 2;
  
  C1 = PrimeRateSingle(A, N);
  C2 = PrimeRateSingle(A2, N);
  
  bad = Set(concat(factor(abs(a))[,1]~, 
         concat(factor(abs(D0))[,1]~, factor(abs(D0p))[,1]~)));
  
  S_corr = 1.0;
  foreach(bad, p,
    if(p != 2,
      wp = 0; wp1 = 0; wp2 = 0;
      for(x = 0, p-1,
        if((a*x^2 + b*x + c) % p == 0, wp1++);
        if((a*x^2 + b*x + c + 2) % p == 0, wp2++);
        if((a*x^2 + b*x + c) % p == 0 || (a*x^2 + b*x + c + 2) % p == 0, wp++);
      );
      num = (1 - wp/p) / (1 - 1/p)^2;
      den = ((1 - wp1/p) / (1 - 1/p)) * ((1 - wp2/p) / (1 - 1/p));
      if(den != 0, S_corr *= num / den);
    );
  );
  
  return(C1 * C2 * S_corr);
};

Функция PrimeRateSingle это обёрнутый в одну функцию набор скриптов HardyLittlewood2.gp из книги
Karim Belabas and Henri Cohen, Numerical Algorithms for Number Theory
скрипты доступны тут: https://www.math.u-bordeaux.fr/~kbelaba ... lgorithms/
Ну а функция PrimeRateTwin запускает функцию PrimeRateSingle два раза, перемножает результаты и вносит всякие поправки.

Запуск:

? PrimeRateTwin(33*x^2+7425*x-1)
time = 15,764 ms.
41.092195190265525928551802759198033134
? PrimeRateTwin(x^2+1)
time = 59 ms.
1.5385570122938306044745759699613277308