кортежи последовательных простых. ключ к 19-252

[Dmitriy40]

Вынос BC из функции довольно прост, чего его показывать. Как и кучу вложенных циклов, тривиально же. Как и next если CalcC() вернула 0. Остальные улучшения не дали большого эффекта. Вот сделаю битовые маски, как с ww[], тогда и покажу, это столь кардинально что пока всё ещё обдумываю.
Треугольник Паскаля выдаётся командой binomial(108)[1..31] (с чистого до максимально грязных). Но реально паттернов сильно меньше, я показываю точное количество.
Но реально в коде он никак не используется, скобка с логарифмами в степени это та самая формула включений-исключений как указал vicvolf, биномиальных коэффициентов в ней не нужно. А погрешность возникает из-за недоучёта загрязнений, видимо если учесть все, то всё хорошо сойдётся и без подгонок как у меня.

[vicvolf]

Мне кажется, что назрел вопрос - какие будут коэффициенты $C(m_1.m_2,...m_k)$ в первой гипотезе Харди-Литтлвуда, в случае, когда кортеж не проходит по какому-то модулю $p$, т.е. решения его сравнения $x(x+m_1)(x+m_2)...(x+m_k)=0(modp)$ образуют полную систему вычетов по модулю $p$.

В этом случае, количество решений сравнения $w(p)=p$ и значение $C(m_1.m_2,...m_k)=\prod_{p \geq 3} \frac{1-w(p)/p}{1-1/p)^k}=0$.

Таким образом, если кортеж не проходит по какому-либо модулю $p$, то количество таких кортежей на любом диапазоне будет нулевым.

[Dmitriy40]

vicvolf
Этот вопрос уже решён ещё когда я показал код программы, там комментарий именно об этом. Не из формулы, а эвристически: если паттерн имеет полную систему вычетов, то такой кортеж может быть только однажды и значит все дальнейшие вероятности равны 0. Что оно на самом 1 (штука) в каких-то ограниченных сверху диапазонах роли не играет, это возможно только в начале числового ряда (по простым более длины паттерна разрешённые остатки есть всегда), где нужных нам решений точно нет (проверяется за секунды, а мы уже годы мучаемся).

-- 25.07.2024, 11:16 --

Но вот то что они могут равны 0 - очень даже важно! Например для C8 можно считать не 352025629371 паттернов, а лишь 165735854, в 2124 раза меньше, т.е. очень намного быстрее.

-- 25.07.2024, 11:20 --

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

А погрешность возникает из-за недоучёта загрязнений, видимо если учесть все, то всё хорошо сойдётся и без подгонок как у меня.

Неа, не сходится, нижнюю границу интегрирования таки надо выбирать разумно, а не просто ставить с запасом от 0 или 2.

[Yadryara]

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

Треугольник Паскаля выдаётся командой binomial(108)[1..31] (с чистого до максимально грязных).

Нет, этой командой выдаётся 109-я строка треугольника.

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

Но реально паттернов сильно меньше, я показываю точное количество.

В курсе, что $6$ меньше чем $440988985822461767289318384$.

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

это та самая формула включений-исключений как указал vicvolf, биномиальных коэффициентов в ней не нужно.

Про формулу включений-исключений тоже в курсе. Ну а номинальное количество проходов вложенных циклов чему равно?

Если Вы покажете код, мне проще будет показать как его ускорить. Чем на пальцах объяснять.

Хотя, попробую.

Вот у нас есть симметричный нечётный паттерн с трёхзначным диаметром $d$.

Допустим, его можно загрязнить 2-кой. А можно ли его загрязнить числом $d-2$ ? Да. При этом количество кортежей для обоих вариантов будет одинаковым.

Допустим, его можно загрязнить 2-кой и 8-кой. А можно ли его загрязнить числами $d-2$ и $d-8$ ? Да. При этом количество кортежей для обоих вариантов будет одинаковым.

Дальше объяснять?

[Dmitriy40]

Собственно код уже показывал, его доработки минимальны:

dd=setminus(vector(v[#v]/2,i,i*2),v); C1=C2=C3=0; C = CalcC(v); print("C =",C); nn=vector(10);
for(a=1,#dd, t = CalcC(concat(dd[a],v)); if(t==0, next); C1 += t; nn[1]++;
   for(b=1,a-1, t=CalcC(concat([dd[a],dd[b]],v)); if(t==0, next); C2 += t; nn[2]++;
      for(d=1,b-1, C3 += CalcC(concat([dd[a],dd[b],dd[d]],v)); nn[3]++; )));
print("C1=",C1); print("C2=",C2); print("C3=",C3); print(concat(1,nn));

Сделать в первом цикле #dd/2 и 2*t конечно можно, C1 будет правильной, а C2 и C3 нет.

-- 25.07.2024, 14:10 --

За 9ч посчиталось C9, получилось 440487592 паттернов:
C9=15994906029871346067126.943725676318187
Доля чистых для 1e25 уменьшилась с 7.1% до 6.9%. Полагаю пока не посчитана C10 можно долю чистых принять за 7%.
Количество чистых кортежей 0.6 до 1e25 и 3.2 до 1e26. Беда. Столько я не просчитаю, это лет 20.

-- 25.07.2024, 14:35 --

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

Сделать в первом цикле #dd/2 и 2*t конечно можно, C1 будет правильной, а C2 и C3 нет.

Изменив второй цикл на for(b=1,#dd, if(b==a, next); t=... можно сделать правильной и C2, но это не ускоряет программу, а замедляет в 1.5 раза.

[Yadryara]

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

Сделать в первом цикле #dd/2 и 2*t конечно можно, C1 будет правильной, а C2 и C3 нет.

Здесь не так просто. Подумаю. Или уже не надо?

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

Доля чистых для 1e25 уменьшилась с 7.1% до 6.9%. Полагаю пока не посчитана C10 можно долю чистых принять за 7%.

Скорее, так и останется на 6.9%. Разве что пару сотых добавит.

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

Количество чистых кортежей 0.6 до 1e25 и 3.2 до 1e26.

Да, это прям по нижним границам моего давнего прогноза с 44-й страницы:

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

В плохом варианте:

$$\frac{52\cdot0.07}{6}\approx 0.607$$

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

В плохом варианте:

$$\frac{246.7\cdot0.08}{6.2}\approx 3.18$$

Иными словами, 1 чистый кортеж на 1.8-1.9 e25. Вы это сейчас можете поточнее посмотреть.

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

Беда. Столько я не просчитаю, это лет 20.

Остаётся коллективный счёт. Хотя здесь это, видимо, оффтоп.

[Dmitriy40]

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

Остаётся коллективный счёт.

Чтобы время снизилось до разумного надо под тысячу вычислительных потоков (и без учёта гипертрейдинга, он слишком малый выигрыш даёт).
Максимум что я могу ещё сделать - ускорить программу раз в 5 и тогда можно обойтись сотней-двумя потоков, но я это с конца января никак не доделаю ...

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

Подумаю. Или уже не надо?

Смотрите сами. Зависит от выигрыша скорости, если вдвое то смысла мало, если на порядок - хорошо бы.
Да и C10 могу прежним способом посчитать (чуть более суток). Если же считать более высокие, то буду делать на битовых масках, а потом уже их на асм/дельфи если понадобится.
В принципе меня и такая точность уже устраивает, можно считать задачу оценки количества кортежей 19-252 выполненной.

-- 25.07.2024, 18:08 --

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

В этом случае, количество решений сравнения $w(p)=p$ и значение $C(m_1.m_2,...m_k)=\prod_{p \geq 3} \frac{1-w(p)/p}{(1-1/p)^k}$.

О! До меня дошло почему в программе вычисления непохожи на формулу 4.2 из вашей статьи, процитированную вами и здесь: в программе это всего лишь правильное вычисление $w(p)$.
Что же, тогда могу посчитать дисперсию. Возьмём посчитанную C=1592669394.7454967047727717796940585531 для всех грязных кортежей и подставим в формулу 4.10 (она же 4.9):

? C=1592669394.7454967047727717796940585531;
? for(po=24,30, print("1e",po,": ",sqrt(C*intnum(t=10^15, 10^po, 1/log(t)^19)-C^2*intnum(t=10^15, 10^po, 1/log(t)^38))))
1e24: 1.3898737007443956264543584586283761667
1e25: 2.9485998036498880172894837364896035915
1e26: 6.3586791580540810367947900699815980154
1e27: 13.920483119568267774553562711600895874
1e28: 30.901005695629591654537037552056588615
1e29: 69.482357971925702784379376511905128366
1e30: 158.11099981301525670893684416378521668

Напомню, это дисперсия вот этих величин:

? for(po=24,30, print("1e",po,": ",C*intnum(t=10^15, 10^po, 1/log(t)^19)))
1e24: 1.9317489040209218077620052012016239896
1e25: 8.6942408020841581689260306347508468737
1e26: 40.432800635071356887023073421493660922
1e27: 193.77985028218509208681211355385055627
1e28: 954.87215300133246363014344258375574022
1e29: 4827.7980693388272613198557464975111254
1e30: 24999.088261871310541976721713820598466

Наверное могу посчитать и дисперсию чистых (с точностью по C9), проверьте правильно ли возвёл в квадрат в L2:

? C =1592669394.7454967047727717796940585531;
? C1=224048612037.93039512469294753628944180;
? C2=15060404457776.501102314463450218573948;
? C3=643975079051969.03166225028027252021444;
? C4=19672864014876028.898119063672064125711;
? C5=457178627398399657.21129392551944562299;
? C6=8402763498659303904.0898362513026471087;
? C7=125377224477963818544.79603392112737254;
? C8=1546929025718799128942.3523933275147810;
? C9=15994906029871346067126.943725676318187;
? {for(po=24,30,
   L=intnum(t=10^15, 10^po, 1/log(t)^19 * (C - C1/log(t) + C2/log(t)^2 - C3/log(t)^3 + C4/log(t)^4 - C5/log(t)^5 + C6/log(t)^6 - C7/log(t)^7 + C8/log(t)^8 - C9/log(t)^9) );
   L2=intnum(t=10^15, 10^po, C^2/log(t)^38 - C1^2/log(t)^40 + C2^2/log(t)^42 - C3^2/log(t)^44 + C4^2/log(t)^46 - C5^2/log(t)^48 + C6^2/log(t)^50 - C7^2/log(t)^52 + C8^2/log(t)^54 - C9^2/log(t)^58);
   print("1e",po,": ",L," +- ",sqrt(L-L2));
);}
1e24: 0.11829355770469481429647974273126314546 +- 0.34393830508493062337047071365040783375
1e25: 0.60411176995236903553222204074424798256 +- 0.77724627368188072615387211369957237360
1e26: 3.1498526059163907923139267256467417295 +- 1.7747824108651716047730631021301682250
1e27: 16.756735512063577967072300620494151094 +- 4.0934991769955906474672652227206090829
1e28: 90.872229370114848187617734915845724712 +- 9.5326926610541078635019494768486222346
1e29: 501.88638627662927717262639734007969388 +- 22.402820944618320034461196549381089396
1e30: 2820.3632400347586454141832210819284909 +- 53.107092182068852614497480434282823882

[vicvolf]

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

Что же, тогда могу посчитать дисперсию. Возьмём посчитанную C=1592669394.7454967047727717796940585531 для всех грязных кортежей и подставим в формулу 4.10

Это точнее не дисперсия, а среднее квадратичное отклонение, т.е. корень квадратный из дисперсии. Но нам оно и нужно.
Например, для интервала до $10 ^{25}$ количество кортежей 19-252 находится с вероятностью $0,99$ в диапазоне:

$8.6942408020841581689260306347508468737-3 \cdot 2.9485998036498880172894837364896035915,$$8.6942408020841581689260306347508468737+3 \cdot 2.9485998036498880172894837364896035915$

Таким образом, с вероятностью 0,99 в диапазоне до $10^{25}$ количество кортежей 19-252 колеблется, грубо говоря, от 0 до 19.

[Dmitriy40]

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

количество кортежей 19-252

Только количество произвольно загрязнённых кортежей. Это очень существенно. Без уточнения "кортежи 19-252" подразумеваем чистые кортежи, которых сильно меньше. Тем более что загрязнённые кортежи формально будут уже не 19-, а больше. Ещё обозначение 19-252 используем как указание на вариант задачи, что речь про эти кортежи.

[Dmitriy40]

Ради интереса посмотрел на оценку кортежей 17-240, их известно 5шт до 1e22, причём первый почти точно на 1e21. Оценка до шестикратного загрязнения:

v=[0, 6, 24, 36, 66, 84, 90, 114, 120, 126, 150, 156, 174, 204, 216, 234, 240]
C =204267977.27052456200777283266142295380
C1=28369943870.637190788534643770043972079
C2=1882883655575.2751806520561345784915500
C3=79503105963393.861809247903002671637974
C4=2398827502988278.8370655540475223328881
C5=55074965123455739.447710661923258904122
C6=1000416540944375388.0187642856124195935
[1, 104, 2952, 41852, 379420, 2471696, 12318076]
10^20: 0.364217
10^21: 1.209901
10^22: 4.437549
10^23: 17.921522
10^24: 78.976390
10^25: 375.024404

Считаю совпадение отличным.
И кстати это довод что меньшего кортежа нет.

[Dmitriy40]

Переделал с векторов на битовые маски, но особо не оптимизировал, ускорение получилось (старое время / новое = ускорение):
C3: 3.60с/1.17с=3.1х
C4: 33.6с/10.4с=3.2х
C5: 5м/1.2м=4.2х
C6: 25м/6.3м=4х
Что-то до порядка не дотянуло.

Код:

v=[0,6,12,30,42,72,90,96,120,126,132,156,162,180,210,222,240,246,252]; print("v=",v);
BC=vector(#v+10,k, prodeulerrat(( p^k - k*p^(k-1) )/(p-1)^k, 1, nextprime(k+1)) );
MC=vector(#v+10,k, x=1.0;forprime(p=3,k,x/=p*(1-1.0/p)^k);forprime(p=k+1,v[#v]/2,x/=p-k);x );
CC=vector(#v+10,k, 2^(k-1) * MC[k] * BC[k]);
dd=setminus(vector(v[#v]/2,i,i*2),v); C1=C2=C3=C4=C5=C6=C7=C8=C9=C10=0; nn=vector(10);
v0=vector(v[#v]/2,p, if(p>2&&isprime(p),setminus(vector(p,i,i-1),Set(-v%p)),[]));
m0=vector(#v0,p, t=0;foreach(v0[p],x, t=bitor(t,2^x););t);
C=CC[#v]; forprime(p=3,#m0, C*=hammingweight(m0[p]); ); printf("C =%f\n",C);
for(a=1,#dd,
   forprime(p=3,#v+1, if(m0[p]==2^(-dd[a]%p), next(2)); );
   m1=m0; t=CC[#v+1]; forprime(p=3,#m1, m1[p]=bitnegimply(m1[p],2^(-dd[a]%p)); t*=hammingweight(m1[p]); ); C1 += t; nn[1]++;
   for(b=1,a-1,
      forprime(p=3,#v+2, if(m1[p]==2^(-dd[b]%p), next(2)); );
      m2=m1; t=CC[#v+2]; forprime(p=3,#m2, m2[p]=bitnegimply(m2[p],2^(-dd[b]%p)); t*=hammingweight(m2[p]); ); C2 += t; nn[2]++;
      for(d=1,b-1,
         forprime(p=3,#v+3, if(m2[p]==2^(-dd[d]%p), next(2)); );
         m3=m2; t=CC[#v+3]; forprime(p=3,#m3, m3[p]=bitnegimply(m3[p],2^(-dd[d]%p)); t*=hammingweight(m3[p]); ); C3 += t; nn[3]++;
         for(e=1,d-1,
            forprime(p=3,#v+4, if(m3[p]==2^(-dd[e]%p), next(2)); );
            m4=m3; t=CC[#v+4]; forprime(p=3,#m4, m4[p]=bitnegimply(m4[p],2^(-dd[e]%p)); t*=hammingweight(m4[p]); ); C4 += t; nn[4]++;
););););
print("C1=",C1); print("C2=",C2); print("C3=",C3); print("C4=",C4); print(concat(1,nn));

Думаю логика добавления вложенных циклов понятна, приводить все не буду.

[vicvolf]

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

Ради интереса посмотрел на оценку кортежей 17-240, их известно 5шт до 1e22, причём первый почти точно на 1e21. Оценка до шестикратного загрязнения:

v=[0, 6, 24, 36, 66, 84, 90, 114, 120, 126, 150, 156, 174, 204, 216, 234, 240]
C =204267977.27052456200777283266142295380
C1=28369943870.637190788534643770043972079
C2=1882883655575.2751806520561345784915500
C3=79503105963393.861809247903002671637974
C4=2398827502988278.8370655540475223328881
C5=55074965123455739.447710661923258904122
C6=1000416540944375388.0187642856124195935
[1, 104, 2952, 41852, 379420, 2471696, 12318076]
10^20: 0.364217
10^21: 1.209901
10^22: 4.437549
10^23: 17.921522
10^24: 78.976390
10^25: 375.024404

Считаю совпадение отличным.
И кстати это довод что меньшего кортежа нет.

Ради интереса подсчитал средне квадратичные отклонения для данного кортежа [0, 6, 24, 36, 66, 84, 90, 114, 120, 126, 150, 156, 174, 204, 216, 234, 240] на том же диапазоне.
C =204267977.27052456200777283266142295380;
for(po=20,25, print("1e",po,": ",sqrt(C*intnum(t=10^15, 10^po, 1/log(t)^17)-C^2*intnum(t=10^15, 10^po, 1/log(t)^34))))

1e20: 1.3250375467236887150272578898712175202
1e21: 2.7267461085413125554224436286046511819
1e22: 5.7290417485650859581276282740212296666
1e23: 12.268690918529430383401940989756943337
1e24: 26.733890778784400798366380695907488922
1e25: 59.186018919027202382429157581189982916

Что-то многовато получилось. За основу взял Вашу программу.

[Yadryara]

Dmitriy40, у меня пока так получилось. Не шибко уклюже.

a=setminus(vector(v[#v]/2,i,i*2),v);

for(i1=1,#a/2,
   forprime(p=3,#v+1, if(m0[p]==2^(-a[i1]%p), next(2)); );
   m1=m0; t=CC[#v+1];
 forprime(p=3,#m1, m1[p]=bitnegimply(m1[p],2^(-a[i1]%p));
 t*=hammingweight(m1[p]); ); C1 +=2*t; nn[1]++;

   for(i2=i1+1,#a,
      forprime(p=3,#v+2, if(m1[p]==2^(-a[i2]%p), next(2)); );
      m2=m1; t=CC[#v+2];
 forprime(p=3,#m2, m2[p]=bitnegimply(m2[p],2^(-a[i2]%p));
 t*=hammingweight(m2[p]); ); C2 +=(2-floor((i2*2-1)/#a))*t; nn[2]++;

      for(i3=i2+1,#a,
         forprime(p=3,#v+3, if(m2[p]==2^(-a[i3]%p), next(2)); );
         m3=m2; t=CC[#v+3];
 forprime(p=3,#m3, m3[p]=bitnegimply(m3[p],2^(-a[i3]%p));
 t*=hammingweight(m3[p]); ); C3 +=(2-floor((i3*2-1)/#a))*t; nn[3]++;

         for(i4=i3+1,#a,
             forprime(p=3,#v+4, if(m3[p]==2^(-a[i4]%p), next(2)); );
            m4=m3; t=CC[#v+4];
 forprime(p=3,#m4, m4[p]=bitnegimply(m4[p],2^(-a[i4]%p));
 t*=hammingweight(m4[p]); ); C4 +=(2-floor((i4*2-1)/#a))*t; nn[4]++;

))));

[Dmitriy40]

Выражение

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

(2-floor((i2*2-1)/#a))*t

может принимать лишь два значения, 2t и t, причём граница между ними ровно #a/2, так что можно заменить на if(i2>#a/2, t, t*2). На скорость правда не влияет, но наверное более понятно.
Насколько понимаю других коррекций нет, только пределы циклов и множитель для t.
Вычисления дали чуть большую погрешность, 2-3 младшие цифры (из 38), плевать.
Скорость фактически не изменилась: 9.2с (была) vs 8.8с (стала) для C4 и 64с vs 63с для C5. Что видно и по количеству просчитанных паттернов для C5: [1, 108, 3124, 44810, 408340, 2658128] vs [1, 54, 2355, 39594, 386756, 2597132]. Т.е. не вижу никакого уменьшения количества паттернов вдвое, ну кроме как для C1.

-- 26.07.2024, 19:00 --

Последней программой посчитал C10 за 8.25ч, получилось 985580430 паттернов:
C10=139968081763087868053405.50716615643237
Доля чистых для 1e25 стала 6.98%. Полагаю 7% достаточно точно для практики, вот какие были доли чистых для C7-C10 (меньшие ненадёжны): 6.105%, 7.142%, 6.946%, 6.979%, следующее не должно уйти ниже 6.946%.

Версия с векторами ещё не завершилась, прошло 21ч, надо ещё несколько (может и больше десятка).

[Yadryara]

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

6.105%, 7.142%, 6.946%, 6.979%, следующее не должно уйти ниже 6.946%.

Ну вот, говорил новый обсчёт пару сотых добавит, добавил 3 сотых. Очень быстро сходится, так что не уйдёт ниже 6.972%.

Вашу последнюю модификацию проверял на 7-ках и 9-ке с 8-9-кратным загрязнением — прекрасно сошлись. Чистые проверять гораздо легче, их в Базах полно.