Вероятность появления взаимно простых чисел

[B3LYP]

Я сел покодить и написал программу, подсчитывающую вероятности перебором. Берутся все числа в интервале от 1 до 100, и двумя вложенными циклами подсчитывается вероятность, что они не являются взаимно-простыми. Потом берутся все числа в интервале от 1 до 200, потом все числа в интервале от 1 до 300 и так далее. Получился такой график (мой компьютер считал это час):



Красная кривая это теоретический предел $6/ \pi^2$. Я не то чтобы допускаю, что в какой-то момент чёрная кривая может уйти ниже красной, но мне хочется чтобы кто-то дал понятное доказательство, почему это невозможно. Хотелось бы доказательство “на пальцах” вроде такого:

Возьмём два ряда

$1+1/2+1/4+1/8+1/16… $

$1+1/2+1/3+1/4+1/5+1/6…$

Как посчитать эти суммы и доказать это? Понятное мне доказательство “на пальцах” такое: для первого ряда, с каждым слагаемым вдвое уменьшается расстояние до 2, соответственно с бесконечным числом слагаемых это расстояние станет бесконечно малым. Второй же ряд расходится, и для доказательства этого можно группировать члены, показав что будет бесконечное число слагаемых, каждое из которых не меньше 1/2.
Просьба дать мне ещё какие-нибудь такие доказательства или что-то ещё по теме.

[wrest]

B3LYP в сообщении #1725972 писал(а):

Как посчитать эти суммы и доказать это?

Первая сумма равна двум. Это очевидно например из такого рассуждения. Берём квадрат площадью два и разрезаем напополам. Будет 1+1=2. Одну из половин разрезаем напополоам, будет 1+1/2+1/2=2, затем разрезаем одну из половин ещё наполовину 1+1/2+1/4+1/4=2 и так до бесконечности. В сумме пока вы режете а часики тикают, на руках всегда будет площадь квадрата (два). Привет вам от Зенона, кстати
Можно и проще, если вы это уже проходили - первый ряд это геометрическая прогрессия со знаменателем 1/2 и первым членом единица.

Второй ряд, как известно, расходится (к логарифму количества членов плюс константа). Этому учат на первом курсе матана, ну можно убедиться и как вы предложили - группировкой.

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

B3LYP в сообщении #1725972 писал(а):

или что-то ещё по теме.

Фихтенгольц, Курс дифференциального и интегрального исчисления в 3 томах, том 2, глава 11.

[wrest]

B3LYP в сообщении #1725972 писал(а):

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

А чему у вас получилась равной вероятность для N=820?

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

B3LYP в сообщении #1725972 писал(а):

Получился такой график (мой компьютер считал это час):

Что-то долго И докуда досчитал, за час?

[Dmitriy40]

N=1276 даёт ещё чуть меньше вероятность.

До N=10000 у меня считалось 13.5с, до N=30000 считалось 126с. За час наверное досчиталось бы до N=160000 (время растёт квадратично). Ну правда я особо и не заморачивался.

[tolstopuz]

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

Да, там самое неинтуитивное, имхо, это то, что
$\sum \frac{1}{n^2}=\prod \dfrac{1}{1-p^{-2}}$
Где слева сумма по всем натуральным $n$, а справа - произведение по всем простым $p$

Это-то как раз довольно интуитивно. Расписываю для B3LYP:

Вероятность того, что оба числа делятся на $2$, равна $1/2^2$, что на $3$ - $1/3^2$, что на $5$ - $1/5^2$ и так далее. Наивно предположим, что все эти вероятности независимы и их можно перемножать. На самом деле предположение довольно разумное, потому что, скажем, чисел, делящихся на $15$, как раз $1/15$ от всех, то есть $1/3\cdot1/5$. Если взять конечный отрезок, на конце будет небольшое расхождение, но при увеличении размера отрезка оно все меньше влияет.

Получим итоговую вероятность, что оба числа совместно не делятся ни на $2$, ни на $3$, ни на $5$, ни вообще ни на какое простое:

$$P=\left(1-\frac1{2^2}\right)\cdot\left(1-\frac1{3^2}\right)\cdot\left(1-\frac1{5^2}\right)\cdot\dotsc.$$

Нам легче работать с $1/P$:

$$\begin{align*}1/P &=\frac1{1-\frac1{2^2}}\cdot\frac1{1-\frac1{3^2}}\cdot\frac1{1-\frac1{5^2}}\cdot\dotsc\\ &=\left(1+\frac1{2^2}+\frac1{(2\cdot2)^2}+\dotsc\right)\cdot\left(1+\frac1{3^2}+\frac1{(3\cdot3)^2}+\dotsc\right)\cdot\left(1+\frac1{5^2}+\frac1{(5\cdot5)^2}+\dotsc\right)\cdot\dotsc\\ &=1+\frac1{2^2}+\frac1{3^2}+\frac1{(2\cdot2)^2}+\frac1{5^2}+\frac1{(2\cdot3)^2}+\dotsc\\ &=1+\frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\frac1{5^2}+\frac1{6^2}+\dotsc, \end{align*}$$

где мы сначала воспользовались формулой суммы бесконечной геометрической прогрессии, а потом тем, что каждое натуральное число единственным способом раскладывается на простые множители, поэтому встретится в сумме ровно один раз.

Вот почему эта сумма равна $\frac{\pi^2}6$ - вопрос, на мой взгляд, намного менее интуитивный. Хотя Эйлер вроде тоже ничего запредельного не делал, просто теорему Вейерштрасса придется принять на веру:
https://ru.wikipedia.org/wiki/%D0%A0%D1 ... 0%B4%D0%B0

[B3LYP]

Мне захотелось поразвлечься и ускорить свой алгоритм: новый алгоритм посчитал предыдущий график в 200 раз быстрее, а за час посчиталось до N=325:



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

[Dmitriy40]

B3LYP
Для каждого следующего N не надо пересчитывать всё заново, достаточно добавить только комбинации именно с этим N в качестве одного из чисел. А учитывая симметрию между числами, то можно сразу за N взять скажем первое число и перебирать только второе, увеличивая счётчик вероятностей сразу на 2.
Так сложность падает с кубической (цикл по N, два цикла по двум числам) до квадратичной (цикл по N, цикл по одному числу).
И смотря какой язык используете, можно взять быструю реализацию gcd(a,b).
И за час реально досчитать до сотни тысяч N.

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

Рекомендую проверить N=820, 1276, 1422, 1926, 2080, 2640, 3160, 3186, 3250, 4446, 4720, 4930.

[mihaild]

А еще НОД тут считать вообще не нужно. Ответ считается более-менее за один проход решетом.

[B3LYP]

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

Для каждого следующего N не надо пересчитывать всё заново, достаточно добавить только комбинации именно с этим N в качестве одного из чисел. А учитывая симметрию между числами, то можно сразу за N взять скажем первое число и перебирать только второе, увеличивая счётчик вероятностей сразу на 2.
Так сложность падает с кубической (цикл по N, два цикла по двум числам) до квадратичной (цикл по N, цикл по одному числу).

Это не очень существенно для того, что мы ищем. То что вы написали - заурядные подходы к оптимизации, а мы ищем математико-алгоритмические подходы. И снижение сложности с кубической до квадратичной - не очень впечатляет (мы же хотим научиться быстро разлагать любое число на простые множители без квантового компьютера, т.е. понизить сложность в гуглы раз). Я позже ещё ускорю свой алгоритм за счёт всяких решёт Эратосфена, но пока предлагаю обсудить, кому интересно, какие подходы я использовал на данный момент.

[Dmitriy40]

B3LYP в сообщении #1726170 писал(а):

(мы же хотим научиться быстро разлагать любое число на простые множители без квантового компьютера, т.е. понизить сложность в гуглы раз)

Вообще-то тема заявлена про конкретный вопрос о вероятности про два случайных числа. При чём здесь разложение на множители?!

[B3LYP]

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

Вообще-то тема заявлена про конкретный вопрос о вероятности про два случайных числа. При чём здесь разложение на множители?!

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

[wrest]

B3LYP в сообщении #1726163 писал(а):

Заодно, если кому интересно, попробуйте догадаться – как работал раньше мой алгоритм и как я его ускорил?

Да кто ж вас знает, у меня до 325 планшет на андроиде считает за 15мс Так что можете ещё на 5-7 порядков ускориться.
pari/gp:
pt(n)=my(s=sum(i=1,n,eulerphi(i)));return((2*s-1)/(n^2))
И это каждое отдельно если считать.
Сразу все 325 считает за ноль миллисекунд, более-менее измеримое время -- все вероятности до n=5000 считает за 15мс
pt_vector(n)=my(v=vector(n),s=0);for(k=1,n,s+=eulerphi(k);v[k]=(2*s-1)/(k^2));return(v);

[Dmitriy40]

wrest
Круто!

А если не делать функцию от N, а сразу считать следующее N из предыдущего, то по идее ещё быстрее:
? m=6/Pi^2; q=0; s=-1; for(n=1,1e8, s+=eulerphi(n)*2; if(s<m*n^2, q++); ); q
time = 1min, 53,069 ms.
%1 = 182255
Количество исключений до N=10^8.

[B3LYP]

Ну, если никому тут не интересно обсуждать общие приёмы алгоритмической оптимизации, жаль...
Сейчас я ещё подускорю свой код с помощью решета Эратосфена, но надо пояснить, что далеко не всякая оптимизация мне вообще нужна.
Сейчас нет смысла заниматься какой-то мало-мальски низкоуровневой оптимизацией, потому что мне хочется понять что-то в математике, чтобы не было “за деревьями не видно леса”. Пример: если мы считаем P(1000), то есть вероятность того что два случайных числа от 1 до 1000 не имеют общего делителя, то это у меня делается перебором в двух вложенных циклах, и ясно что можно ускорить этот перебор в два раза, начиная второй счётчик не с 1 а с первого счётчика (использовать знание, что факт того что A и B не имеют общего делителя совпадает с аналогичным фактом что B и A не имеют общего делителя). Но такая оптимизация дополнительно усложняет код, соответственно увеличивает вероятность ошибок в чём-то более важном.
Мой код ужасно медленный просто потому, что там многократно дублируются вычисления. Если код считает сначала P(100), потом P(200), потом P(300) и так далее, то на тысячной итерации, т.е. когда доходит до P(100000), все предыдущие итерации оказываются ненужными, и суммарно график строится почти в тысячу раз медленнее, чем если сразу посчитать P(100000). Можно это ускорить за счёт того, что когда мы посчитали P(99900) и далее считаем P(100000), можно перебирать только пары чисел от 99901 до 100000, и далее правильно всё просуммировать. Но зачем мне это делать, если это увеличивает сложность кода и соответственно вероятность ошибок?
Поэтому я сейчас свой код даже замедлил: написал функцию CalcPNonCoprime(N:integer), которая каждый раз инициирует и заполняет все массивы. Пусть она запускается 10 раз; за 22 минуты мой код с предыдущим алгоритмом дал такой график:



Сейчас я его ускорю решетом Эратосфена (или может это как-то ещё называется).

[Null]

Можно посчитать $F(N)=\left [\frac{N}{1^2}\right ]-\left [\frac{N}{2^2}\right ]-\left [\frac{N}{3^2}\right ]-\left [\frac{N}{5^2}\right ]+\left [\frac{N}{6^2}\right ]\dots$ - формула включений-исключений. Сложность $\sqrt N$. Для вычисления функции Мебиуса используем Решето Эратосфена и еще массив (в начале все 1), числа кратные очередному $p$ умножаем на $-1$, а кратные $p^2$ зануляем.