Articolele lui Cismaru Eric-Dimitrie (Eric_278)

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;
}

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

Vom prezenta o soluție alternativă a problemei xOy, însă mai ineficientă din punct de vedere al timpului de execuție (de complexitate \(\mathcal{O}(\texttt{VMAX})\), unde \(\texttt{VMAX} = 10^6\)).

Introducere

În problema #xOy , se cere determinarea poziției unui punct \(P\) de coordonate întregi, fiind date punctele \(A\), \(B\), \(C\) (de coordonate întregi) și distanțele \(PA^2\), \(PB^2\), \(PC^2\) (exprimate tot prin numere întregi).

Soluție

Pasul 1

În rezolvarea problemei, vom folosi următoarea proprietate clasică:

Relația lui Leibniz

În triunghiul \(\Delta{}ABC\), cu centrul de greutate în \(G\), pentru orice punct \(P \in \mathcal{P}\) din plan are loc egalitatea:
\[ PG^2 = \frac{1}{3}\left(PA^2 + PB^2 + PC^2\right) – \frac{1}{9}\left(AB^2 + BC^2 + CA^2\right). \]
Demonstrație: Deoarece \(G\) este centrul de greutate al \(\Delta{}ABC\), pentru orice punct \(P \in \mathcal{P}\):
\[ \overrightarrow{PG} = \frac{\overrightarrow{PA} + \overrightarrow{PB} + \overrightarrow{PC}}{3}. \]
Prin ridicare la pătrat, vom avea
\begin{align*}
PG^2 = \overrightarrow{PG} \cdot \overrightarrow{PG} &= \frac{1}{9} (\overrightarrow{PA} + \overrightarrow{PB} + \overrightarrow{PC})^2 \\ &= \frac{1}{9} \left(\sum PA^2 + 2 \sum \overrightarrow{PA} \cdot \overrightarrow{PB}\right) \\ &= \frac{1}{9} \left(\sum PA^2 + \sum(PA^2 + PB^2 – AB^2)\right) \\ &= \frac{1}{3} (PA^2 + PB^2 + PC^2) – \frac{1}{9} (AB^2 + BC^2 + CA^2),
\end{align*}
ceea ce încheie demonstrația.

Pasul 2

Deoarece distanțele \(PA^2\), \(PB^2\), \(PC^2\) sunt date, iar
\[PG^2 = \frac{1}{3}\left(PA^2 + PB^2 + PC^2\right) – \frac{1}{9}\left(AB^2 + BC^2 + CA^2\right),\]
deducem că \(PG^2 = \text{const} \Leftrightarrow PG = \text{const}\). Așadar, din punct de vedere geometric, \(P\) se află pe un cerc \(\mathcal{C}\) de centru \(G\) și rază \(\sqrt{k}\), unde am notat
\[k := PG^2 = \frac{1}{3}\left(PA^2 + PB^2 + PC^2\right) – \frac{1}{9}\left(AB^2 + BC^2 + CA^2\right).\]

Pasul 3

Ne folosi acum de restricțiile problemei, anume că toate cele trei distanțe date și punctele \(A\), \(B\), \(C\) au coordonate întregi. Astfel, numărul
\[k’ := 9k = 3(PA^2 + PB^2 + PC^2) – (AB^2 + BC^2 + CA^2) \in \mathbb{N}.\]
De asemenea, centrul de greutate al triunghiului \(\Delta{}ABC\) are coordonatele:
\begin{cases}
x_G = \dfrac{x_A + x_B + x_C}{3} \in \mathbb{Q}, \\
y_G = \dfrac{y_A + y_B + y_C}{3} \in \mathbb{Q}.
\end{cases}
Apoi \[PG^2 = k \Leftrightarrow (x_P – x_G)^2 + (y_P – y_G)^2 = k \Leftrightarrow (3x_P – x_{G’})^2 + (3y_P – y_{G’})^2 = k’^2,\]
unde am efectuat notațiile:
\begin{cases}
x_{G’} := 3x_G = x_A + x_B + x_C \in \mathbb{Z}, \\
y_{G’} := 3y_G = y_A + y_B + y_C \in \mathbb{Z}.
\end{cases}
Totodată, luând \(x := 3x_P – x_{G’}\) și \(y := 3y_P – y_{G’}\), avem \(x^2 + y^2 = k’\), cu \(x\), \(y\), \(k’ \in \mathbb{Z}\).

Pasul 4

Vom determina o soluție \((x, y) \in \mathbb{Z} \times \mathbb{Z}\) a ecuației \(x^2 + y^2 = k’^2\). Revenind la notațiile efectuate,
\begin{cases}
x_P = \dfrac{x + x_{G’}}{3} \in \mathbb{Z}, \\
y_P = \dfrac{y + y_{G’}}{3} \in \mathbb{Z}.
\end{cases}
Am determinat așadar un posibil punct laticial \(P(x_P, y_P) \in \mathcal{C}\), însă trebuie să mai verificăm dacă distanțele de la \(P\) la punctele \(A\), \(B\) și \(C\) sunt cele date în problemă:
\[PA^2 = D_{PA}^2, PB^2 = D_{PB}^2, PC^2 = D_{PC}^2.\]
Să observăm că este necesar ca \(x + x_{G’} \equiv 0 (\text{mod } 3)\) și \(y + y_{G’} \equiv 0 (\text{mod } 3)\).

Detalii de implementare

Vom folosi o structură Point, prin care vom memora coordonatele unui punct:

struct Point
{
    int x, y;
};

Distanța dintre două puncte

Pentru două puncte \(p(p_x, p_y)\) și \(q(q_x, q_y)\), definim o funcție Dist(p, q) care să calculeze pătratul distanței dintre acestea, dată de formula:
\[(p_x – q_x)^2 + (p_y – q_y)^2 = dx^2 + dy^2,\]
unde: \begin{cases}
\displaystyle{dx = p_x – q_x}, \\
\displaystyle{dy = p_y – q_y}.
\end{cases}

long long Dist(Point &p, Point &q)
{
    int dx = p.x - q.x,
        dy = p.y - q.y;
    return (long long)dx * dx +
           (long long)dy * dy;
}

Funcția de verificare

De asemenea, vom implementa și o funcție TryPoint(x, y) care verifică dacă perechea \((x, y)\) are proprietățile menționate anterior.

void TryPoint(int x, int y)
{
    if((x + G.x) % 3 == 0 && (y + G.y) % 3 == 0)
    {
        P = Point((x + G.x) / 3, (y + G.y) / 3);
        if(Dist(P, A) == PA2 && Dist(P, B) == PB2 && Dist(P, C) == PC2)
        {
            cout << P << '\n';
            exit(0);
        }
    }
}

Găsirea soluției

Vom itera cu x de la \(0\) la \([\sqrt{k}]\) și vom verifica dacă \(k – x^2\) este pătrat perfect. De menționat că va trebui să probăm același lucru și pentru perechile \((-x, -y)\), \((-x, y)\), \((x, -y)\), motiv pentru care vom declara doi vectori dx și dy (având similar cu vectorii de deplasare). La pasul \(0 \leqslant i < 4\), vom valida perechea \(x \cdot \text{dx}[i], y \cdot \text{dy}[i]\).

Cod sursă

#include <iostream>
#include <cstdlib>
#include <cmath>

using namespace std;

const int dx[] = {1, 1, -1, -1},
          dy[] = {1, -1, 1, -1};

struct Point
{
    int x, y;
    Point(int xx = 0, int yy = 0): x(xx), y(yy) {}
    friend istream& operator >> (istream &in, Point &p)
    {
        in >> p.x >> p.y;
        return in;
    }
    friend ostream& operator << (ostream &out, const Point &p)
    {
        out << p.x << ' ' << p.y;
        return out;
    }
};

Point A, B, C, P, G;
int PA2, PB2, PC2;
long long k;

/// Pătratul distanței dintre punctele p și q
long long Dist(Point &p, Point &q)
{
    int dx = p.x - q.x,
        dy = p.y - q.y;
    return (long long)dx * dx +
           (long long)dy * dy;
}

/// Verifică dacă perechea (x, y) convine.
/// Dacă sunt verificate toate condițiile, se afișează
/// coordonatele lui P și se termină execuția programului.
void TryPoint(int x, int y)
{
    if((x + G.x) % 3 == 0 && (y + G.y) % 3 == 0)
    {
        P = Point((x + G.x) / 3, (y + G.y) / 3);
        if(Dist(P, A) == PA2 && Dist(P, B) == PB2 && Dist(P, C) == PC2)
        {
            cout << P << '\n';
            exit(0);
        }
    }
}

int main()
{
    cin >> A >> B >> C
        >> PA2 >> PB2 >> PC2;
    k = 3LL * (PA2 + PB2 + PC2) - (Dist(A, B) + Dist(B, C) + Dist(C, A));
    G = Point(A.x + B.x + C.x, A.y + B.y + C.y);
    for(int x = 0; 1LL * x * x <= k; x++)
    {
        int y = (int)sqrt(k - 1LL * x * x);
        if(1LL * y * y != k - 1LL * x * x)
            continue;
        for(int i = 0; i < 4; i++)
            TryPoint(x * dx[i], y * dy[i]);
    }
    cout << P;
    return 0;
}

Soluția oficială poate fi accesată aici.