Detailed mathematical approach

Problem Summary

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.

Mathematical Approach

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.

Step 1: Möbius Inversion for Squarefreeness

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.

Step 2: Roots Modulo \(p^2\) and the Chinese Remainder Theorem

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\).

Step 3: Counting the Small-Divisor Contribution

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.

Step 4: Rewriting the Large-Divisor Part as a Negative Pell Equation

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.

Step 5: Generating All Large Solutions from the Fundamental One

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\).

Step 6: Worked Example

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\).

How the Code Works

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.

Complexity Analysis

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.

Footnotes and References

1. Problem page: Project Euler 864
2. Squarefree integers: Wikipedia — Squarefree integer
3. Möbius function: Wikipedia — Möbius function
4. Chinese remainder theorem: Wikipedia — Chinese remainder theorem
5. Pell's equation: Wikipedia — Pell's equation
6. Hensel lifting: Wikipedia — Hensel's lemma