Let
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ is squarefree}\right\}.$$
Problem 864 asks for \(C(123567101113)\). The difficulty is that squarefreeness is a global condition: for every prime \(p\), we must avoid \(p^2\mid x^2+1\). The implementations therefore replace direct testing of each \(x\) by an inclusion-exclusion sum over square divisors, then split that sum into a small-divisor range handled by congruence counting and a large-divisor range handled by a negative Pell equation.
The key object is the indicator of squarefreeness. Once that indicator is expanded arithmetically, the problem becomes counting solutions to \(x^2\equiv -1\) modulo squares.
For any positive integer \(m\), the standard identity
$$\mathbf{1}_{m\text{ squarefree}}=\sum_{d^2\mid m}\mu(d)$$
holds, where \(\mu\) is the Möbius function. Applying this to \(m=x^2+1\) gives
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
with
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
So the whole problem reduces to understanding which \(d\) admit roots of \(x^2\equiv -1\pmod{d^2}\), and how many such roots there are.
If an odd prime \(p\) divides \(x^2+1\), then \(-1\) must be a quadratic residue modulo \(p\), which happens exactly for \(p\equiv 1\pmod 4\). Therefore any prime \(p\equiv 3\pmod 4\) can never appear inside a contributing divisor \(d\). Also, \(\mu(d)=0\) whenever \(d\) is not squarefree, so only squarefree products of primes \(p\equiv 1\pmod 4\) matter.
For one such prime, choose a root \(r_0\) with
$$r_0^2\equiv -1\pmod p.$$
Write \(r_0^2+1=mp\), and look for a lift \(r=r_0+kp\). Then
$$r^2+1\equiv r_0^2+1+2r_0kp \equiv p\left(m+2r_0k\right)\pmod{p^2}.$$
Thus we need
$$m+2r_0k\equiv 0\pmod p,$$
so
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
This produces one root modulo \(p^2\), and its negative gives the second one. Hence every prime \(p\equiv 1\pmod 4\) contributes exactly two roots modulo \(p^2\).
If
$$d=\prod_{i=1}^{t}p_i$$
is squarefree with all \(p_i\equiv 1\pmod 4\), then the Chinese remainder theorem combines the independent choices at each prime. Therefore \(x^2\equiv -1\pmod{d^2}\) has exactly
$$2^{\omega(d)}$$
solutions modulo \(d^2\), where \(\omega(d)\) is the number of distinct prime factors of \(d\).
Fix a split point \(B\). For squarefree \(d\le B\) made only of primes \(p\equiv 1\pmod 4\), let \(\mathcal{R}_d\) be the set of roots modulo \(d^2\). Then
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
Since \(\#\mathcal{R}_d=2^{\omega(d)}\), each admissible \(d\) contributes a completely explicit term \(\mu(d)A_d(n)\). The small-divisor part of the computation is therefore
$$\sum_{\substack{d\le B\\ d\text{ squarefree}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
The implementations enumerate these \(d\) recursively. Every time a new prime is appended, the modulus changes from \(d^2\) to \((dp)^2\), and each existing root splits into two new roots by CRT. The Möbius sign alternates with the number of prime factors, so the recursion itself is exactly the inclusion-exclusion expansion.
For the remaining terms with \(d>B\), congruence enumeration becomes wasteful. Instead, write
$$x^2+1=d^2k$$
and rename \(y=d\). Then the condition becomes
$$x^2-ky^2=-1,$$
with
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
So every large square divisor turns into a solution of a negative Pell equation. The parameter \(k\) is small because \(y\) is large.
There is also a sharp arithmetic filter on \(k\). If an odd prime \(q\equiv 3\pmod 4\) divided \(k\), then \(q\mid x^2+1\), which is impossible. Hence every admissible \(k\) avoids primes \(q\equiv 3\pmod 4\). Perfect-square values of \(k\) are excluded as well, because the Pell equation must be genuinely irrational.
For a fixed admissible nonsquare \(k\), either the equation
$$x^2-ky^2=-1$$
has no solution, or it has a least positive solution \((x_0,y_0)\). The continued-fraction expansion of \(\sqrt{k}\) detects which case occurs and supplies that minimal solution.
From \((x_0,y_0)\), form
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
Then \(u^2-kv^2=1\), so multiplying by \(u+v\sqrt{k}\) preserves the equation \(x^2-ky^2=-1\). All positive solutions are generated by
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
which is equivalent to
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
Whenever \(y_m>B\), that solution contributes \(\mu(y_m)\) to the total. If \(y_m\) is not squarefree, then \(\mu(y_m)=0\), so it contributes nothing. This exactly matches the original inclusion-exclusion sum over large divisors \(d\).
The checkpoint \(C(10)=9\) is easy to see directly. For \(x=1,\dots,10\), the values of \(x^2+1\) are
$$2,5,10,17,26,37,50,65,82,101.$$
Only \(50\) is not squarefree, because \(50=2\cdot 5^2\). Therefore \(C(10)=9\).
The inclusion-exclusion formula sees the same thing. The term \(d=1\) contributes \(10\). The only nonzero correction comes from \(d=5\), because \(x^2\equiv -1\pmod{25}\) has roots \(x\equiv 7,18\pmod{25}\), and only \(x=7\) lies in \([1,10]\). Thus
$$C(10)=10-\!1=9.$$
The same point also appears in the Pell formulation:
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1.$$
So the failing value \(x=7\) corresponds to the negative Pell equation with \(k=2\) and \(y=5\).
The C++, Python, and Java implementations all follow the same two-part plan. First they sieve primes up to the split point and keep only primes \(p\equiv 1\pmod 4\). For each such prime they compute the two lifted roots of \(x^2\equiv -1\pmod{p^2}\). The small-divisor phase then enumerates squarefree products of these primes, combines the root sets with the Chinese remainder theorem, counts how many roots fall in \([1,n]\), and adds the result with the correct Möbius sign. The C++ and Java versions parallelize this enumeration over independent starting branches; the Python version performs the same logic serially.
The second phase sieves admissible values of \(k\) by removing multiples of primes \(q\equiv 3\pmod 4\). For each remaining nonsquare \(k\), the implementation uses continued fractions of \(\sqrt{k}\) to test solvability of the negative Pell equation and to recover the fundamental solution when one exists. It then generates the whole solution chain until \(x>n\). Each time the corresponding \(y\) exceeds the split point, the implementation evaluates \(\mu(y)\) by trial division and adds that weight.
Finally, the base term \(d=1\) contributes \(n\) immediately, and the program checks the machinery against the known values \(C(10)=9\) and \(C(1000)=895\) before evaluating the full target.
Let
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Sieving primes up to \(B\) and admissible \(k\) values up to \(K\) costs \(O(B\log\log B+K\log\log K)\) time and \(O(B+K)\) memory. The small-divisor phase is proportional to the number of admissible squarefree products \(d\le B\) together with their root sets; every extra prime doubles the number of roots, but the product itself also grows quickly, so the recursion depth stays small in practice. The large-divisor phase processes admissible nonsquare \(k\le K\); for each one, the continued-fraction search is short and the Pell solutions grow exponentially, so only a modest number of solutions survive below \(x\le n\). Memory usage is dominated by the sieve arrays and the current root buffers rather than by the Pell phase itself.
Wir definieren
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ ist quadratfrei}\right\}.$$
In Problem 864 ist \(C(123567101113)\) gesucht. Die Schwierigkeit liegt darin, dass Quadratfreiheit gleichzeitig alle Primquadrate betrifft. Deshalb ersetzen die Implementierungen den direkten Test jedes einzelnen \(x\) durch eine Inklusions-Exklusions-Summe über quadratische Teiler und teilen diese Summe anschließend in einen kleinen Bereich für Kongruenzen und einen großen Bereich für eine negative Pell-Gleichung auf.
Der Ausgangspunkt ist die arithmetische Kennfunktion der Quadratfreiheit. Danach bleibt nur noch das Zählen von Lösungen der Kongruenz \(x^2\equiv -1\) modulo Quadraten.
Für jede positive ganze Zahl \(m\) gilt
$$\mathbf{1}_{m\text{ quadratfrei}}=\sum_{d^2\mid m}\mu(d).$$
Setzt man \(m=x^2+1\), so folgt
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
wobei
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
Damit ist das Problem auf zwei Fragen reduziert: Für welche \(d\) existieren überhaupt Lösungen, und wie viele Lösungen gibt es modulo \(d^2\)?
Teilt ein ungerades Primelement \(p\) den Ausdruck \(x^2+1\), dann muss \(-1\) ein quadratischer Rest modulo \(p\) sein. Das geschieht genau für \(p\equiv 1\pmod 4\). Also können Primzahlen \(p\equiv 3\pmod 4\) in keinem beitragenden Divisor \(d\) vorkommen. Außerdem verschwindet \(\mu(d)\), sobald \(d\) nicht quadratfrei ist. Es bleiben also nur quadratfreie Produkte von Primzahlen \(p\equiv 1\pmod 4\).
Für eine solche Primzahl wähle \(r_0\) mit
$$r_0^2\equiv -1\pmod p.$$
Schreibe \(r_0^2+1=mp\) und suche einen Lift \(r=r_0+kp\). Dann gilt
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2},$$
also muss
$$m+2r_0k\equiv 0\pmod p$$
sein, also
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
So erhält man eine Lösung modulo \(p^2\), und ihr Negatives liefert die zweite. Jede Primzahl \(p\equiv 1\pmod 4\) trägt also genau zwei Lösungen modulo \(p^2\) bei.
Ist nun
$$d=\prod_{i=1}^{t}p_i$$
quadratfrei mit \(p_i\equiv 1\pmod 4\), dann kombiniert der chinesische Restsatz die lokalen Entscheidungen unabhängig. Daher besitzt \(x^2\equiv -1\pmod{d^2}\) genau
$$2^{\omega(d)}$$
Lösungen modulo \(d^2\).
Sei \(B\) ein Trennwert. Für quadratfreie \(d\le B\), die nur aus Primzahlen \(p\equiv 1\pmod 4\) bestehen, bezeichne \(\mathcal{R}_d\) die Lösungsmenge modulo \(d^2\). Dann ist
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
Da \(\#\mathcal{R}_d=2^{\omega(d)}\) gilt, ist jeder Summand \(\mu(d)A_d(n)\) vollständig explizit. Der kleine Teil der Summe lautet also
$$\sum_{\substack{d\le B\\ d\text{ quadratfrei}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
Die Implementierungen erzeugen diese \(d\) rekursiv. Fügt man eine neue Primzahl hinzu, wächst der Modulus von \(d^2\) auf \((dp)^2\), und jede vorhandene Restklasse verzweigt durch den chinesischen Restsatz in zwei neue. Das wechselnde Vorzeichen ist genau das Möbius-Vorzeichen der Inklusions-Exklusion.
Für \(d>B\) wird die Kongruenzsicht unpraktisch. Stattdessen schreibt man
$$x^2+1=d^2k$$
und setzt \(y=d\). Dann erhält man
$$x^2-ky^2=-1,$$
mit
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Jeder große quadratische Teiler entspricht also einer Lösung einer negativen Pell-Gleichung. Der Parameter \(k\) bleibt klein, weil \(y\) groß ist.
Zusätzlich gibt es einen starken Filter: Würde eine ungerade Primzahl \(q\equiv 3\pmod 4\) den Wert \(k\) teilen, dann müsste \(q\mid x^2+1\) gelten, was unmöglich ist. Zulässige \(k\) dürfen also keine solchen Primfaktoren besitzen. Quadratische \(k\) werden ebenfalls ausgeschlossen, da die Pell-Gleichung echt irrational sein muss.
Für ein festes zulässiges, nichtquadratisches \(k\) ist die Gleichung
$$x^2-ky^2=-1$$
entweder unlösbar oder besitzt eine kleinste positive Lösung \((x_0,y_0)\). Die Kettenbruchentwicklung von \(\sqrt{k}\) entscheidet dies und liefert im lösbaren Fall die Fundamentallösung.
Aus \((x_0,y_0)\) bildet man
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
Dann erfüllt \(u^2-kv^2=1\), und die Multiplikation mit \(u+v\sqrt{k}\) erhält die Form \(x^2-ky^2=-1\). Alle positiven Lösungen entstehen durch
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
also durch
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
Sobald \(y_m>B\) ist, trägt diese Lösung den Wert \(\mu(y_m)\) zur Summe bei. Ist \(y_m\) nicht quadratfrei, dann ist \(\mu(y_m)=0\), und es gibt keinen Beitrag. Damit ist der große Bereich exakt dieselbe Inklusions-Exklusions-Summe wie zuvor, nur anders parametrisiert.
Die Testgröße \(C(10)=9\) sieht man sofort. Für \(x=1,\dots,10\) erhält man
$$2,5,10,17,26,37,50,65,82,101.$$
Nur \(50\) ist nicht quadratfrei, denn \(50=2\cdot 5^2\). Also ist \(C(10)=9\).
Die Inklusions-Exklusions-Formel beschreibt genau dasselbe. Der Term \(d=1\) liefert \(10\). Die einzige nichttriviale Korrektur kommt von \(d=5\), denn \(x^2\equiv -1\pmod{25}\) hat die Lösungen \(x\equiv 7,18\pmod{25}\), und im Bereich \([1,10]\) liegt nur \(x=7\). Daher
$$C(10)=10-\!1=9.$$
Dasselbe Beispiel erscheint in der Pell-Sicht als
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1.$$
Die C++-, Python- und Java-Implementierungen folgen demselben Zweiphasenplan. Zuerst werden Primzahlen bis zum Trennwert gesiebt, und es bleiben nur \(p\equiv 1\pmod 4\) übrig. Für jede davon werden die beiden gelifteten Lösungen von \(x^2\equiv -1\pmod{p^2}\) berechnet. Im kleinen Bereich werden dann quadratfreie Primprodukte rekursiv aufgebaut, ihre Restklassen mit dem chinesischen Restsatz kombiniert, die Treffer in \([1,n]\) gezählt und mit dem passenden Möbius-Vorzeichen aufsummiert. Die C++- und Java-Versionen parallelisieren diese Enumeration über unabhängige Startäste; die Python-Version rechnet dieselbe Logik seriell.
Im großen Bereich werden zulässige Werte von \(k\) per Sieb erzeugt, indem Vielfache von Primzahlen \(q\equiv 3\pmod 4\) entfernt werden. Für jedes verbleibende nichtquadratische \(k\) verwendet die Implementierung Kettenbrüche von \(\sqrt{k}\), um die Lösbarkeit der negativen Pell-Gleichung zu testen und im positiven Fall die Fundamentallösung zu finden. Anschließend wird die gesamte Lösungsfolge bis \(x>n\) erzeugt. Immer wenn das zugehörige \(y\) größer als der Trennwert ist, wird \(\mu(y)\) per Probeteilung bestimmt und als Gewicht addiert.
Der Basisterm \(d=1\) liefert sofort den Beitrag \(n\). Vor dem eigentlichen Zielwert prüfen die Implementierungen die Methode außerdem an \(C(10)=9\) und \(C(1000)=895\).
Sei
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Das Sieben der Primzahlen bis \(B\) und der zulässigen \(k\)-Werte bis \(K\) kostet \(O(B\log\log B+K\log\log K)\) Zeit und \(O(B+K)\) Speicher. Der kleine Bereich ist proportional zur Anzahl der zulässigen quadratfreien Produkte \(d\le B\) zusammen mit ihren Restklassensystemen; jede zusätzliche Primzahl verdoppelt zwar die Zahl der Wurzeln, aber das Produkt wächst ebenfalls schnell, sodass die Rekursionstiefe praktisch klein bleibt. Der große Bereich durchläuft zulässige nichtquadratische \(k\le K\); die Kettenbruchsuche ist kurz, und die Pell-Lösungen wachsen exponentiell, sodass unter \(x\le n\) nur wenige Folgenglieder auftreten. Der Speicherverbrauch wird vor allem von den Siebtabellen und den aktuellen Wurzelpuffern bestimmt.
Şunu tanımlayalım:
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ karefree'dir}\right\}.$$
Problem 864, \(C(123567101113)\) değerini istiyor. Zorluk şu: karefree olma koşulu bütün asal kareleri aynı anda dışlamayı gerektirir. Bu yüzden uygulamalar her \(x\) için doğrudan test yapmak yerine önce kare bölenler üzerinde dahil et-çıkar kuruyor, sonra da bu toplamı küçük bölenler için kongruans sayımına, büyük bölenler için ise negatif Pell denklemine ayırıyor.
Anahtar nesne karefree göstergesidir. Bu gösterge aritmetik olarak açılınca problem, \(x^2\equiv -1\) denkleminin kare modüller altındaki çözümlerini saymaya dönüşür.
Her pozitif \(m\) için
$$\mathbf{1}_{m\text{ karefree}}=\sum_{d^2\mid m}\mu(d)$$
özdeşliği geçerlidir. Burada \(\mu\) Möbius fonksiyonudur. Bunu \(m=x^2+1\) için uygularsak
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n)$$
elde edilir. Burada
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}$$
olarak tanımlanır. Böylece problem, hangi \(d\) değerleri için çözüm olduğunu ve her \(d^2\) modülünde kaç çözüm bulunduğunu anlamaya iner.
Eğer tek bir asal \(p\), \(x^2+1\)'i bölüyorsa, o zaman \(-1\) sayısı mod \(p\) altında kuadratik kalıntı olmalıdır. Bu da tam olarak \(p\equiv 1\pmod 4\) olduğunda mümkündür. Demek ki \(p\equiv 3\pmod 4\) olan asallar katkı veren hiçbir \(d\) içinde yer alamaz. Ayrıca \(d\) karefree değilse zaten \(\mu(d)=0\) olur. Sonuç olarak yalnızca \(p\equiv 1\pmod 4\) asal çarpanlarından oluşan karefree \(d\) değerleri kalır.
Böyle bir asal için
$$r_0^2\equiv -1\pmod p$$
sağlayan bir \(r_0\) seçelim. \(r_0^2+1=mp\) yazıp \(r=r_0+kp\) biçiminde bir lift arayalım. O zaman
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2}$$
olur. Dolayısıyla
$$m+2r_0k\equiv 0\pmod p$$
yani
$$k\equiv -m(2r_0)^{-1}\pmod p$$
gerekir. Böylece mod \(p^2\) altında bir kök elde edilir; negatif değeri ikinci kökü verir. Sonuç olarak her \(p\equiv 1\pmod 4\) asalı mod \(p^2\) altında tam iki kök üretir.
Eğer
$$d=\prod_{i=1}^{t}p_i$$
karefree ise ve tüm \(p_i\equiv 1\pmod 4\) koşulunu sağlıyorsa, Çin kalan teoremi bu yerel seçimleri bağımsız biçimde birleştirir. Bu yüzden \(x^2\equiv -1\pmod{d^2}\) denkleminin çözüm sayısı tam olarak
$$2^{\omega(d)}$$
olur.
Bir ayırma eşiği \(B\) seçelim. Sadece \(p\equiv 1\pmod 4\) asallarından oluşan karefree \(d\le B\) için, mod \(d^2\) kök kümesini \(\mathcal{R}_d\) ile gösterelim. O zaman
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}$$
olur. \(\#\mathcal{R}_d=2^{\omega(d)}\) olduğu için her \(\mu(d)A_d(n)\) terimi doğrudan hesaplanabilir. Küçük bölenler kısmı böylece
$$\sum_{\substack{d\le B\\ d\text{ karefree}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n)$$
şeklindedir. Uygulamalar bu \(d\) değerlerini özyinelemeli olarak üretir. Her yeni asal eklendiğinde modül \(d^2\)'den \((dp)^2\)'ye çıkar ve var olan her kök CRT ile iki yeni köke ayrılır. İşaretin değişmesi de tam olarak Möbius işaretidir; yani recursion doğrudan dahil et-çıkar açılımını yürütür.
\(d>B\) olduğunda kongruans tabanlı tarama verimsizleşir. Bunun yerine
$$x^2+1=d^2k$$
yazıp \(y=d\) diyelim. Böylece
$$x^2-ky^2=-1$$
denklemi ortaya çıkar ve
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor$$
koşulları geçerli olur. Yani büyük kare bölenler, parametresi küçük olan negatif Pell denklemlerinin çözümlerine karşılık gelir.
\(k\) üzerinde de güçlü bir filtre vardır. Eğer tek bir asal \(q\equiv 3\pmod 4\), \(k\)'yı bölseydi, \(q\mid x^2+1\) olurdu; bu mümkün değildir. Demek ki uygun \(k\) değerleri bu tür asal çarpanlar içermez. Ayrıca tam kare \(k\) değerleri de dışarı atılır; aksi halde denklem gerçek anlamda Pell tipi olmaz.
Sabit bir uygun kare olmayan \(k\) için
$$x^2-ky^2=-1$$
denklemi ya hiç çözümsüzdür ya da en küçük pozitif çözümü \((x_0,y_0)\) vardır. \(\sqrt{k}\)'nın sürekli kesir açılımı hangi durumun geçerli olduğunu belirler ve çözümlü durumda temel çifti verir.
Bu temel çözümden
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}$$
tanımlanır. Burada \(u^2-kv^2=1\) olduğundan, \(u+v\sqrt{k}\) ile çarpmak denklemi korur. Tüm pozitif çözümler
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k})$$
eşitliğiyle, yani
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m$$
özyinelemesiyle üretilir. \(y_m>B\) olduğunda bu çözüm toplamda \(\mu(y_m)\) kadar katkı yapar. \(y_m\) karefree değilse \(\mu(y_m)=0\) olur ve katkı yoktur. Böylece büyük bölgedeki toplam, orijinal dahil et-çıkar toplamıyla tam olarak aynı kalır.
Kontrol değeri \(C(10)=9\) kolayca görülür. \(x=1,\dots,10\) için \(x^2+1\) değerleri
$$2,5,10,17,26,37,50,65,82,101$$
şeklindedir. Bunların içinde yalnızca \(50\) karefree değildir; çünkü \(50=2\cdot 5^2\). Demek ki \(C(10)=9\).
Dahil et-çıkar formülü de aynısını söyler. \(d=1\) terimi \(10\) verir. Tek düzeltme \(d=5\)'ten gelir; çünkü \(x^2\equiv -1\pmod{25}\) denkleminin kökleri \(x\equiv 7,18\pmod{25}\)'tir ve \([1,10]\) aralığında yalnızca \(x=7\) vardır. Bu yüzden
$$C(10)=10-\!1=9$$
olur. Aynı durum Pell biçiminde de
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1$$
olarak görünür.
C++, Python ve Java uygulamaları aynı iki aşamalı planı izler. Önce ayırma eşiğine kadar asal süzgeci kurulur ve yalnızca \(p\equiv 1\pmod 4\) olan asallar tutulur. Her biri için \(x^2\equiv -1\pmod{p^2}\) denkleminin iki lift edilmiş kökü hesaplanır. Küçük bölen aşamasında karefree asal çarpımları özyinelemeli biçimde gezilir, kök kümeleri CRT ile birleştirilir, \([1,n]\) aralığındaki uygun sayılar bölüm-kalan formülüyle sayılır ve Möbius işaretiyle toplanır. C++ ve Java sürümleri bağımsız başlangıç dallarını paralel çalıştırır; Python sürümü aynı matematiği seri biçimde uygular.
İkinci aşama, \(q\equiv 3\pmod 4\) asal katlarını eleyerek uygun \(k\) değerlerini süzer. Kalan her kare olmayan \(k\) için uygulama, \(\sqrt{k}\)'nın sürekli kesirini kullanarak negatif Pell denkleminin çözülebilir olup olmadığını anlar ve varsa temel çözümü elde eder. Sonra çözüm zinciri \(x>n\) olana kadar üretilir. İlgili \(y\) değeri eşikten büyük olduğunda \(\mu(y)\) deneme bölmesiyle hesaplanır ve bu ağırlık toplamda kullanılır.
Son olarak \(d=1\) taban terimi doğrudan \(n\) katkısı verir. Uygulamalar tam hedefe geçmeden önce \(C(10)=9\) ve \(C(1000)=895\) kontrollerini de doğrular.
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor$$
olsun. \(B\)'ye kadar asal süzmek ve \(K\)'ye kadar uygun \(k\) değerlerini işaretlemek \(O(B\log\log B+K\log\log K)\) zaman ve \(O(B+K)\) bellek gerektirir. Küçük bölen aşamasının maliyeti, uygun karefree \(d\le B\) değerlerinin sayısına ve bunların kök kümelerine orantılıdır; her yeni asal kök sayısını ikiye katlasa da çarpım hızla büyüdüğü için özyineleme derinliği pratikte küçüktür. Büyük bölen aşaması uygun kare olmayan \(k\le K\) değerlerini işler; sürekli kesir araması kısa sürer ve Pell çözümleri üstel büyüdüğünden \(x\le n\) altında az sayıda çözüm kalır. Belleği asıl tüketen kısım Pell tarafı değil, süzgeç dizileri ile geçici kök tamponlarıdır.
Definimos
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ es libre de cuadrados}\right\}.$$
El problema 864 pide \(C(123567101113)\). La dificultad es que la condición de ser libre de cuadrados obliga a excluir simultáneamente todos los cuadrados primos. Por eso las implementaciones sustituyen la comprobación directa de cada \(x\) por una suma de inclusión-exclusión sobre divisores cuadrados, y luego dividen esa suma en un rango pequeño tratado con congruencias y un rango grande tratado con una ecuación de Pell negativa.
La pieza central es el indicador aritmético de la propiedad libre de cuadrados. Una vez expandido, todo se reduce a contar soluciones de \(x^2\equiv -1\) módulo cuadrados.
Para cualquier entero positivo \(m\) se cumple
$$\mathbf{1}_{m\text{ libre de cuadrados}}=\sum_{d^2\mid m}\mu(d).$$
Aplicando esto a \(m=x^2+1\), obtenemos
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
donde
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
Así, el problema pasa a ser: qué valores de \(d\) admiten raíces y cuántas raíces hay módulo \(d^2\).
Si un primo impar \(p\) divide \(x^2+1\), entonces \(-1\) debe ser un residuo cuadrático módulo \(p\). Eso ocurre exactamente cuando \(p\equiv 1\pmod 4\). Por tanto, ningún primo \(p\equiv 3\pmod 4\) puede aparecer en un divisor \(d\) que contribuya. Además, si \(d\) no es libre de cuadrados, entonces \(\mu(d)=0\). Solo sobreviven los productos libres de cuadrados de primos \(p\equiv 1\pmod 4\).
Para un primo así, escogemos \(r_0\) con
$$r_0^2\equiv -1\pmod p.$$
Escribimos \(r_0^2+1=mp\) y buscamos una elevación \(r=r_0+kp\). Entonces
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2},$$
de modo que necesitamos
$$m+2r_0k\equiv 0\pmod p,$$
es decir,
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
Esto produce una raíz módulo \(p^2\), y su opuesta produce la segunda. En consecuencia, cada primo \(p\equiv 1\pmod 4\) aporta exactamente dos raíces módulo \(p^2\).
Si
$$d=\prod_{i=1}^{t}p_i$$
es libre de cuadrados y todos los \(p_i\equiv 1\pmod 4\), el teorema chino del resto combina las elecciones locales de manera independiente. Por ello, \(x^2\equiv -1\pmod{d^2}\) tiene exactamente
$$2^{\omega(d)}$$
soluciones módulo \(d^2\).
Fijamos un umbral \(B\). Para cada \(d\le B\) libre de cuadrados formado solo por primos \(p\equiv 1\pmod 4\), llamamos \(\mathcal{R}_d\) al conjunto de raíces módulo \(d^2\). Entonces
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
Como \(\#\mathcal{R}_d=2^{\omega(d)}\), cada término \(\mu(d)A_d(n)\) es explícito. La parte pequeña queda
$$\sum_{\substack{d\le B\\ d\text{ libre de cuadrados}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
Las implementaciones generan estos \(d\) de forma recursiva. Cada vez que se añade un nuevo primo, el módulo pasa de \(d^2\) a \((dp)^2\) y cada raíz existente se divide en dos nuevas raíces mediante el teorema chino del resto. El signo alternante es exactamente el signo de Möbius, así que la recursión coincide con la inclusión-exclusión.
Cuando \(d>B\), enumerar congruencias deja de ser eficiente. Entonces escribimos
$$x^2+1=d^2k$$
y renombramos \(y=d\). Eso da
$$x^2-ky^2=-1,$$
con
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Así, cada divisor cuadrado grande corresponde a una solución de una ecuación de Pell negativa, y el parámetro \(k\) es pequeño porque \(y\) es grande.
Además, \(k\) debe pasar un filtro aritmético fuerte. Si un primo impar \(q\equiv 3\pmod 4\) dividiera a \(k\), entonces tendría que dividir a \(x^2+1\), lo cual es imposible. Por eso los \(k\) admisibles no contienen esos primos. También se descartan los \(k\) cuadrados perfectos, porque la ecuación de Pell debe ser genuinamente no trivial.
Para un \(k\) admisible y no cuadrado, la ecuación
$$x^2-ky^2=-1$$
o bien no tiene soluciones o bien tiene una solución positiva mínima \((x_0,y_0)\). La fracción continua de \(\sqrt{k}\) detecta el caso y proporciona la solución fundamental cuando existe.
A partir de \((x_0,y_0)\), definimos
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
Entonces \(u^2-kv^2=1\), y multiplicar por \(u+v\sqrt{k}\) preserva la ecuación \(x^2-ky^2=-1\). Todas las soluciones positivas se generan mediante
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
equivalentemente
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
Cada vez que \(y_m>B\), esa solución aporta \(\mu(y_m)\) al total. Si \(y_m\) no es libre de cuadrados, entonces \(\mu(y_m)=0\) y no aporta nada. De este modo, la parte grande reproduce exactamente la suma original de inclusión-exclusión sobre divisores \(d\).
El valor de control \(C(10)=9\) se ve enseguida. Para \(x=1,\dots,10\), los valores de \(x^2+1\) son
$$2,5,10,17,26,37,50,65,82,101.$$
Solo \(50\) no es libre de cuadrados, porque \(50=2\cdot 5^2\). Por tanto \(C(10)=9\).
La fórmula de inclusión-exclusión dice exactamente lo mismo. El término \(d=1\) aporta \(10\). La única corrección no trivial viene de \(d=5\), porque \(x^2\equiv -1\pmod{25}\) tiene raíces \(x\equiv 7,18\pmod{25}\), y en \([1,10]\) solo aparece \(x=7\). Luego
$$C(10)=10-\!1=9.$$
El mismo punto reaparece en la formulación de Pell:
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1.$$
Las implementaciones en C++, Python y Java siguen el mismo plan en dos fases. Primero hacen una criba de primos hasta el umbral y conservan solo los \(p\equiv 1\pmod 4\). Para cada uno calculan las dos raíces elevadas de \(x^2\equiv -1\pmod{p^2}\). En la fase de divisores pequeños se recorren recursivamente los productos libres de cuadrados, se combinan sus clases residuales con el teorema chino del resto, se cuentan los aciertos en \([1,n]\) con la fórmula cociente-resto y se acumula todo con el signo correcto de Möbius. Las versiones en C++ y Java paralelizan ramas independientes; la versión en Python implementa la misma lógica en serie.
La segunda fase criba los valores admisibles de \(k\) eliminando múltiplos de primos \(q\equiv 3\pmod 4\). Para cada \(k\) no cuadrado que sobrevive, la implementación usa fracciones continuas de \(\sqrt{k}\) para decidir si la ecuación de Pell negativa es resoluble y, en caso afirmativo, recuperar la solución fundamental. Después genera toda la cadena de soluciones hasta que \(x>n\). Cada vez que el correspondiente \(y\) supera el umbral, se calcula \(\mu(y)\) por división de prueba y se añade ese peso.
Finalmente, el término base \(d=1\) aporta inmediatamente \(n\), y las implementaciones verifican el método con \(C(10)=9\) y \(C(1000)=895\) antes de evaluar el objetivo completo.
Sea
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Cribar primos hasta \(B\) y los valores admisibles de \(k\) hasta \(K\) cuesta \(O(B\log\log B+K\log\log K)\) tiempo y \(O(B+K)\) memoria. La fase de divisores pequeños es proporcional al número de productos libres de cuadrados \(d\le B\) y a sus conjuntos de raíces; cada primo extra duplica el número de raíces, pero el producto crece con rapidez, de modo que la profundidad de la recursión sigue siendo pequeña en la práctica. La fase de divisores grandes recorre los \(k\) admisibles y no cuadrados; la búsqueda por fracciones continuas es breve y las soluciones de Pell crecen exponencialmente, así que solo sobreviven pocos términos con \(x\le n\). La memoria está dominada por las cribas y los buffers de raíces, no por la parte de Pell.
定义
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ 是平方自由数}\right\}.$$
第 864 题要求的是 \(C(123567101113)\)。难点在于“平方自由”不是一个局部条件，而是要求对所有素数 \(p\) 都同时避免 \(p^2\mid x^2+1\)。因此实现不能逐个检查每个 \(x\)，而是先把平方自由指标写成 Möbius 反演形式，再把总和分成两部分：较小的平方因子用同余类计数，较大的平方因子改写成负 Pell 方程。
核心在于把“\(x^2+1\) 是平方自由数”转写为一个可以求和的算术指标。这样一来，问题就变成了统计 \(x^2\equiv -1\) 在各种平方模数下的解。
对任意正整数 \(m\)，都有经典恒等式
$$\mathbf{1}_{m\text{ 是平方自由数}}=\sum_{d^2\mid m}\mu(d),$$
其中 \(\mu\) 是 Möbius 函数。把 \(m\) 换成 \(x^2+1\)，得到
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
其中
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
也就是说，整个问题已经变成两个子问题：哪些 \(d\) 会产生非零贡献，以及对每个这样的 \(d\)，模 \(d^2\) 下到底有多少个根。
如果一个奇素数 \(p\) 整除 \(x^2+1\)，那么 \(-1\) 必须是模 \(p\) 的二次剩余。这只会在 \(p\equiv 1\pmod 4\) 时发生。因此任何 \(p\equiv 3\pmod 4\) 的素数都不可能出现在有贡献的 \(d\) 里。另外，只要 \(d\) 不是平方自由数，\(\mu(d)=0\)，于是我们最终只需要考虑由素数 \(p\equiv 1\pmod 4\) 构成的平方自由 \(d\)。
对这样一个素数，先选一个模 \(p\) 的根 \(r_0\)，满足
$$r_0^2\equiv -1\pmod p.$$
把 \(r_0^2+1\) 写成 \(mp\)，然后寻找模 \(p^2\) 的提升 \(r=r_0+kp\)。代入后有
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2}.$$
因此必须满足
$$m+2r_0k\equiv 0\pmod p,$$
也就是
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
这样就得到一个模 \(p^2\) 的根，而它的相反数给出第二个根。所以每个 \(p\equiv 1\pmod 4\) 的素数都恰好贡献两个模 \(p^2\) 的根。
如果
$$d=\prod_{i=1}^{t}p_i$$
是平方自由数，且所有 \(p_i\equiv 1\pmod 4\)，那么这些局部根可以通过中国剩余定理独立组合。因此 \(x^2\equiv -1\pmod{d^2}\) 恰好有
$$2^{\omega(d)}$$
个解，其中 \(\omega(d)\) 表示 \(d\) 的不同素因子个数。
取一个分界点 \(B\)。对每个 \(d\le B\) 的平方自由数，且其素因子都满足 \(p\equiv 1\pmod 4\)，记模 \(d^2\) 的根集为 \(\mathcal{R}_d\)。那么
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
由于 \(\#\mathcal{R}_d=2^{\omega(d)}\)，每一项 \(\mu(d)A_d(n)\) 都可以显式计算。于是小因子部分就是
$$\sum_{\substack{d\le B\\ d\text{ 是平方自由数}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
实现中并不是把所有 \(d\) 预先列出来，而是递归地枚举这些平方自由乘积。每加入一个新素数，模数就从 \(d^2\) 变成 \((dp)^2\)，原来的每个根通过中国剩余定理分裂成两个新根。符号也随着素因子个数奇偶交替变化，这正是 Möbius 反演的包含-排除结构。
当 \(d>B\) 时，继续枚举同余类就不经济了。这时把条件写成
$$x^2+1=d^2k$$
并记 \(y=d\)，就得到
$$x^2-ky^2=-1,$$
同时满足
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
也就是说，所有“大平方因子”都被改写成负 Pell 方程的整数解。由于 \(y\) 很大，所以新的参数 \(k\) 反而比较小。
这里还有一个重要筛选条件：如果某个奇素数 \(q\equiv 3\pmod 4\) 整除 \(k\)，那么它也会整除 \(x^2+1\)，而这不可能发生。因此可行的 \(k\) 一定不含这类素因子。实现还会排除完全平方数 \(k\)，因为 Pell 方程只有在 \(k\) 非平方时才是我们需要的那类情形。
对固定的可行非平方 \(k\)，方程
$$x^2-ky^2=-1$$
要么无解，要么存在最小的正整数解 \((x_0,y_0)\)。\(\sqrt{k}\) 的连分数展开恰好可以检测这一点，并在有解时给出这个最小解。
得到 \((x_0,y_0)\) 后，构造
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
于是 \(u^2-kv^2=1\)。把任意一个负 Pell 解乘上这个单位，就会得到下一个负 Pell 解。因此全部正整数解可以递推为
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
也就是
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
当某个解的 \(y_m>B\) 时，它对应原问题里的一个大平方因子，此时要把 \(\mu(y_m)\) 加入总和。如果 \(y_m\) 不是平方自由数，则 \(\mu(y_m)=0\)，自然不会贡献。这样就和原始的 Möbius 求和完全对应上了。
实现中的检查值 \(C(10)=9\) 很适合作为例子。对 \(x=1,\dots,10\)，有
$$2,5,10,17,26,37,50,65,82,101.$$
其中只有 \(50\) 不是平方自由数，因为 \(50=2\cdot 5^2\)。所以答案确实是 \(9\)。
从包含-排除的角度看，\(d=1\) 先给出 \(10\) 个数。唯一的修正来自 \(d=5\)，因为 \(x^2\equiv -1\pmod{25}\) 的根是 \(x\equiv 7,18\pmod{25}\)，在区间 \([1,10]\) 内只有 \(x=7\)。于是
$$C(10)=10-\!1=9.$$
同一个现象在 Pell 形式下写成
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1,$$
说明 \(x=7\) 这个失败点正好对应 \(k=2\) 的负 Pell 解。
C++、Python 和 Java 三个实现都遵循相同的两阶段策略。第一阶段先筛出分界点以内的素数，只保留 \(p\equiv 1\pmod 4\) 的部分，并为每个这样的素数求出 \(x^2\equiv -1\pmod{p^2}\) 的两个提升根。然后递归枚举平方自由素数乘积，用中国剩余定理维护当前模 \(d^2\) 下的全部根，利用“整段周期数加上尾部剩余”的公式统计区间 \([1,n]\) 中的解，并按 Möbius 符号累加。C++ 和 Java 版本把不同起始分支分配到多个线程；Python 版本采用同样的数学流程，但串行执行。
第二阶段先通过筛法去掉所有含有素数 \(q\equiv 3\pmod 4\) 因子的 \(k\)。对每个剩下的非平方 \(k\)，实现利用 \(\sqrt{k}\) 的连分数展开判断负 Pell 方程是否有解，并在有解时得到最小正解。之后按照 Pell 递推不断生成更大的解，直到 \(x>n\)。每当相应的 \(y\) 超过分界点，就用试除法计算 \(\mu(y)\)，并把这个权重加入答案。
此外，基准项 \(d=1\) 会直接贡献 \(n\)。在计算最终目标之前，程序还会先验证 \(C(10)=9\) 和 \(C(1000)=895\) 这两个检查点。
记
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
筛到 \(B\) 的素数以及筛到 \(K\) 的可行 \(k\)，总时间为 \(O(B\log\log B+K\log\log K)\)，空间为 \(O(B+K)\)。小因子阶段的成本与可行平方自由数 \(d\le B\) 的个数及其根集大小成正比；虽然每多一个素因子根数会翻倍，但乘积本身也会迅速增大，因此递归深度在实际参数下很小。大因子阶段遍历所有可行的非平方 \(k\le K\)；每个 \(k\) 的连分数搜索通常较短，而 Pell 解按指数速度增长，所以落在 \(x\le n\) 范围内的项不会很多。整体内存主要由筛数组和当前根缓冲区占据，而不是由 Pell 递推本身占据。
Обозначим
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ свободно от квадратов}\right\}.$$
В задаче 864 требуется найти \(C(123567101113)\). Главная трудность в том, что условие квадратсвободности одновременно запрещает делимость на \(p^2\) для всех простых \(p\). Поэтому реализации не проверяют каждое \(x\) напрямую, а сначала раскладывают индикатор квадратсвободности через Möbius, а затем делят итоговую сумму на малую часть, считаемую по сравнениям, и большую часть, переводимую в отрицательное уравнение Пелля.
Ключевая идея состоит в том, чтобы заменить условие «\(x^2+1\) квадратсвободно» арифметической формулой, после чего остаётся считать решения сравнения \(x^2\equiv -1\) по модулям вида \(d^2\).
Для любого положительного целого \(m\) верно тождество
$$\mathbf{1}_{m\text{ свободно от квадратов}}=\sum_{d^2\mid m}\mu(d).$$
Подставляя \(m=x^2+1\), получаем
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
где
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
Тем самым задача сводится к описанию всех допустимых \(d\) и количества решений по модулю \(d^2\).
Если нечётное простое \(p\) делит \(x^2+1\), то число \(-1\) обязано быть квадратичным вычетом по модулю \(p\). Это происходит ровно при \(p\equiv 1\pmod 4\). Следовательно, простые \(p\equiv 3\pmod 4\) никогда не входят в делители \(d\), дающие вклад. Кроме того, если \(d\) не квадратсвободно, то \(\mu(d)=0\). Значит, достаточно рассматривать только квадратсвободные произведения простых \(p\equiv 1\pmod 4\).
Для такого простого выберем \(r_0\), удовлетворяющее
$$r_0^2\equiv -1\pmod p.$$
Пусть \(r_0^2+1=mp\), и будем искать подъём \(r=r_0+kp\). Тогда
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2},$$
поэтому требуется
$$m+2r_0k\equiv 0\pmod p,$$
то есть
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
Это даёт один корень по модулю \(p^2\), а его отрицание даёт второй. Значит, каждое простое \(p\equiv 1\pmod 4\) создаёт ровно два корня по модулю \(p^2\).
Если
$$d=\prod_{i=1}^{t}p_i$$
квадратсвободно и все \(p_i\equiv 1\pmod 4\), то локальные выборы объединяются независимо по китайской теореме об остатках. Поэтому сравнение \(x^2\equiv -1\pmod{d^2}\) имеет ровно
$$2^{\omega(d)}$$
решений по модулю \(d^2\).
Выберем границу \(B\). Для квадратсвободных \(d\le B\), состоящих только из простых \(p\equiv 1\pmod 4\), обозначим через \(\mathcal{R}_d\) множество корней по модулю \(d^2\). Тогда
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
Так как \(\#\mathcal{R}_d=2^{\omega(d)}\), каждый член \(\mu(d)A_d(n)\) вычисляется явно. Малая часть суммы имеет вид
$$\sum_{\substack{d\le B\\ d\text{ квадратсвободно}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
Реализация перечисляет такие \(d\) рекурсивно. При добавлении нового простого модуль меняется с \(d^2\) на \((dp)^2\), а каждый уже найденный корень разветвляется в два новых через китайскую теорему об остатках. Чередующийся знак полностью совпадает со знаком функции Мёбиуса, то есть сама рекурсия и есть инклюзионно-эксклюзионное разложение.
Когда \(d>B\), прямое перечисление сравнений становится невыгодным. Тогда записываем
$$x^2+1=d^2k$$
и обозначаем \(y=d\). Получаем
$$x^2-ky^2=-1,$$
причём
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Таким образом, каждый большой квадратный делитель соответствует решению отрицательного уравнения Пелля. Параметр \(k\) остаётся сравнительно малым, потому что \(y\) велико.
На \(k\) накладывается сильное арифметическое ограничение. Если бы нечётное простое \(q\equiv 3\pmod 4\) делило \(k\), то оно делило бы и \(x^2+1\), что невозможно. Поэтому допустимые \(k\) не содержат таких простых множителей. Полные квадраты \(k\) тоже исключаются: уравнение должно быть действительно пеллевским.
Для фиксированного допустимого неквадратного \(k\) уравнение
$$x^2-ky^2=-1$$
либо не имеет решений, либо имеет наименьшее положительное решение \((x_0,y_0)\). Разложение \(\sqrt{k}\) в цепную дробь позволяет определить, какой случай реализуется, и в разрешимом случае восстановить это минимальное решение.
Из \((x_0,y_0)\) строим
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
Тогда \(u^2-kv^2=1\), и умножение на \(u+v\sqrt{k}\) сохраняет уравнение \(x^2-ky^2=-1\). Все положительные решения получаются по рекурсии
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
то есть
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
Как только \(y_m>B\), соответствующее решение даёт вклад \(\mu(y_m)\). Если \(y_m\) не квадратсвободно, то \(\mu(y_m)=0\), и вклад исчезает. Тем самым большая часть суммы в точности совпадает с исходной суммой по большим делителям \(d\).
Контрольное значение \(C(10)=9\) видно сразу. Для \(x=1,\dots,10\) величины \(x^2+1\) равны
$$2,5,10,17,26,37,50,65,82,101.$$
Только число \(50\) не квадратсвободно, потому что \(50=2\cdot 5^2\). Следовательно, \(C(10)=9\).
Формула включения-исключения видит то же самое. Термин \(d=1\) даёт \(10\). Единственная ненулевая поправка идёт от \(d=5\), потому что сравнение \(x^2\equiv -1\pmod{25}\) имеет корни \(x\equiv 7,18\pmod{25}\), и в интервале \([1,10]\) присутствует только \(x=7\). Поэтому
$$C(10)=10-\!1=9.$$
Тот же факт в форме Пелля выглядит так:
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1.$$
Реализации на C++, Python и Java следуют одному и тому же двухэтапному плану. Сначала строится решето простых до выбранной границы, после чего сохраняются только простые \(p\equiv 1\pmod 4\). Для каждого такого простого вычисляются два поднятых корня сравнения \(x^2\equiv -1\pmod{p^2}\). Затем малая часть перечисляет квадратсвободные произведения этих простых, объединяет классы остатков с помощью китайской теоремы об остатках, считает число попаданий в \([1,n]\) по формуле «полные периоды плюс хвост» и накапливает сумму с правильным знаком Мёбиуса. Версии на C++ и Java распараллеливают независимые стартовые ветви; версия на Python делает то же самое последовательно.
Во второй фазе решетом отбрасываются все \(k\), кратные простым \(q\equiv 3\pmod 4\). Для каждого оставшегося неквадратного \(k\) реализация использует цепные дроби для \(\sqrt{k}\), чтобы проверить разрешимость отрицательного уравнения Пелля и, если оно разрешимо, найти фундаментальное решение. После этого генерируется вся цепочка решений до тех пор, пока \(x\le n\). Каждый раз, когда соответствующее \(y\) превышает границу, значение \(\mu(y)\) вычисляется пробным делением и добавляется в ответ.
Отдельно базовый член \(d=1\) сразу даёт вклад \(n\). Перед вычислением целевого значения реализации также проверяют себя на контрольных точках \(C(10)=9\) и \(C(1000)=895\).
Пусть
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
Решето простых до \(B\) и допустимых значений \(k\) до \(K\) требует \(O(B\log\log B+K\log\log K)\) времени и \(O(B+K)\) памяти. Малая фаза пропорциональна числу допустимых квадратсвободных произведений \(d\le B\) и размерам их множеств корней; каждый новый простой удваивает число корней, но само произведение быстро растёт, поэтому глубина рекурсии на практике остаётся небольшой. Большая фаза перебирает допустимые неквадратные \(k\le K\); поиск по цепным дробям короток, а решения Пелля растут экспоненциально, поэтому ниже границы \(x\le n\) остаётся лишь немного членов. Основная память расходуется на решета и текущие буферы корней, а не на саму часть с уравнением Пелля.
لنعرّف
$$C(n)=\#\left\{x\in\mathbb{Z}:1\le x\le n,\ x^2+1\text{ squarefree}\right\}.$$
المطلوب في المسألة 864 هو حساب \(C(123567101113)\). الصعوبة أن شرط الخلو من المربعات ليس شرطًا محليًا بسيطًا، بل يعني وجوب استبعاد القسمة على \(p^2\) لكل عدد أولي \(p\). لذلك لا تعتمد التطبيقات على فحص كل \(x\) مباشرة، بل تعيد كتابة الشرط بواسطة دالة Möbius، ثم تقسّم المجموع إلى جزء ذي قواسم مربعة صغيرة يعالج بالتطابقات، وجزء ذي قواسم كبيرة يعالج عبر معادلة Pell سالبة.
الفكرة الأساسية هي تحويل عبارة «\(x^2+1\) خالٍ من المربعات» إلى مؤشر حسابي يمكن جمعه. بعد ذلك يصبح المطلوب عدّ حلول التطابق \(x^2\equiv -1\) بترديدات من الشكل \(d^2\).
لكل عدد صحيح موجب \(m\) لدينا الهوية
$$\mathbf{1}_{m\text{ squarefree}}=\sum_{d^2\mid m}\mu(d).$$
وبتطبيقها على \(m=x^2+1\) نحصل على
$$C(n)=\sum_{x=1}^{n}\sum_{d^2\mid x^2+1}\mu(d)=\sum_{d\ge 1}\mu(d)\,A_d(n),$$
حيث
$$A_d(n)=\#\left\{x:1\le x\le n,\ x^2\equiv -1 \pmod{d^2}\right\}.$$
إذًا صارت المسألة: ما هي قيم \(d\) التي تعطي حلولًا غير صفرية، وكم عدد الحلول لكل منها modulo \(d^2\).
إذا كان عدد أولي فردي \(p\) يقسم \(x^2+1\)، فلا بد أن تكون \(-1\) بقايا تربيعية modulo \(p\)، وهذا يحدث تمامًا عندما \(p\equiv 1\pmod 4\). لذلك لا يمكن لأي عدد أولي \(p\equiv 3\pmod 4\) أن يظهر داخل أي \(d\) ذي مساهمة. كذلك إذا لم يكن \(d\) خاليًا من المربعات فإن \(\mu(d)=0\). لذا لا يبقى إلا القيم \(d\) الخالية من المربعات والمكوّنة من أعداد أولية تحقق \(p\equiv 1\pmod 4\).
لعدد أولي من هذا النوع نختار \(r_0\) بحيث
$$r_0^2\equiv -1\pmod p.$$
نكتب \(r_0^2+1=mp\)، ثم نبحث عن رفع من الشكل \(r=r_0+kp\). عندئذ
$$r^2+1\equiv p\left(m+2r_0k\right)\pmod{p^2},$$
ومن ثم يجب أن يكون
$$m+2r_0k\equiv 0\pmod p,$$
أي
$$k\equiv -m(2r_0)^{-1}\pmod p.$$
وهكذا نحصل على جذر واحد modulo \(p^2\)، ويعطي السالب الجذر الثاني. إذن كل عدد أولي \(p\equiv 1\pmod 4\) يولّد بالضبط جذرين modulo \(p^2\).
إذا كان
$$d=\prod_{i=1}^{t}p_i$$
خالياً من المربعات وجميع عوامله الأولية تحقق \(p_i\equiv 1\pmod 4\)، فإن مبرهنة الباقي الصيني تجمع الاختيارات المحلية بصورة مستقلة. لذلك يكون عدد حلول \(x^2\equiv -1\pmod{d^2}\) مساويًا تمامًا لـ
$$2^{\omega(d)}$$
حيث \(\omega(d)\) هو عدد القواسم الأولية المختلفة لـ \(d\).
نختار حدًا فاصلًا \(B\). ولكل \(d\le B\) خالٍ من المربعات ومكوّن فقط من أعداد أولية تحقق \(p\equiv 1\pmod 4\)، نرمز إلى مجموعة الجذور modulo \(d^2\) بالرمز \(\mathcal{R}_d\). عندئذ
$$A_d(n)=\left\lfloor\frac{n}{d^2}\right\rfloor\#\mathcal{R}_d+\#\left\{r\in\mathcal{R}_d:r\le n\bmod d^2\right\}.$$
وبما أن \(\#\mathcal{R}_d=2^{\omega(d)}\)، فإن كل حد \(\mu(d)A_d(n)\) يصبح صريحًا. ومن ثم يكون الجزء الصغير من المجموع هو
$$\sum_{\substack{d\le B\\ d\text{ squarefree}\\ p\mid d\Rightarrow p\equiv 1\ (\mathrm{mod}\ 4)}}\mu(d)\,A_d(n).$$
التنفيذ يولّد هذه القيم \(d\) توليدًا递归يًا. فكلما أضيف عدد أولي جديد، يتحول المقياس من \(d^2\) إلى \((dp)^2\)، وكل جذر قديم ينقسم إلى جذرين جديدين بواسطة مبرهنة الباقي الصيني. كما أن تبدل الإشارة يطابق تمامًا إشارة دالة Möbius، أي إن هذا التوليد نفسه هو توسيع التضمين والاستبعاد.
عندما يكون \(d>B\)، يصبح عدّ فئات التطابق أقل كفاءة. لذلك نكتب
$$x^2+1=d^2k$$
ونسمي \(y=d\). فنحصل على
$$x^2-ky^2=-1,$$
مع الشروط
$$1\le y,\qquad y>B,\qquad 1\le x\le n,\qquad 1\le k\le \left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
إذًا كل قاسم مربّع كبير يقابل حلًا صحيحًا لمعادلة Pell سالبة، بينما يبقى \(k\) صغيرًا لأن \(y\) كبير.
وهنا يوجد ترشيح حسابي حاد: لو كان عدد أولي فردي \(q\equiv 3\pmod 4\) يقسم \(k\)، لكان يقسم أيضًا \(x^2+1\)، وهذا مستحيل. لذلك لا يجوز للقيم المسموح بها من \(k\) أن تحتوي على مثل هذه العوامل. كما تُستبعد القيم المربعة الكاملة لـ \(k\)، لأن المعادلة يجب أن تكون من نوع Pell الحقيقي غير التافه.
لـ \(k\) ثابت مقبول وغير مربع، فإن المعادلة
$$x^2-ky^2=-1$$
إما بلا حلول، أو لها أصغر حل موجب \((x_0,y_0)\). يحدد التوسع إلى كسور مستمرة لـ \(\sqrt{k}\) أي الحالتين تقع، وإذا كانت المعادلة قابلة للحل فإنه يعطينا هذا الحل الأساسي.
بعد ذلك نكوّن
$$u+v\sqrt{k}=(x_0+y_0\sqrt{k})^2=(x_0^2+ky_0^2)+2x_0y_0\sqrt{k}.$$
عندئذ يتحقق \(u^2-kv^2=1\)، وضرب أي حل في \(u+v\sqrt{k}\) يحافظ على المعادلة \(x^2-ky^2=-1\). ومن هنا تُولد جميع الحلول الموجبة بالعلاقة
$$x_{m+1}+y_{m+1}\sqrt{k}=(x_m+y_m\sqrt{k})(u+v\sqrt{k}),$$
أي
$$x_{m+1}=ux_m+kvy_m,\qquad y_{m+1}=vx_m+uy_m.$$
كلما تحقق \(y_m>B\)، أعطى هذا الحل مساهمة مقدارها \(\mu(y_m)\). وإذا لم يكن \(y_m\) خاليًا من المربعات فإن \(\mu(y_m)=0\) فلا يضيف شيئًا. وبهذا تبقى المساهمة الكبيرة مطابقة تمامًا للمجموع الأصلي على القواسم \(d\).
قيمة الفحص \(C(10)=9\) توضح الفكرة بوضوح. من أجل \(x=1,\dots,10\) نحصل على
$$2,5,10,17,26,37,50,65,82,101.$$
العدد الوحيد غير الخالي من المربعات هو \(50\)، لأن \(50=2\cdot 5^2\). لذا يكون \(C(10)=9\).
وصيغة التضمين والاستبعاد تعطي الشيء نفسه. الحد \(d=1\) يساهم بـ \(10\). والتصحيح الوحيد غير الصفري يأتي من \(d=5\)، لأن التطابق \(x^2\equiv -1\pmod{25}\) له الجذور \(x\equiv 7,18\pmod{25}\)، ومن بينها يقع \(x=7\) فقط داخل المجال \([1,10]\). لذلك
$$C(10)=10-\!1=9.$$
ويظهر المثال نفسه بصيغة Pell على الشكل
$$7^2+1=2\cdot 5^2\qquad\Longleftrightarrow\qquad 7^2-2\cdot 5^2=-1.$$
تتبع تطبيقات C++ وPython وJava الخطة نفسها ذات المرحلتين. في المرحلة الأولى تُجرى غربلة للأعداد الأولية حتى الحد الفاصل، ثم يُحتفظ فقط بالأعداد الأولية التي تحقق \(p\equiv 1\pmod 4\). ولكل منها تُحسب جذور المعادلة \(x^2\equiv -1\pmod{p^2}\) بعد الرفع. ثم يجري توليد جميع المضروبات الخالية من المربعات بصورة تكرارية، وتُدمج مجموعات الجذور باستعمال مبرهنة الباقي الصيني، ويُحسب عدد الحلول داخل \([1,n]\) بصيغة «عدد الدورات الكاملة زائد الذيل»، ثم تُجمع القيم مع إشارة Möbius المناسبة. نسختا C++ وJava توزعان الفروع الابتدائية المستقلة على عدة خيوط، أما نسخة Python فتطبق الرياضيات نفسها بشكل تسلسلي.
في المرحلة الثانية تُغربل القيم الممكنة لـ \(k\) بحذف جميع المضاعفات التي تحتوي على عامل أولي \(q\equiv 3\pmod 4\). ولكل \(k\) متبقٍ وغير مربع، تستخدم الشيفرة الكسور المستمرة لـ \(\sqrt{k}\) لتحديد قابلية حل معادلة Pell السالبة، ثم لاستخراج الحل الأساسي إن وجد. بعد ذلك تُولَّد سلسلة الحلول حتى يصبح \(x>n\). وكلما تجاوز \(y\) الحد الفاصل تُحسب \(\mu(y)\) بالقسمة التجريبية ويضاف هذا الوزن إلى المجموع.
وأخيرًا، يساهم الحد الأساسي \(d=1\) مباشرة بمقدار \(n\). وقبل الانتقال إلى الهدف النهائي، تتحقق التطبيقات أيضًا من النقطتين \(C(10)=9\) و\(C(1000)=895\).
لنضع
$$K=\left\lfloor\frac{n^2+1}{B^2}\right\rfloor.$$
تكلفة غربلة الأعداد الأولية حتى \(B\) وغربلة قيم \(k\) المسموح بها حتى \(K\) هي \(O(B\log\log B+K\log\log K)\) زمنًا و\(O(B+K)\) ذاكرة. أما مرحلة القواسم الصغيرة فتتناسب مع عدد القيم \(d\le B\) الخالية من المربعات ومع أحجام مجموعات جذورها؛ فكل عدد أولي إضافي يضاعف عدد الجذور، لكن حاصل الضرب نفسه ينمو بسرعة، ولذلك يبقى عمق الاستدعاء التكراري صغيرًا عمليًا. ومرحلة القواسم الكبيرة تعالج القيم غير المربعة المسموح بها من \(k\le K\)؛ البحث بالكسور المستمرة قصير عادة، وحلول Pell تنمو أسيًا، لذا فإن عدد الحلول التي تبقى تحت الشرط \(x\le n\) يظل محدودًا. يهيمن على الذاكرة استخدام مصفوفات الغربلة ومخازن الجذور الحالية أكثر من مرحلة Pell نفسها.