Rezolvare PrimeIntreEle1

Rezumat

Articolul de față prezintă două rezolvări alternative ale problemei \(\texttt{PrimeIntreEle1}\), de complexități \(\mathcal{O}(n \log n)\), una bazată pe formula de inversiune Möbius, iar cealaltă utilizând indicatorul lui Euler și proprietățile acestei funcții.

Să considerăm problema #PrimeIntreEle1 :

Se citește un număr natural \(n\), \(n > 1\). Să se determine câte perechi \((a, b)\), \(1 \leqslant a \leqslant b \leqslant n\) de numere naturale sunt prime între ele.

Soluție brute-force

O primă soluție constă în numărarea efectivă a perechilor \((a, b)\) cu \(1 \leqslant a \leqslant b \leqslant n\) de numere naturale prime între ele. Această abordare are complexitatea \(\mathcal{O}(n^2 \log n)\), deoarece calculul celui mai mare divizor comun \(\mathrm{cmmdc}(a, b)\) necesită \(\mathcal{O}(\log n)\) pași (complexitatea algoritmului lui Euclid – pentru descrierea acestui algoritm, vezi [5]). Putem descrie cu ajutorul pseudocodului această idee:
───────────────────────────────
Soluție brute-force (\(\mathcal{O}(n^2 \log n)\)).
───────────────────────────────
\(cnt \gets 0\);

\(\texttt{pentru}\;\;a = \overline{1, n}\;\texttt{ execută}\)
│ \(\texttt{pentru}\;\;b = \overline{a, n}\;\texttt{execută}\)
│ │ \(\texttt{dacă}\;\;\text{cmmdc}(a, b) = 1\;\texttt{atunci}\)
│ │ │ \(cnt \gets cnt + 1\);
│ │ └■
│ └■
└■

\(\texttt{scrie } cnt\);
───────────────────────────────

Implementare:

#include <iostream>

using namespace std;

int Cmmdc(int x, int y)
{
    int r;
    while(y != 0)
    {
        r = x % y;
        x = y;
        y = r;
    }
    return x;
}

int NumarPerechi(int n)
{
    int res = 0;
    for(int a = 1; a <= n; a++)
        for(int b = a; b <= n; b++)
            if(Cmmdc(a, b) == 1)
                res++;
    return res;
}

int main()
{
    int n;
    cin >> n;
    cout << NumarPerechi(n);
    return 0;
}

Soluție îmbunătățită

Vom prezenta acum și o altă soluție, mult mai eficientă, bazată pe formula de inversiune a lui Möbius.

Fie \(N = |\{(a, b) \; | \; a, b \in \{1, 2, \ldots, n\}, \text{cmmdc}(a, b) = 1\}|\).

Răspunsul problemei va fi atunci \(\frac{N + 1}{2}\). Într-adevăr, perechea \((a, b)\) cu \(a < b\) va fi numărată de două ori: \((a, b)\) și \((b, a)\), excepție făcând perechea \((1, 1)\) (perechile \((a, a)\) nu vor fi numărate pentru \(a > 1\)).

Funcția Möbius, \(\mu : \mathbb{N}^{\star} \to \mathbb{Z}\) este definită astfel:
\[
\mu(n) =
\begin{cases}
1, \text{ dacă } n = 1 \\
0, \text{ dacă există un număr prim } p \text{ astfel încât } p^2 | n \\
(-1)^k, \text{ dacă } n = p_1 \cdot \ldots \cdot p_k, \text{ unde } p_1, \ldots, p_k \text{ sunt prime distincte}
\end{cases}
\]

De asemenea, definim și funcția unitate \(\varepsilon : \mathbb{N}^{\star} \to \mathbb{N}\),
\[
\varepsilon(n) =
\begin{cases}
1, \text{ dacă } n = 1 \\
0, \text{ altfel }
\end{cases}
\]

Definiție. O funcție \(f : \mathbb{N}^{\star} \to \mathbb{C}\) se numește funcție aritmetică.

Teoremă (Formula de inversiune Möbius).

Dacă \(f\), \(g\) sunt două funcții aritmetice, atunci \[ f(n) = \sum_{d | n} \mu(d) g\left(\frac{n}{d}\right) \Leftrightarrow g(n) = \sum_{d | n} f(d). \]

Pentru demonstrație, se poate consulta [2].

În particular, luând \(f = \mu\), avem \(\mu(n) = f(n) = \sum_{d | n} \mu(d) \varepsilon\left(\frac{n}{d}\right)\), deci \(g = \varepsilon\), i.e.
\begin{equation}
\sum_{d | n} \mu(d) = \varepsilon(n). \;\;\;\; (1)
\label{eq1}
\end{equation}

Revenind la problemă, avem:
\begin{align*}
N &= |\{(a, b) \; | \; a, b \in \{1, 2, \ldots, n\}, \text{cmmdc}(a, b) = 1\}| \\ &= \sum_{1 \leqslant a, b \leqslant n} \varepsilon(\text{cmmdc}(a, b)) \\ &= \sum_{1 \leqslant a \leqslant n} \sum_{1 \leqslant b \leqslant n} \varepsilon(\text{cmmdc}(a, b)).
\end{align*}

Punând \(n \to \text{cmmdc}(a, b)\) în relația \((1)\), obținem:
\begin{align*}
N &= \sum_{1 \leqslant a \leqslant n} \sum_{1 \leqslant b \leqslant n} \varepsilon(\text{cmmdc}(a, b)) \\ &= \sum_{1 \leqslant a \leqslant n} \sum_{1 \leqslant b \leqslant n} \sum_{d | \text{cmmdc}(a, b)} \mu(d)
\end{align*}

Folosim acum o proprietate cunoscută, anume că \(d | \text{cmmdc}(a, b)\) dacă și numai dacă \(d | a\) și \(d | b\). Deducem că:
\begin{align*}
N &= \sum_{1 \leqslant a \leqslant n} \sum_{1 \leqslant b \leqslant n} \sum_{\substack{d | a \\ d | b}} \mu(d)
\end{align*}

Schimbând ordinea de sumare,
\begin{align*}
N &= \sum_{1 \leqslant a \leqslant n} \sum_{1 \leqslant b \leqslant n} \sum_{\substack{d | a \\ d | b}} \mu(d) \\ &= \sum_{d = 1}^n \mu(d) \sum_{\substack{1 \leqslant a \leqslant n \\ d | a}} \sum_{\substack{1 \leqslant b \leqslant n \\ d | b}} 1 \\ &= \sum_{d = 1}^n \mu(d) \left(\sum_{\substack{1 \leqslant a \leqslant n \\ d | a}} 1\right)^2 \\ &= \sum_{d = 1}^n \mu(d) \left\lfloor\frac{n}{d}\right\rfloor^2.
\end{align*}

Acest lucru este posibil deoarece:
\[\{(a, b, d) | 1 \leqslant a \leqslant n, 1 \leqslant b \leqslant n, d | a, d | b\} = \{(a, b, d) | 1 \leqslant d \leqslant n, 1 \leqslant a \leqslant n, d | a, 1 \leqslant b \leqslant n, d | b\}.\]
De remarcat că am folosit că \[|\{a \;|\; 1 \leqslant a \leqslant n, d | a\}| = |\{k \cdot d \;|\; 1 \leqslant k \cdot d \leqslant n\}| = \left\lfloor\frac{n}{d}\right\rfloor\] pentru un număr \(d \in \{1, 2, \ldots, n\}\) fixat.

În consecință, numărul căutat este
\begin{equation}
N = \sum_{d = 1}^n \mu(d) \left\lfloor\frac{n}{d}\right\rfloor^2. \;\;\;\; (2)
\label{eq2}
\end{equation}

Ecuația \((2)\) ne oferă un mod de calcul pentru \(N\) în \(\mathcal{O}(n)\). În funcție de cum se precalculează funcția \(\mu\), complexitatea variază de la \(\mathcal{O}(n)\) până la \(\mathcal{O}(n \log n)\).

O variantă de calcul a funcției Möbius în timp \(\mathcal{O}(n \log n)\) se bazează pe faptul că \(\mu(1) = 1\) și pe relația \((1)\), ce poate fi rescrisă
\begin{equation}
\mu(n) = \varepsilon(n)\; – \sum_{\substack{d | n \\ d < n}} \mu(d).
\label{eq3}
\end{equation}

Astfel, obținem următoarea schemă pentru calculul funcției \(\mu\):

───────────────────────────────
Precalcularea funcției Möbius \(\mu\)
───────────────────────────────
\(\texttt{mobius}(1) \gets 1\);
\(\texttt{pentru}\;\;i = \overline{1, n}\;\texttt{execută}\)
│ \(\texttt{pentru}\;\;j = \overline{2i, n, i}\;\texttt{execută}\)
│ │ \(\texttt{mobius}(j) \gets \texttt{mobius}(j)\; – \;\texttt{mobius}(i)\);
│ └■
└■
───────────────────────────────

Pentru exemplificare, să luăm cazul \(n = 6\). Avem \(12\) perechi de numere prime între ele: \((1, 1), (1, 2), \ldots, (1, 6), (2, 3), (2, 5), (3, 4), (3, 5), (4, 5), (5, 6)\).
Într-adevăr, aplicând formula de mai sus,
\begin{align*}
N &= \sum_{d = 1}^{6} \mu(d) \cdot \left\lfloor\frac{6}{d}\right\rfloor^2 \\ &= \mu(1) \cdot \left\lfloor\frac{6}{1}\right\rfloor^2 + \mu(2) \cdot \left\lfloor\frac{6}{2}\right\rfloor^2 + \mu(3) \cdot \left\lfloor\frac{6}{3}\right\rfloor^2 + \mu(4) \cdot \left\lfloor\frac{6}{4}\right\rfloor^2 + \mu(5) \cdot \left\lfloor\frac{6}{5}\right\rfloor^2 + \mu(6) \cdot \left\lfloor\frac{6}{6}\right\rfloor^2 \\ &= 36 – 9 – 4 – 1 + 1 \\ &= 23,
\end{align*}
iar răspunsul este dat de \(\frac{N + 1}{2} = \frac{23 + 1}{2} = 12\). (Avem \(\mu(1) = \mu(6) = 1\), \(\mu(2) = \mu(3) = \mu(5) = -1\), \(\mu(4) = \mu(2^2) = 0\).)

Implementare:

#include <iostream>

using namespace std;

const int NMAX = 1000;

int mobius[NMAX + 1];

void Precalc()
{
    mobius[1] = 1;
    for(int i = 1; i <= NMAX; i++)
        for(int j = i << 1; j <= NMAX; j += i)
            mobius[j] -= mobius[i];
}

int NumarPerechi(int n)
{
    int res = 0;
    for(int d = 1; d <= n; d++)
        res += mobius[d] * (n / d) * (n / d);
    return (res + 1) >> 1;
}

int main()
{
    Precalc();
    int n;
    cin >> n;
    cout << NumarPerechi(n);
    return 0;
}

Soluție cu indicatorul lui Euler

Pentru această problemă, mai există încă o soluție în care intervine indicatorul lui Euler. Pentru a calcula \(N := |\{ (a, b) | 1 \leqslant a \leqslant b \leqslant n, \text{cmmdc}(a, b) = 1 \}|\), putem fixa valoarea \(b \in \{1, 2, \ldots n\}\). Atunci:
\[ N = \sum_{b = 1}^{n} |\{ (a, b) | 1 \leqslant a \leqslant b, \text{cmmdc}(a, b) = 1\}|. \]

Fie \(\varphi : \mathbb{N}^{\star} \to \mathbb{N}^{\star}\), \(\varphi(n) = |\{x | 1 \leqslant x \leqslant n, \text{cmmdc}(x, n) = 1\}|\). Cu alte cuvinte, \(\varphi(n)\) reprezintă numărul de numere prime cu \(n\), mai mici sau egale cu \(n\). (Avem \(\varphi(1) = 1\).) Deducem că:
\[ N = \sum_{x = 1}^n \varphi(x) \]

Pentru precalcularea indicatorului lui Euler, vom folosi o proprietate datorată lui Gauss:
\[ \sum_{d | n} \varphi(d) = n \Leftrightarrow \varphi(n) = n \;-\; \sum_{\substack{d | n \\ d < n}} \varphi(d). \]

───────────────────────────────
Precalcularea funcției \(\varphi\)
───────────────────────────────
\(\texttt{phi}(1) \gets 1\);

\(\texttt{pentru } i = \overline{2, n} \texttt{ execută }\)
│ \(\texttt{phi}(i) \gets i \;-\; 1\);
└■
\(\texttt{pentru } i = \overline{2, n} \texttt{ execută } \)
│ \(\texttt{pentru } j = \overline{2i, n, i} \texttt{ execută }\)
│ │ \(\texttt{phi}(j) \gets \texttt{phi}(j) \;-\; \texttt{phi}(i)\);
│ └■
└■
───────────────────────────────

Complexitatea este tot \(\mathcal{O}(n \log n)\).

Remarcă: Pentru precalcularea funcțiilor \(\mu\) și \(\varphi\), complexitatea este \(\mathcal{O}(n \log n)\), deoarece, pentru \(i \in \{2, 3, \ldots, n\}\), \(j\) poate lua cel mult \(\frac{n}{i}\) valori. În total, vor fi efectuate maxim \(\displaystyle{\sum_{k = 1}^n} \frac{n}{k} = n \cdot H_n\) operații, unde \(H_n = \displaystyle{\sum_{k = 1}^n \frac{1}{k}}\).

Observație. Avem \(H_n = \Theta(\log n)\).

Demonstrație. Vom arăta că are loc următoarea inegalitate: \[\frac{n}{2} + 1 \leqslant H_{2^n} < n + 1, \text{pentru orice } n \in \mathbb{N^{\star}}.\]

Într-adevăr, să observăm că: \[\{1, 2, \ldots, 2^n\} = [1, 2^n] \cap \mathbb{N^{\star}} = \{1\} \cup \bigcup_{k = 1}^n \{2^{k – 1} + 1, \ldots, 2^k\}.\]
De aici rezultă \[H_{2^n} = 1 + \sum_{k = 1}^n \sum_{2^{k – 1} < j \leqslant 2^k} \frac{1}{j}.\]

Dar \(2^{k – 1} < j \leqslant 2^k \Leftrightarrow \frac{1}{2^k} \leqslant \frac{1}{j} < \frac{1}{2^{k – 1}}\) și, prin urmare,
\[\frac{1}{2} = \sum_{2^{k – 1} < j \leqslant 2^k} \frac{1}{2^k} \leqslant \sum_{2^{k – 1} < j \leqslant 2^k} \frac{1}{j} < \sum_{2^{k – 1} < j \leqslant 2^k} \frac{1}{2^{k – 1}} = 1.\]

În consecință, \[\frac{n}{2} = \sum_{k = 1}^n \frac{1}{2} \leqslant \sum_{k = 1}^n \sum_{2^{k – 1} < j \leqslant 2^k} \frac{1}{j} < \sum_{k = 1}^n 1 = n,\]
adică \(\frac{n}{2} + 1 \leqslant H_{2^n} < n + 1\), pentru orice \(n \in \mathbb{N^{\star}}\).

Totodată,
\begin{align*}
H_n &= \sum_{k = 1}^n \frac{1}{k} \\ &= \sum_{k = 1}^{n – 1} \frac{1}{k} + \frac{1}{n} \\ &= H_{n – 1} + \frac{1}{n}
\end{align*}
Din această egalitate reiese că \(H_n > H_{n – 1}\), pentru orice \(n \geqslant 2\), adică șirul \((H_n)_{n \geqslant 1}\) este monoton strict crescător.

Fie acum un număr natural \(n \in \mathbb{N^{\star}}\) arbitrar, fixat. Notăm \(k := \lfloor \log_2 n \rfloor\), deci \(k \leqslant \log_2 n < k + 1 \Leftrightarrow 2^k \leqslant n < 2^{k + 1}\).

Din monotonia șirului \((H_n)_{n \geqslant 1}\), obținem \(H_{2^k} \leqslant H_n < H_{2^{k + 1}}\). Conform inegalității precedente, deducem că \(H_{2^k} \geqslant \frac{k}{2} + 1\) și \(H_{2^{k + 1}} < k + 2\), deci \(\frac{k}{2} + 1 \leqslant H_n < k + 2\), i.e. \(H_n = \Theta(\log n)\).

Implementare:

#include <iostream>

using namespace std;

const int NMAX = 1000;

int phi[NMAX + 1];

void Precalc()
{
    phi[1] = 1;
    for(int i = 2; i <= NMAX; i++)
        phi[i] = i - 1;
    for(int i = 2; i <= NMAX; i++)
        for(int j = i << 1; j <= NMAX; j += i)
            phi[j] -= phi[i];
}

int NumarPerechi(int n)
{
    int res = 0;
    for(int i = 1; i <= n; i++)
        res += phi[i];
    return res;
}

int main()
{
    Precalc();
    int n;
    cin >> n;
    cout << NumarPerechi(n);
    return 0;
}

Aplicație

O altă problemă în care se cere calcularea numărului de perechi de numere naturale mai mici decât un număr \(n\) este #divq .

Avem de răspuns la două întrebări (interogări), și anume:

• Tipul 1: O cameră cu suprafața \(N\) trebuie împărțită în zone cu aceeași suprafață, număr natural. În câte moduri se poate face această împărțire?
• Tipul 2: Două camere cu suprafețele \(a\) și \(b\) pot fi vecine doar dacă \(a\) și \(b\) au un divizor comun mai mare decât \(1\). Câte perechi numere naturale din intervalul \([1, N]\) pot fi suprafețele unor camere vecine?

Pentru primul tip de întrebare, se observă că răspunsul căutat este chiar \(\tau(N)\), adică numărul divizorilor naturali ai lui \(N\). Deoarece avem \(Q \leqslant 10^6\) întrebări, va trebui să precalculăm, pentru toate numerele \(\leqslant \texttt{NMAX}\), numărul de divizori al acestora. Notăm cu \[\texttt{nrDiv}[n] = \tau(n) \text{ (numărul de divizori ai lui } n \text{)}.\]
Vom proceda asemănător cu algoritmul de la ciurul lui Eratostene. Parcurgem cu \(i\) toate numerele de la \(1\) la \(n\) și cu \(j\) toți multiplii lui \(i\) – de forma \(j = k \cdot i, k \in \mathbb{N}^{\star}\) – mai mici decât \(n\). Pentru fiecare astfel de număr \(j\), știm că \(i | j\), prin urmare incrementăm numărul de divizori al lui \(j\).

───────────────────────────────
Precalcularea funcției \(\tau\)
───────────────────────────────
\(\texttt{pentru } i = \overline{1, n} \texttt{ execută } \)
│ \(\texttt{pentru } j = \overline{i, n, i} \texttt{ execută }\)
│ │ \(\texttt{nrDiv}[j] \gets \texttt{nrDiv}[j] + 1\);
│ └■
└■
───────────────────────────────

Pentru cel de-al doilea tip de întrebare, numărul căutat este
\[\frac{N \cdot (N – 1)}{2} \;-\; \sum_{x = 1}^{N} \varphi(x) \;+\; 1.\]
Altfel spus, din numărul total de perechi de numere \((a, b)\), egal cu \(\binom{N}{2} = \frac{N \cdot (N – 1)}{2}\), scădem numărul tuturor perechilor coprime (vezi problema precedentă) și adăugăm \(1\) (deoarece nu ne interesează perechea \((1, 1)\) care ar fi contorizată prin \(\varphi(1)\) din suma de mai sus).

Codul sursă.

#include <fstream>

using namespace std;

ifstream f("divq.in");
ofstream g("divq.out");

const int NMAX = 1500000;

int nrDiv[NMAX + 1];
long long phi[NMAX + 1];
int Q, T, N;

void CalcNrDiv()
{
    for(int i = 1; i <= NMAX; i++)
        for(int j = i; j <= NMAX; j += i)
            nrDiv[j]++; /// j = k * i este multiplu de i
}

void CalcPhi()
{
    for(int i = 1; i <= NMAX; i++)
        phi[i] = i;
    for(int i = 1; i <= NMAX; i++)
        for(int j = i << 1; j <= NMAX; j += i)
            phi[j] -= phi[i];
    for(int i = 1; i <= NMAX; i++)
        phi[i] += phi[i - 1]; /// sume partiale
}

int main()
{
    CalcNrDiv();
    CalcPhi();
    f >> Q;
    while(Q--)
    {
        f >> T >> N;
        if(T == 1)
            g << nrDiv[N];
        else
            g << (long long)N * (N - 1) / 2 - phi[N] + 1;
        g << '\n';
    }
    f.close();
    g.close();
    return 0;
}

Bibliografie

[1] Math note — Möbius inversion
[2] Titu Andreescu, Dorin Andrica, NUMBER THEORY – Structures, Examples, and Problems
[3] Eratostene și alte ciururi
[4] Secvențe Farey
[5] CMMDC și CMMMC. Algoritmul lui Euclid