Felix' Math Place (Posts about coprime.)
https://math.fontein.de/tag/coprime.atom
2024-01-05T07:05:27Z
Felix Fontein
Nikola
The Probability That Two Numbers Are Coprime.
https://math.fontein.de/2012/07/10/the-probability-that-two-numbers-are-coprime/
2012-07-10T23:30:36+02:00
2012-07-10T23:30:36+02:00
Felix Fontein
<div>
<div>
<p>
Let <span class="inline-formula"><img class="img-inline-formula img-formula" width="11" height="8" src="https://math.fontein.de/formulae/CiIJDoNXXhwwshmAknaOy.cbqWs.Z_qmDZe21A.svgz" alt="n" title="n"></span> be a natural number. Denote by <span class="inline-formula"><img class="img-inline-formula img-formula" width="20" height="15" src="https://math.fontein.de/formulae/x5nPIESCq4Gi.JgNYJUNJR.NY6G3n00fWAbwNw.svgz" alt="S_n" title="S_n"></span> the set
</p>
<div class="display-formula">
<img class="img-display-formula img-formula" width="315" height="18" src="https://math.fontein.de/formulae/Cw_Yx7GtxzVMgEONwddOXM0deXwKFUhs.yygfg.svgz" alt="\{ (x, y) \in \N \mid 1 \le x, y \le n, \; \gcd(x, y) = 1 \}" title="\{ (x, y) \in \N \mid 1 \le x, y \le n, \; \gcd(x, y) = 1 \}">
</div>
<p>
of coprime pairs <span class="inline-formula"><img class="img-inline-formula img-formula" width="41" height="18" src="https://math.fontein.de/formulae/p1xOHJn1uc3KaP04NNve3_vm0Zr_NQ5jV8MqiQ.svgz" alt="(x, y)" title="(x, y)"></span> in <span class="inline-formula"><img class="img-inline-formula img-formula" width="45" height="19" src="https://math.fontein.de/formulae/7OB5dmmCc89RV663vm4wxxrsWlZlKGUHzpAuWg.svgz" alt="[1, n]^2" title="[1, n]^2"></span>, and let <span class="inline-formula"><img class="img-inline-formula img-formula" width="75" height="18" src="https://math.fontein.de/formulae/EKPb4Iu6Svbj8cDySnu__2QMpVwwPSSNiR_8dA.svgz" alt="P_n = |S_n|" title="P_n = |S_n|"></span> be their number. It is <a href="https://en.wikipedia.org/wiki/Coprime#Probabilities">well-known</a> that <span class="inline-formula"><img class="img-inline-formula img-formula" width="221" height="22" src="https://math.fontein.de/formulae/alRjB9CZyi0ejj0WbcjGHUEPWwtm9akWd_qmLw.svgz" alt="\lim_{n\to\infty} \frac{P_n}{n^2} = \frac{6}{\pi^2} \approx 0.607927" title="\lim_{n\to\infty} \frac{P_n}{n^2} = \frac{6}{\pi^2} \approx 0.607927"></span>; this was first proven by <a href="https://en.wikipedia.org/wiki/Ernesto_Ces%25C3%25A0ro">Ernesto Cesàro</a>.
</p>
<p>
We are interested in <span class="inline-formula"><img class="img-inline-formula img-formula" width="21" height="15" src="https://math.fontein.de/formulae/Ro04tdQRs1hB1P1Rurtcyti_hAC0M5KBw4sdEw.svgz" alt="P_n" title="P_n"></span> itself. Even though we know that <span class="inline-formula"><img class="img-inline-formula img-formula" width="120" height="24" src="https://math.fontein.de/formulae/VGfRvr7XoRwDKAUK.NexqmQQFlnJthqYdtMGNA.svgz" alt="P_n \ge \frac{6 n^2}{\pi^2} - \varepsilon n^2" title="P_n \ge \frac{6 n^2}{\pi^2} - \varepsilon n^2"></span> for some <span class="inline-formula"><img class="img-inline-formula img-formula" width="41" height="12" src="https://math.fontein.de/formulae/ly6Ap7KW4SRUC02VaqDollTKyMXjq8.lKwPx8g.svgz" alt="\varepsilon > 0" title="\varepsilon > 0"></span>, we do not know how this <span class="inline-formula"><img class="img-inline-formula img-formula" width="8" height="8" src="https://math.fontein.de/formulae/XzsXnVu3.wlKUnnE7YgOn8OXuRsBN5WnE9qjKg.svgz" alt="\varepsilon" title="\varepsilon"></span> can be chosen.
</p>
<p>
First note that <span class="inline-formula"><img class="img-inline-formula img-formula" width="105" height="18" src="https://math.fontein.de/formulae/_syo1vvlo9M4jPwuV4VkWbatsgSkSwjLTJ6wVQ.svgz" alt="P_n = 2 P_n' - 1" title="P_n = 2 P_n' - 1"></span>, where
</p>
<div class="display-formula">
<img class="img-display-formula img-formula" width="390" height="19" src="https://math.fontein.de/formulae/ZfX077I_ENFm1rvu8.CyLOBr1rKtpyY5.RAImA.svgz" alt="P_n' = |\{ (x, y) \in \N \mid 1 \le x \le y \le n, \; \gcd(x, y) = 1 \}|," title="P_n' = |\{ (x, y) \in \N \mid 1 \le x \le y \le n, \; \gcd(x, y) = 1 \}|,">
</div>
<p>
and that <span class="inline-formula"><img class="img-inline-formula img-formula" width="127" height="22" src="https://math.fontein.de/formulae/Ja6ki2vok9utpcAqpbgiG1cZHlNRYKFUg6VDfA.svgz" alt="P_n' = \sum_{y=1}^n \varphi(y)" title="P_n' = \sum_{y=1}^n \varphi(y)"></span>, where <span class="inline-formula"><img class="img-inline-formula img-formula" width="12" height="11" src="https://math.fontein.de/formulae/G0SX86eTASn_9M49WF7HzDK85n7NoD6NrqM3ew.svgz" alt="\varphi" title="\varphi"></span> is <a href="https://en.wikipedia.org/wiki/Euler%2527s_totient_function">Euler's <span class="inline-formula"><img class="img-inline-formula img-formula" width="12" height="11" src="https://math.fontein.de/formulae/G0SX86eTASn_9M49WF7HzDK85n7NoD6NrqM3ew.svgz" alt="\varphi" title="\varphi"></span> function</a>. We want to find an explicit lower bound for <span class="inline-formula"><img class="img-inline-formula img-formula" width="21" height="18" src="https://math.fontein.de/formulae/JwK.L_wIJtqPB6QIsXtyDdRZqo2Zw9fDT4LkXw.svgz" alt="P_n'" title="P_n'"></span>, which allows us to find an explicit lower bound for <span class="inline-formula"><img class="img-inline-formula img-formula" width="21" height="15" src="https://math.fontein.de/formulae/Ro04tdQRs1hB1P1Rurtcyti_hAC0M5KBw4sdEw.svgz" alt="P_n" title="P_n"></span> itself.
</p>
<p>
As in the proof of Theorem 330 in <a href="https://en.wikipedia.org/wiki/An_Introduction_to_the_Theory_of_Numbers">Hardy and Wright</a>, we have
</p>
<div class="display-formula">
<img class="img-display-formula img-formula" width="296" height="52" src="https://math.fontein.de/formulae/3CO80_VKpPyLPtr9OV8XhrZpXSj..VomyIipQg.svgz" alt="\sum_{d=1}^n \varphi(d) = \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d} \Bigr\rfloor^2 + \Bigl\lfloor \frac{n}{d} \Bigr\rfloor \biggr)," title="\sum_{d=1}^n \varphi(d) = \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d} \Bigr\rfloor^2 + \Bigl\lfloor \frac{n}{d} \Bigr\rfloor \biggr),">
</div>
<p>
where <span class="inline-formula"><img class="img-inline-formula img-formula" width="11" height="11" src="https://math.fontein.de/formulae/QGPIFfjz60I1ozjZngbI3x6aeMfrAlc1UlQkzg.svgz" alt="\mu" title="\mu"></span> is the <a href="https://en.wikipedia.org/wiki/M%25C3%25B6bius_function">Möbius function</a>. Hardy and Wright compare that expression to <span class="inline-formula"><img class="img-inline-formula img-formula" width="227" height="27" src="https://math.fontein.de/formulae/9v0N6HIQfL_tdcowKiQTJIXQe_7w7H8EFmuTkQ.svgz" alt="\frac{1}{2} \sum_{d=1}^\infty \mu(d) \frac{n^2}{d^2} = \frac{n^2}{2 \zeta(2)} = \frac{3 n^2}{\pi^2}" title="\frac{1}{2} \sum_{d=1}^\infty \mu(d) \frac{n^2}{d^2} = \frac{n^2}{2 \zeta(2)} = \frac{3 n^2}{\pi^2}"></span> and show that the error is <span class="inline-formula"><img class="img-inline-formula img-formula" width="78" height="18" src="https://math.fontein.de/formulae/cn08pcrLZA3eE9bId4jHE6SPb_y1OQRsq3.Rsw.svgz" alt="O(n \log n)" title="O(n \log n)"></span>.
</p>
<p>
More precisely, it is
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="295" height="169" src="https://math.fontein.de/formulae/gWJlMNlAshvn89imjqy_8yrutoN_XiYdgIYihg.svgz" alt="P_n' ={} & \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d} \Bigr\rfloor^2 + \Bigl\lfloor \frac{n}{d} \Bigr\rfloor \biggr) \\
{}={} & \frac{1}{2} \sum_{d=1}^\infty \mu(d) \frac{n^2}{d^2} - \frac{1}{2} \sum_{d=n+1}^\infty \mu(d) \frac{n^2}{d^2} \\
{}+{} & \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d}\Bigr\rfloor^2 + \Bigl\lfloor\frac{n}{d}\Bigr\rfloor - \frac{n^2}{d^2} \biggr)." title="P_n' ={} & \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d} \Bigr\rfloor^2 + \Bigl\lfloor \frac{n}{d} \Bigr\rfloor \biggr) \\
{}={} & \frac{1}{2} \sum_{d=1}^\infty \mu(d) \frac{n^2}{d^2} - \frac{1}{2} \sum_{d=n+1}^\infty \mu(d) \frac{n^2}{d^2} \\
{}+{} & \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor \frac{n}{d}\Bigr\rfloor^2 + \Bigl\lfloor\frac{n}{d}\Bigr\rfloor - \frac{n^2}{d^2} \biggr).">
</div>
<p>
First,
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="389" height="98" src="https://math.fontein.de/formulae/4BQnwMx0unSL7Gow5AbEnMjDVCHDyRH6qgqs9w.svgz" alt="\biggl| \frac{1}{2} \sum_{d=n+1}^\infty \mu(d) \frac{n^2}{d^2} \biggr| \le{} & \frac{n^2}{2} \sum_{d=n+1}^\infty \frac{1}{d^2} \le \frac{n^2}{2} \int_n^\infty x^{-2} \; dx \\
{}={} & \frac{n^2}{2} \Bigl[ -\tfrac{1}{3} x^{-3} \Bigr]_n^\infty = \frac{1}{6 n}." title="\biggl| \frac{1}{2} \sum_{d=n+1}^\infty \mu(d) \frac{n^2}{d^2} \biggr| \le{} & \frac{n^2}{2} \sum_{d=n+1}^\infty \frac{1}{d^2} \le \frac{n^2}{2} \int_n^\infty x^{-2} \; dx \\
{}={} & \frac{n^2}{2} \Bigl[ -\tfrac{1}{3} x^{-3} \Bigr]_n^\infty = \frac{1}{6 n}.">
</div>
<p>
Next,
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="499" height="227" src="https://math.fontein.de/formulae/NQByPkuAJv_8oEZVqgDzbJ7dlV5nuQuniAny5w.svgz" alt="& \biggl| \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor\frac{n}{d}\Bigr\rfloor^2 + \Bigl\lfloor\frac{n}{d}\Bigr\rfloor - \frac{n^2}{d^2} \biggr) \biggr| \\
{}\le{} & \frac{1}{2} \sum_{d=1}^n \biggl| \frac{(n - (n \mod d))^2}{d^2} + \frac{n - (n \mod d)}{d} - \frac{n^2}{d^2} \biggr| \\
{}={} & \frac{1}{2} \sum_{d=1}^n \biggl| \frac{-2 n (n \mod d) + (n \mod d)^2}{d^2} + \frac{n - (n \mod d)}{d} \biggr| \\
{}\le{} & \frac{1}{2} \sum_{d=1}^n \Bigl( \frac{2 n d}{d^2} + \frac{d^2}{d^2} + \frac{n}{d} + \frac{d}{d} \Bigr) = \frac{1}{2} \sum_{d=1}^n \Bigl( \frac{3 n}{d} + 2 \Bigr) = \frac{3}{2} n \sum_{d=1}^n \frac{1}{d} + n." title="& \biggl| \frac{1}{2} \sum_{d=1}^n \mu(d) \biggl( \Bigl\lfloor\frac{n}{d}\Bigr\rfloor^2 + \Bigl\lfloor\frac{n}{d}\Bigr\rfloor - \frac{n^2}{d^2} \biggr) \biggr| \\
{}\le{} & \frac{1}{2} \sum_{d=1}^n \biggl| \frac{(n - (n \mod d))^2}{d^2} + \frac{n - (n \mod d)}{d} - \frac{n^2}{d^2} \biggr| \\
{}={} & \frac{1}{2} \sum_{d=1}^n \biggl| \frac{-2 n (n \mod d) + (n \mod d)^2}{d^2} + \frac{n - (n \mod d)}{d} \biggr| \\
{}\le{} & \frac{1}{2} \sum_{d=1}^n \Bigl( \frac{2 n d}{d^2} + \frac{d^2}{d^2} + \frac{n}{d} + \frac{d}{d} \Bigr) = \frac{1}{2} \sum_{d=1}^n \Bigl( \frac{3 n}{d} + 2 \Bigr) = \frac{3}{2} n \sum_{d=1}^n \frac{1}{d} + n.">
</div>
<p>
Since
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="260" height="49" src="https://math.fontein.de/formulae/eM_QliDET0WC1IPXP9JFFS6iKSu72F143GI8qQ.svgz" alt="\sum_{d=1}^n \frac{1}{d} \le 1 + \int_1^n \frac{1}{x} \; dx = 1 + \log n," title="\sum_{d=1}^n \frac{1}{d} \le 1 + \int_1^n \frac{1}{x} \; dx = 1 + \log n,">
</div>
<p>
we finally obtain
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="479" height="39" src="https://math.fontein.de/formulae/ATq64I4jlbL8eVY6pQXsxPxD4O95rI0zuW2E1g.svgz" alt="P_n' \ge \frac{3 n^2}{\pi^2} - \frac{1}{6 n} - \frac{3}{2} n (1 + \log n) - n = \frac{3 }{\pi^2} n^2 - \frac{8}{3} n - \frac{3}{2} n \log n." title="P_n' \ge \frac{3 n^2}{\pi^2} - \frac{1}{6 n} - \frac{3}{2} n (1 + \log n) - n = \frac{3 }{\pi^2} n^2 - \frac{8}{3} n - \frac{3}{2} n \log n.">
</div>
<p>
With a simple computer program, I verified that
</p>
<div class="display-formula">
<img class="img-display-formula img-formula" width="401" height="36" src="https://math.fontein.de/formulae/T9zKRcopOUSAuvBKlSxeUmFo..2WGsTLDXc8dw.svgz" alt="\frac{P_n}{n^2} \ge \frac{2 \cdot 494866 - 1}{1276^2} = \frac{989731}{1628176} \approx 0.6078771582433" title="\frac{P_n}{n^2} \ge \frac{2 \cdot 494866 - 1}{1276^2} = \frac{989731}{1628176} \approx 0.6078771582433">
</div>
<p>
for all <span class="inline-formula"><img class="img-inline-formula img-formula" width="87" height="12" src="https://math.fontein.de/formulae/sH0BbXd3MP8k8_K_CG38.rV8Q51ScxwwwWF8UA.svgz" alt="n < 932453" title="n < 932453"></span>, and the above bound shows
</p>
<div class="align-formula">
<img class="img-align-formula img-formula" width="411" height="102" src="https://math.fontein.de/formulae/iGYmK1C7uLi.Ysj9hkkeyFcdxgME2QezULVtXQ.svgz" alt="P_n ={} & 2 P_n' - 1 \ge \frac{6}{\pi^2} n^2 - \frac{16}{3} n - 3 n \log n - 1 \\
{}\ge{} & \Bigl( \frac{6}{\pi^2} - \frac{16}{3 \cdot 932453} - \frac{3 \log 932453}{932453} - \frac{1}{932453^2} \Bigr) n^2 \\
{}>{} & 0.607877158257419 n^2" title="P_n ={} & 2 P_n' - 1 \ge \frac{6}{\pi^2} n^2 - \frac{16}{3} n - 3 n \log n - 1 \\
{}\ge{} & \Bigl( \frac{6}{\pi^2} - \frac{16}{3 \cdot 932453} - \frac{3 \log 932453}{932453} - \frac{1}{932453^2} \Bigr) n^2 \\
{}>{} & 0.607877158257419 n^2">
</div>
<p>
for <span class="inline-formula"><img class="img-inline-formula img-formula" width="87" height="14" src="https://math.fontein.de/formulae/SNEBf9hGuV4jgubX1BpS3vVkrlg6h_gRe_ejVA.svgz" alt="n \ge 932453" title="n \ge 932453"></span>. The computer verification took less than one second. We therefore obtain:
</p>
<div class="theorem-environment theorem-theorem-environment qed">
<div class="theorem-header theorem-theorem-header">
Theorem.
</div>
<div class="theorem-content theorem-theorem-content">
<p>
We have <span class="inline-formula"><img class="img-inline-formula img-formula" width="288" height="21" src="https://math.fontein.de/formulae/nhOGPMEfb0tsaERVRKSimApsJ.Bzz2Pwos0kWA.svgz" alt="P_n \ge \frac{989731}{1628176} n^2 > 0.6078771582433 n^2" title="P_n \ge \frac{989731}{1628176} n^2 > 0.6078771582433 n^2"></span> for all <span class="inline-formula"><img class="img-inline-formula img-formula" width="45" height="13" src="https://math.fontein.de/formulae/UeChwhh6h7Bw5vi122OTLvfxRX.k6iU0MpIILg.svgz" alt="n \in \N" title="n \in \N"></span>. Equality holds if and only if <span class="inline-formula"><img class="img-inline-formula img-formula" width="70" height="11" src="https://math.fontein.de/formulae/uLyG2CMaAgxyR2sYSdHH7KqSoVPIwgR89d_y5g.svgz" alt="n = 1276" title="n = 1276"></span>.
</p>
</div>
<div class="qed-block"><span class="qed-sign"></span></div>
</div>
<p>
The program for verification of the bound is the following C++ program:
</p>
</div>
<div class="code-c++"><pre class="code literal-block"><span></span><span class="linenos"> 1</span><span class="cp">#include</span><span class="w"> </span><span class="cpf"><iostream></span>
<span class="linenos"> 2</span><span class="cp">#include</span><span class="w"> </span><span class="cpf"><iomanip></span>
<span class="linenos"> 3</span><span class="cp">#include</span><span class="w"> </span><span class="cpf"><cmath></span>
<span class="linenos"> 4</span><span class="cp">#include</span><span class="w"> </span><span class="cpf"><vector></span>
<span class="linenos"> 5</span>
<span class="linenos"> 6</span><span class="c1">// Simple implementation of Euclid's algorithm</span>
<span class="linenos"> 7</span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="nf">gcd</span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">a</span><span class="p">,</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">b</span><span class="p">)</span>
<span class="linenos"> 8</span><span class="p">{</span>
<span class="linenos"> 9</span><span class="w"> </span><span class="k">while</span><span class="w"> </span><span class="p">(</span><span class="nb">true</span><span class="p">)</span>
<span class="linenos">10</span><span class="w"> </span><span class="p">{</span>
<span class="linenos">11</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">a</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="linenos">12</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">b</span><span class="p">;</span>
<span class="linenos">13</span><span class="w"> </span><span class="n">b</span><span class="w"> </span><span class="o">%=</span><span class="w"> </span><span class="n">a</span><span class="p">;</span>
<span class="linenos">14</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">b</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="linenos">15</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">a</span><span class="p">;</span>
<span class="linenos">16</span><span class="w"> </span><span class="n">a</span><span class="w"> </span><span class="o">%=</span><span class="w"> </span><span class="n">b</span><span class="p">;</span>
<span class="linenos">17</span><span class="w"> </span><span class="p">}</span>
<span class="linenos">18</span><span class="p">}</span>
<span class="linenos">19</span>
<span class="linenos">20</span><span class="c1">// Compute all primes up to a bound</span>
<span class="linenos">21</span>
<span class="linenos">22</span><span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o"><</span><span class="kt">unsigned</span><span class="o">></span><span class="w"> </span><span class="n">primes</span><span class="p">;</span>
<span class="linenos">23</span>
<span class="linenos">24</span><span class="kt">void</span><span class="w"> </span><span class="nf">create_primes</span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">max</span><span class="p">)</span>
<span class="linenos">25</span><span class="p">{</span>
<span class="linenos">26</span><span class="w"> </span><span class="c1">// TODO: replace with sieve of Eratosthenes</span>
<span class="linenos">27</span><span class="w"> </span><span class="n">primes</span><span class="p">.</span><span class="n">push_back</span><span class="p">(</span><span class="mi">2</span><span class="p">);</span>
<span class="linenos">28</span><span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">p</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">3</span><span class="p">;</span><span class="w"> </span><span class="n">p</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="n">max</span><span class="p">;</span><span class="w"> </span><span class="n">p</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="mi">2</span><span class="p">)</span>
<span class="linenos">29</span><span class="w"> </span><span class="p">{</span>
<span class="linenos">30</span><span class="w"> </span><span class="kt">bool</span><span class="w"> </span><span class="n">prime</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">true</span><span class="p">;</span>
<span class="linenos">31</span><span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">primes</span><span class="p">.</span><span class="n">size</span><span class="p">();</span><span class="w"> </span><span class="o">++</span><span class="n">i</span><span class="p">)</span>
<span class="linenos">32</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">((</span><span class="n">p</span><span class="w"> </span><span class="o">%</span><span class="w"> </span><span class="n">primes</span><span class="p">[</span><span class="n">i</span><span class="p">])</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="linenos">33</span><span class="w"> </span><span class="p">{</span>
<span class="linenos">34</span><span class="w"> </span><span class="n">prime</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">false</span><span class="p">;</span>
<span class="linenos">35</span><span class="w"> </span><span class="k">break</span><span class="p">;</span>
<span class="linenos">36</span><span class="w"> </span><span class="p">}</span>
<span class="linenos">37</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">prime</span><span class="p">)</span>
<span class="linenos">38</span><span class="w"> </span><span class="n">primes</span><span class="p">.</span><span class="n">push_back</span><span class="p">(</span><span class="n">p</span><span class="p">);</span>
<span class="linenos">39</span><span class="w"> </span><span class="p">}</span>
<span class="linenos">40</span><span class="p">}</span>
<span class="linenos">41</span>
<span class="linenos">42</span><span class="c1">// Compute the Euler phi function (for not too large arguments)</span>
<span class="linenos">43</span><span class="kt">unsigned</span><span class="w"> </span><span class="nf">phi</span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">n</span><span class="p">)</span>
<span class="linenos">44</span><span class="p">{</span>
<span class="linenos">45</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">n</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="mi">2</span><span class="p">)</span>
<span class="linenos">46</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="linenos">47</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">n</span><span class="p">;</span>
<span class="linenos">48</span><span class="w"> </span><span class="kt">bool</span><span class="w"> </span><span class="n">found_divisor</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">false</span><span class="p">;</span>
<span class="linenos">49</span><span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">0</span><span class="p">;</span><span class="w"> </span><span class="n">i</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">primes</span><span class="p">.</span><span class="n">size</span><span class="p">();</span><span class="w"> </span><span class="o">++</span><span class="n">i</span><span class="p">)</span>
<span class="linenos">50</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">((</span><span class="n">n</span><span class="w"> </span><span class="o">%</span><span class="w"> </span><span class="n">primes</span><span class="p">[</span><span class="n">i</span><span class="p">])</span><span class="w"> </span><span class="o">==</span><span class="w"> </span><span class="mi">0</span><span class="p">)</span>
<span class="linenos">51</span><span class="w"> </span><span class="p">{</span>
<span class="linenos">52</span><span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">r</span><span class="w"> </span><span class="o">/</span><span class="w"> </span><span class="n">primes</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="p">(</span><span class="n">primes</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="mi">1</span><span class="p">);</span>
<span class="linenos">53</span><span class="w"> </span><span class="n">found_divisor</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">true</span><span class="p">;</span>
<span class="linenos">54</span><span class="w"> </span><span class="p">}</span>
<span class="linenos">55</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="o">!</span><span class="n">found_divisor</span><span class="p">)</span>
<span class="linenos">56</span><span class="w"> </span><span class="o">--</span><span class="n">r</span><span class="p">;</span>
<span class="linenos">57</span><span class="w"> </span><span class="k">return</span><span class="w"> </span><span class="n">r</span><span class="p">;</span>
<span class="linenos">58</span><span class="p">}</span>
<span class="linenos">59</span>
<span class="linenos">60</span><span class="kt">int</span><span class="w"> </span><span class="nf">main</span><span class="p">()</span>
<span class="linenos">61</span><span class="p">{</span>
<span class="linenos">62</span><span class="w"> </span><span class="c1">// Compute probabilities up to:</span>
<span class="linenos">63</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="n">N</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">932453</span><span class="p">;</span>
<span class="linenos">64</span><span class="w"> </span><span class="c1">// unsigned N = 10000000;</span>
<span class="linenos">65</span>
<span class="linenos">66</span><span class="w"> </span><span class="c1">// Compute enough primes</span>
<span class="linenos">67</span><span class="w"> </span><span class="n">create_primes</span><span class="p">(</span><span class="n">sqrt</span><span class="p">(</span><span class="n">N</span><span class="p">));</span>
<span class="linenos">68</span><span class="w"> </span><span class="n">std</span><span class="o">::</span><span class="n">cout</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">"Prepared "</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">primes</span><span class="p">.</span><span class="n">size</span><span class="p">()</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">" primes</span><span class="se">\n</span><span class="s">"</span><span class="p">;</span>
<span class="linenos">69</span>
<span class="linenos">70</span><span class="w"> </span><span class="c1">// Compute probabilities</span>
<span class="linenos">71</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">C</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="linenos">72</span><span class="w"> </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="n">n</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">2</span><span class="p">;</span><span class="w"> </span><span class="n">n</span><span class="w"> </span><span class="o"><</span><span class="w"> </span><span class="n">N</span><span class="p">;</span><span class="w"> </span><span class="o">++</span><span class="n">n</span><span class="p">)</span>
<span class="linenos">73</span><span class="w"> </span><span class="p">{</span>
<span class="linenos">74</span><span class="w"> </span><span class="n">C</span><span class="w"> </span><span class="o">+=</span><span class="w"> </span><span class="n">phi</span><span class="p">(</span><span class="n">n</span><span class="p">);</span><span class="w"> </span><span class="c1">// number of pairs (x, y) with 1 <= x <= y <= n such that gcd(x, y) == 1</span>
<span class="linenos">75</span><span class="w"> </span><span class="c1">// Compute rational number representing probability</span>
<span class="linenos">76</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">num</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="mi">2</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">C</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="mi">1</span><span class="p">;</span>
<span class="linenos">77</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">den</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">n</span><span class="p">;</span>
<span class="linenos">78</span><span class="w"> </span><span class="n">den</span><span class="w"> </span><span class="o">*=</span><span class="w"> </span><span class="n">n</span><span class="p">;</span>
<span class="linenos">79</span><span class="w"> </span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="n">g</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">gcd</span><span class="p">(</span><span class="n">num</span><span class="p">,</span><span class="w"> </span><span class="n">den</span><span class="p">);</span>
<span class="linenos">80</span><span class="w"> </span><span class="n">num</span><span class="w"> </span><span class="o">/=</span><span class="w"> </span><span class="n">g</span><span class="p">;</span>
<span class="linenos">81</span><span class="w"> </span><span class="n">den</span><span class="w"> </span><span class="o">/=</span><span class="w"> </span><span class="n">g</span><span class="p">;</span>
<span class="linenos">82</span><span class="w"> </span><span class="c1">// Compute floating point approximation of probability</span>
<span class="linenos">83</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">double</span><span class="w"> </span><span class="n">P</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="p">(</span><span class="kt">long</span><span class="w"> </span><span class="kt">double</span><span class="p">)</span><span class="n">num</span><span class="w"> </span><span class="o">/</span><span class="w"> </span><span class="p">(</span><span class="kt">long</span><span class="w"> </span><span class="kt">double</span><span class="p">)</span><span class="n">den</span><span class="p">;</span>
<span class="linenos">84</span><span class="w"> </span><span class="c1">// Output if probability is <= 0.6079</span>
<span class="linenos">85</span><span class="w"> </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">num</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="p">(</span><span class="kt">unsigned</span><span class="w"> </span><span class="kt">long</span><span class="w"> </span><span class="kt">long</span><span class="p">)((</span><span class="kt">long</span><span class="w"> </span><span class="kt">double</span><span class="p">)</span><span class="mf">0.6079</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">den</span><span class="p">))</span>
<span class="linenos">86</span><span class="w"> </span><span class="n">std</span><span class="o">::</span><span class="n">cout</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">n</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">": "</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">C</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">", "</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">num</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">"/"</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">den</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">" </span><span class="se">\\</span><span class="s">approx "</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">std</span><span class="o">::</span><span class="n">setprecision</span><span class="p">(</span><span class="mi">20</span><span class="p">)</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="n">P</span><span class="w"> </span><span class="o"><<</span><span class="w"> </span><span class="s">"</span><span class="se">\n</span><span class="s">"</span><span class="p">;</span>
<span class="linenos">87</span><span class="w"> </span><span class="p">}</span>
<span class="linenos">88</span><span class="p">}</span>
</pre>
</div> <div>
<p>
On my laptop, it runs in less than one second. For <span class="inline-formula"><img class="img-inline-formula img-formula" width="117" height="12" src="https://math.fontein.de/formulae/aSUoEpwu0SVfTm4AN.Bdq9I_bVqYFFWE1aNrPQ.svgz" alt="N = 10\,000\,000" title="N = 10\,000\,000"></span>, it requires less than 22 seconds.
</p>
</div>
</div>