Problem 864: Square + 1 = Squarefree
View on Project EulerProject Euler Problem 864 Solution
EulerSolve provides an optimized solution for Project Euler Problem 864, Square + 1 = Squarefree, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: Project Euler 864
- Squarefree integers: Wikipedia — Squarefree integer
- Möbius function: Wikipedia — Möbius function
- Chinese remainder theorem: Wikipedia — Chinese remainder theorem
- Pell's equation: Wikipedia — Pell's equation
- Hensel lifting: Wikipedia — Hensel's lemma
Problem 864 source code
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <vector>
using std::int64_t;
namespace {
using i128 = __int128_t;
constexpr int64_t kTarget = 123567101113LL;
// Balance point between DFS (d <= B) and Pell (d > B). Must be > sqrt(n).
constexpr int64_t kBGlobal = 60000000LL;
constexpr int kMaxRoots = 128; // Enough for B=6e7 (max omega is 6 -> 64 roots).
struct PrimeRoot {
int p;
int64_t p2;
int64_t r1;
int64_t r2;
};
struct Roots {
int size = 0;
int64_t vals[kMaxRoots]{};
};
std::vector<PrimeRoot> g_primes;
std::vector<int> g_small_primes;
int64_t mod_pow(int64_t base, int64_t exp, int64_t mod) {
int64_t res = 1 % mod;
base %= mod;
while (exp > 0) {
if (exp & 1) res = (int64_t)((i128)res * base % mod);
base = (int64_t)((i128)base * base % mod);
exp >>= 1;
}
return res;
}
int64_t egcd(int64_t a, int64_t b, int64_t &x, int64_t &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
int64_t x1 = 0, y1 = 0;
int64_t g = egcd(b, a % b, x1, y1);
x = y1;
y = x1 - (int64_t)((i128)y1 * (a / b));
return g;
}
int64_t mod_inv(int64_t a, int64_t mod) {
int64_t x = 0, y = 0;
int64_t g = egcd(a, mod, x, y);
if (g != 1) return 0;
x %= mod;
if (x < 0) x += mod;
return x;
}
int64_t sqrt_minus_one_mod_p(int64_t p) {
// Find any quadratic non-residue b, then b^((p-1)/4) is a sqrt(-1).
for (int64_t b = 2; b < p; ++b) {
if (mod_pow(b, (p - 1) / 2, p) == p - 1) {
return mod_pow(b, (p - 1) / 4, p);
}
}
return -1;
}
PrimeRoot build_prime_root(int p) {
int64_t r0 = sqrt_minus_one_mod_p(p);
int64_t p2 = (int64_t)p * p;
int64_t m = (r0 * r0 + 1) / p;
int64_t inv = mod_pow((2 * r0) % p, p - 2, p);
int64_t k = (p - (m % p)) % p;
k = (int64_t)((i128)k * inv % p);
int64_t r = r0 + (int64_t)p * k;
return PrimeRoot{p, p2, r, p2 - r};
}
void init_prime_roots(int64_t limit) {
std::vector<bool> is_prime(limit + 1, true);
is_prime[0] = is_prime[1] = false;
for (int64_t i = 2; i * i <= limit; ++i) {
if (is_prime[i]) {
for (int64_t j = i * i; j <= limit; j += i) {
is_prime[(size_t)j] = false;
}
}
}
for (int p = 5; p <= limit; ++p) {
if (is_prime[p] && (p % 4 == 1)) {
g_primes.push_back(build_prime_root(p));
}
}
}
void init_small_primes(int limit) {
std::vector<bool> is_prime(limit + 1, true);
is_prime[0] = is_prime[1] = false;
for (int i = 2; i * i <= limit; ++i) {
if (is_prime[i]) {
for (int j = i * i; j <= limit; j += i) {
is_prime[j] = false;
}
}
}
for (int i = 2; i <= limit; ++i) {
if (is_prime[i]) g_small_primes.push_back(i);
}
}
int64_t count_solutions(int64_t n, int64_t mod, const Roots &roots) {
int64_t q = n / mod;
int64_t s = n - q * mod;
int64_t count = q * roots.size;
for (int i = 0; i < roots.size; ++i) {
if (roots.vals[i] <= s) ++count;
}
return count;
}
void combine_roots(const Roots &roots, const PrimeRoot &pr, int64_t mod_old, Roots &out) {
int64_t mod_new = pr.p2;
int64_t inv = mod_inv(mod_old % mod_new, mod_new);
out.size = roots.size * 2;
if (out.size > kMaxRoots) {
std::cerr << "Root buffer too small.\n";
std::exit(1);
}
int idx = 0;
for (int i = 0; i < roots.size; ++i) {
int64_t ro = roots.vals[i];
int64_t ro_mod = ro % mod_new;
int64_t r_list[2] = {pr.r1, pr.r2};
for (int j = 0; j < 2; ++j) {
int64_t diff = r_list[j] - ro_mod;
if (diff < 0) diff += mod_new;
int64_t t = (int64_t)((i128)diff * inv % mod_new);
int64_t new_r = ro + (int64_t)((i128)mod_old * t);
out.vals[idx++] = new_r;
}
}
}
void dfs_part1(int start_idx,
int64_t d,
int64_t d2,
const Roots &roots,
int sign,
int64_t n,
int64_t B,
int64_t &sum) {
for (int i = start_idx; i < (int)g_primes.size(); ++i) {
const PrimeRoot &pr = g_primes[i];
if (pr.p > B) break;
if (d > B / pr.p) break;
int64_t next_d = d * pr.p;
int64_t next_d2 = next_d * next_d;
Roots next_roots;
combine_roots(roots, pr, d2, next_roots);
int64_t term = count_solutions(n, next_d2, next_roots);
int next_sign = -sign;
sum += (int64_t)next_sign * term;
dfs_part1(i + 1, next_d, next_d2, next_roots, next_sign, n, B, sum);
}
}
int64_t count_part1(int64_t n, int64_t B) {
if (B < 5) return 0;
int threads = (int)std::thread::hardware_concurrency();
if (threads <= 0) threads = 4;
std::vector<int64_t> sums(threads, 0);
std::vector<std::thread> workers;
for (int t = 0; t < threads; ++t) {
workers.emplace_back([&, t]() {
int64_t local = 0;
for (int idx = t; idx < (int)g_primes.size(); idx += threads) {
const PrimeRoot &pr = g_primes[idx];
if (pr.p > B) break;
int64_t d = pr.p;
int64_t d2 = d * d;
Roots roots;
roots.size = 2;
roots.vals[0] = pr.r1;
roots.vals[1] = pr.r2;
int64_t term = count_solutions(n, d2, roots);
local -= term; // mu(p) = -1
dfs_part1(idx + 1, d, d2, roots, -1, n, B, local);
}
sums[t] = local;
});
}
for (auto &th : workers) th.join();
int64_t total = 0;
for (int64_t v : sums) total += v;
return total;
}
bool solve_negative_pell(int64_t k, int64_t n, int64_t &x, int64_t &y) {
int64_t a0 = (int64_t)std::sqrt((long double)k);
while ((a0 + 1) * (a0 + 1) <= k) ++a0;
while (a0 * a0 > k) --a0;
if (a0 * a0 == k) return false;
int64_t m = 0, d = 1, a = a0;
i128 p_prev1 = 1;
i128 q_prev1 = 0;
i128 p_curr = a0, q_curr = 1;
if (p_curr > n || q_curr > n) return false;
for (int step = 1; ; ++step) {
m = d * a - m;
d = (k - m * m) / d;
a = (a0 + m) / d;
i128 p_next = (i128)a * p_curr + p_prev1;
i128 q_next = (i128)a * q_curr + q_prev1;
p_prev1 = p_curr;
q_prev1 = q_curr;
p_curr = p_next;
q_curr = q_next;
if (d == 1 && a == 2 * a0) {
if (step % 2 == 1 && p_prev1 <= n && q_prev1 <= n) {
x = (int64_t)p_prev1;
y = (int64_t)q_prev1;
return true;
}
return false;
}
if (p_curr > n || q_curr > n) return false;
}
}
int mobius_squarefree(int64_t y) {
int mu = 1;
int64_t tmp = y;
for (int p : g_small_primes) {
int64_t pp = (int64_t)p * p;
if (pp > tmp) break;
if (tmp % p == 0) {
tmp /= p;
if (tmp % p == 0) return 0;
mu = -mu;
}
}
if (tmp > 1) mu = -mu;
return mu;
}
std::vector<char> build_valid_k(int64_t k_max) {
std::vector<char> valid(k_max + 1, 1);
if (k_max >= 0) valid[0] = 0;
std::vector<bool> is_prime(k_max + 1, true);
if (k_max >= 0) is_prime[0] = false;
if (k_max >= 1) is_prime[1] = false;
for (int64_t i = 2; i * i <= k_max; ++i) {
if (is_prime[i]) {
for (int64_t j = i * i; j <= k_max; j += i) is_prime[(size_t)j] = false;
}
}
for (int64_t p = 2; p <= k_max; ++p) {
if (!is_prime[p]) continue;
if (p % 4 == 3) {
for (int64_t j = p; j <= k_max; j += p) valid[(size_t)j] = 0;
}
}
return valid;
}
int64_t count_part2(int64_t n, int64_t B) {
int64_t k_max = (int64_t)(((i128)n * n + 1) / ((i128)B * B));
if (k_max <= 0) return 0;
std::vector<char> valid = build_valid_k(k_max);
int threads = (int)std::thread::hardware_concurrency();
if (threads <= 0) threads = 4;
std::vector<int64_t> sums(threads, 0);
std::vector<std::thread> workers;
for (int t = 0; t < threads; ++t) {
workers.emplace_back([&, t]() {
int64_t local = 0;
for (int64_t k = 1 + t; k <= k_max; k += threads) {
if (!valid[(size_t)k]) continue;
int64_t s = (int64_t)std::sqrt((long double)k);
if (s * s == k) continue;
int64_t x0 = 0, y0 = 0;
if (!solve_negative_pell(k, n, x0, y0)) continue;
i128 X_mul = (i128)x0 * x0 + (i128)k * y0 * y0;
i128 Y_mul = 2 * (i128)x0 * y0;
i128 Y_mul_k = Y_mul * k;
int64_t x = x0;
int64_t y = y0;
while (x <= n) {
if (y > B) {
int mu = mobius_squarefree(y);
if (mu != 0) local += mu;
}
if (X_mul > n || Y_mul_k > n) break;
i128 x_next = (i128)x * X_mul + (i128)y * Y_mul_k;
i128 y_next = (i128)x * Y_mul + (i128)y * X_mul;
if (x_next > n) break;
x = (int64_t)x_next;
y = (int64_t)y_next;
}
}
sums[t] = local;
});
}
for (auto &th : workers) th.join();
int64_t total = 0;
for (int64_t v : sums) total += v;
return total;
}
int64_t compute_C(int64_t n) {
int64_t B = std::min(kBGlobal, n);
int64_t total = n; // d=1
total += count_part1(n, B);
total += count_part2(n, B);
return total;
}
int64_t brute_force_C(int64_t n) {
int64_t cnt = 0;
for (int64_t x = 1; x <= n; ++x) {
int64_t v = x * x + 1;
bool squarefree = true;
for (int p : g_small_primes) {
int64_t pp = (int64_t)p * p;
if (pp > v) break;
if (v % pp == 0) {
squarefree = false;
break;
}
}
if (squarefree) ++cnt;
}
return cnt;
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
std::cout << "Building prime roots...\n";
init_prime_roots(kBGlobal);
init_small_primes(400000);
std::cout << "Validation checks...\n";
int64_t c10 = compute_C(10);
int64_t c10_bf = brute_force_C(10);
std::cout << "C(10) = " << c10 << " (expected 9, brute " << c10_bf << ")\n";
if (c10 != 9 || c10_bf != 9) {
std::cerr << "Validation failed for C(10).\n";
return 1;
}
int64_t c1000 = compute_C(1000);
int64_t c1000_bf = brute_force_C(1000);
std::cout << "C(1000) = " << c1000 << " (expected 895, brute " << c1000_bf << ")\n";
if (c1000 != 895 || c1000_bf != 895) {
std::cerr << "Validation failed for C(1000).\n";
return 1;
}
std::cout << "Computing target...\n";
int64_t result = compute_C(kTarget);
std::cout << "C(" << kTarget << ") = " << result << "\n";
return 0;
}
Python
import math
def solve():
TARGET = 123567101113; B = 60000000
def mod_pow(base, exp, mod):
r = 1; base %= mod
while exp > 0:
if exp & 1: r = r * base % mod
base = base * base % mod; exp >>= 1
return r
def egcd(a, b):
if b == 0: return a, 1, 0
g, x1, y1 = egcd(b, a%b)
return g, y1, x1 - (a//b)*y1
def mod_inv(a, m):
g, x, _ = egcd(a % m, m)
return x % m if g == 1 else 0
def sqrt_m1_mod_p(p):
for b in range(2, p):
if mod_pow(b, (p-1)//2, p) == p-1:
return mod_pow(b, (p-1)//4, p)
return -1
# Build prime roots for p ≡ 1 (mod 4) up to B
sieve = bytearray(b'\x01')*(B+1); sieve[0] = sieve[1] = 0
for i in range(2, int(B**0.5)+1):
if sieve[i]:
for j in range(i*i, B+1, i): sieve[j] = 0
prime_roots = []
for p in range(5, B+1):
if sieve[p] and p % 4 == 1:
r0 = sqrt_m1_mod_p(p); p2 = p*p
m = (r0*r0+1)//p; inv = mod_pow(2*r0%p, p-2, p)
k = (p - m%p) % p; k = k*inv%p
r = r0 + p*k
prime_roots.append((p, p2, r, p2-r))
# Small primes for Mobius
slp = []; sp_sieve = bytearray(b'\x01')*400001; sp_sieve[0] = sp_sieve[1] = 0
for i in range(2, 400001):
if sp_sieve[i]:
slp.append(i)
for j in range(i*i, 400001, i): sp_sieve[j] = 0
def count_sols(n, mod, roots):
q = n // mod; s = n - q*mod
return q*len(roots) + sum(1 for r in roots if r <= s)
def combine(roots, pr, mod_old):
_, p2, r1, r2 = pr; inv = mod_inv(mod_old % p2, p2)
out = []
for ro in roots:
rm = ro % p2
for rn in (r1, r2):
diff = (rn - rm) % p2
t = diff * inv % p2
out.append(ro + mod_old * t)
return out
def dfs1(si, d, d2, roots, sign, n, total):
for i in range(si, len(prime_roots)):
p, p2, r1, r2 = prime_roots[i]
if p > B or d > B // p: break
nd = d*p; nd2 = nd*nd
nr = combine(roots, prime_roots[i], d2)
term = count_sols(n, nd2, nr)
ns = -sign; total[0] += ns * term
dfs1(i+1, nd, nd2, nr, ns, n, total)
total1 = [0]
for i in range(len(prime_roots)):
p, p2, r1, r2 = prime_roots[i]
if p > B: break
roots = [r1, r2]; term = count_sols(TARGET, p2, roots)
total1[0] -= term
dfs1(i+1, p, p2, roots, -1, TARGET, total1)
# Part 2: Pell equation
def solve_neg_pell(k, n):
a0 = int(k**0.5)
while (a0+1)**2 <= k: a0 += 1
while a0*a0 > k: a0 -= 1
if a0*a0 == k: return None
m = 0; d = 1; a = a0
pp1, qp1 = 1, 0; pc, qc = a0, 1
if pc > n: return None
for step in range(1, 10000):
m = d*a - m; d = (k - m*m)//d; a = (a0 + m)//d
pn = a*pc + pp1; qn = a*qc + qp1
pp1, qp1 = pc, qc; pc, qc = pn, qn
if d == 1 and a == 2*a0:
if step%2 == 1 and pp1 <= n and qp1 <= n: return pp1, qp1
return None
if pc > n: return None
return None
def mobius_sf(y):
mu = 1; tmp = y
for p in slp:
if p*p > tmp: break
if tmp%p == 0:
tmp //= p
if tmp%p == 0: return 0
mu = -mu
if tmp > 1: mu = -mu
return mu
k_max = (TARGET*TARGET + 1) // (B*B)
valid = bytearray(b'\x01')*(k_max+1); valid[0] = 0
ksv = bytearray(b'\x01')*(k_max+1); ksv[0] = ksv[1] = 0
for i in range(2, k_max+1):
if ksv[i]:
if i%4 == 3:
for j in range(i, k_max+1, i): valid[j] = 0
for j in range(i*i, k_max+1, i): ksv[j] = 0
total2 = 0
for k in range(1, k_max+1):
if not valid[k]: continue
s = int(k**0.5)
if s*s == k: continue
res = solve_neg_pell(k, TARGET)
if res is None: continue
x0, y0 = res
X_mul = x0*x0 + k*y0*y0; Y_mul = 2*x0*y0
x, y = x0, y0
while x <= TARGET:
if y > B:
mu = mobius_sf(y)
if mu != 0: total2 += mu
xn = x*X_mul + y*Y_mul*k; yn = x*Y_mul + y*X_mul
if xn > TARGET: break
x, y = xn, yn
return str(TARGET + total1[0] + total2)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.math.BigInteger;
public class Euler864 {
static final long kTarget = 123567101113L;
static final long kBGlobal = 60000000L;
static final int kMaxRoots = 128;
static class PrimeRoot {
int p;
long p2;
long r1;
long r2;
PrimeRoot(int p, long p2, long r1, long r2) {
this.p = p;
this.p2 = p2;
this.r1 = r1;
this.r2 = r2;
}
}
static class Roots {
int size = 0;
long[] vals = new long[kMaxRoots];
}
static ArrayList<PrimeRoot> gPrimes = new ArrayList<>();
static ArrayList<Integer> gSmallPrimes = new ArrayList<>();
static long modPow(long base, long exp, long mod) {
long res = 1 % mod;
base %= mod;
while (exp > 0) {
if ((exp & 1) == 1)
res = BigInteger.valueOf(res).multiply(BigInteger.valueOf(base)).mod(BigInteger.valueOf(mod))
.longValue();
base = BigInteger.valueOf(base).multiply(BigInteger.valueOf(base)).mod(BigInteger.valueOf(mod)).longValue();
exp >>= 1;
}
return res;
}
static long[] egcd(long a, long b) {
if (b == 0)
return new long[] { a, 1, 0 };
long[] res = egcd(b, a % b);
long g = res[0];
long x1 = res[1];
long y1 = res[2];
long x = y1;
long y = x1 - BigInteger.valueOf(y1).multiply(BigInteger.valueOf(a / b)).longValue();
return new long[] { g, x, y };
}
static long modInv(long a, long mod) {
long[] res = egcd(a, mod);
if (res[0] != 1)
return 0;
long x = res[1] % mod;
if (x < 0)
x += mod;
return x;
}
static long sqrtMinusOneModP(long p) {
for (long b = 2; b < p; ++b) {
if (modPow(b, (p - 1) / 2, p) == p - 1) {
return modPow(b, (p - 1) / 4, p);
}
}
return -1;
}
static PrimeRoot buildPrimeRoot(int p) {
long r0 = sqrtMinusOneModP(p);
long p2 = (long) p * p;
long m = (r0 * r0 + 1) / p;
long inv = modPow((2 * r0) % p, p - 2, p);
long k = (p - (m % p)) % p;
k = BigInteger.valueOf(k).multiply(BigInteger.valueOf(inv)).mod(BigInteger.valueOf(p)).longValue();
long r = r0 + (long) p * k;
return new PrimeRoot(p, p2, r, p2 - r);
}
static void initPrimeRoots(long limit) {
boolean[] isPrime = new boolean[(int) limit + 1];
java.util.Arrays.fill(isPrime, true);
isPrime[0] = isPrime[1] = false;
for (long i = 2; i * i <= limit; ++i) {
if (isPrime[(int) i]) {
for (long j = i * i; j <= limit; j += i) {
isPrime[(int) j] = false;
}
}
}
for (int p = 5; p <= limit; ++p) {
if (isPrime[p] && (p % 4 == 1)) {
gPrimes.add(buildPrimeRoot(p));
}
}
}
static void initSmallPrimes(int limit) {
boolean[] isPrime = new boolean[limit + 1];
java.util.Arrays.fill(isPrime, true);
isPrime[0] = isPrime[1] = false;
for (int i = 2; i * i <= limit; ++i) {
if (isPrime[i]) {
for (int j = i * i; j <= limit; j += i) {
isPrime[j] = false;
}
}
}
for (int i = 2; i <= limit; ++i) {
if (isPrime[i])
gSmallPrimes.add(i);
}
}
static long countSolutions(long n, long mod, Roots roots) {
long q = n / mod;
long s = n - q * mod;
long count = q * roots.size;
for (int i = 0; i < roots.size; ++i) {
if (roots.vals[i] <= s)
count++;
}
return count;
}
static void combineRoots(Roots roots, PrimeRoot pr, long modOld, Roots out) {
long modNew = pr.p2;
long inv = modInv(modOld % modNew, modNew);
out.size = roots.size * 2;
int idx = 0;
for (int i = 0; i < roots.size; ++i) {
long ro = roots.vals[i];
long roMod = ro % modNew;
long[] rList = { pr.r1, pr.r2 };
for (int j = 0; j < 2; ++j) {
long diff = rList[j] - roMod;
if (diff < 0)
diff += modNew;
long t = BigInteger.valueOf(diff).multiply(BigInteger.valueOf(inv)).mod(BigInteger.valueOf(modNew))
.longValue();
long newR = ro + BigInteger.valueOf(modOld).multiply(BigInteger.valueOf(t)).longValue();
out.vals[idx++] = newR;
}
}
}
static void dfsPart1(int startIdx, long d, long d2, Roots roots, int sign, long n, long B, long[] sum) {
for (int i = startIdx; i < gPrimes.size(); ++i) {
PrimeRoot pr = gPrimes.get(i);
if (pr.p > B)
break;
if (d > B / pr.p)
break;
long nextD = d * pr.p;
long nextD2 = nextD * nextD;
Roots nextRoots = new Roots();
combineRoots(roots, pr, d2, nextRoots);
long term = countSolutions(n, nextD2, nextRoots);
int nextSign = -sign;
sum[0] += (long) nextSign * term;
dfsPart1(i + 1, nextD, nextD2, nextRoots, nextSign, n, B, sum);
}
}
static long countPart1(long n, long B) {
if (B < 5)
return 0;
int threads = Runtime.getRuntime().availableProcessors();
if (threads <= 0)
threads = 4;
long[] sums = new long[threads];
Thread[] workers = new Thread[threads];
for (int t = 0; t < threads; ++t) {
final int tId = t;
final int thr = threads;
workers[t] = new Thread(() -> {
long local = 0;
for (int idx = tId; idx < gPrimes.size(); idx += thr) {
PrimeRoot pr = gPrimes.get(idx);
if (pr.p > B)
break;
long d = pr.p;
long d2 = d * d;
Roots roots = new Roots();
roots.size = 2;
roots.vals[0] = pr.r1;
roots.vals[1] = pr.r2;
long term = countSolutions(n, d2, roots);
local -= term;
long[] localArr = { local };
dfsPart1(idx + 1, d, d2, roots, -1, n, B, localArr);
local = localArr[0];
}
sums[tId] = local;
});
workers[t].start();
}
for (Thread th : workers) {
try {
th.join();
} catch (InterruptedException e) {
}
}
long total = 0;
for (long v : sums)
total += v;
return total;
}
static boolean solveNegativePell(long k, long n, long[] xy) {
long a0 = (long) Math.sqrt(k);
while ((a0 + 1) * (a0 + 1) <= k)
a0++;
while (a0 * a0 > k)
a0--;
if (a0 * a0 == k)
return false;
long m = 0, d = 1, a = a0;
BigInteger pPrev1 = BigInteger.ONE;
BigInteger qPrev1 = BigInteger.ZERO;
BigInteger pCurr = BigInteger.valueOf(a0);
BigInteger qCurr = BigInteger.ONE;
BigInteger nBig = BigInteger.valueOf(n);
if (pCurr.compareTo(nBig) > 0 || qCurr.compareTo(nBig) > 0)
return false;
for (int step = 1;; ++step) {
m = d * a - m;
d = (k - m * m) / d;
a = (a0 + m) / d;
BigInteger aBig = BigInteger.valueOf(a);
BigInteger pNext = aBig.multiply(pCurr).add(pPrev1);
BigInteger qNext = aBig.multiply(qCurr).add(qPrev1);
pPrev1 = pCurr;
qPrev1 = qCurr;
pCurr = pNext;
qCurr = qNext;
if (d == 1 && a == 2 * a0) {
if (step % 2 == 1 && pPrev1.compareTo(nBig) <= 0 && qPrev1.compareTo(nBig) <= 0) {
xy[0] = pPrev1.longValue();
xy[1] = qPrev1.longValue();
return true;
}
return false;
}
if (pCurr.compareTo(nBig) > 0 || qCurr.compareTo(nBig) > 0)
return false;
}
}
static int mobiusSquarefree(long y) {
int mu = 1;
long tmp = y;
for (int p : gSmallPrimes) {
long pp = (long) p * p;
if (pp > tmp)
break;
if (tmp % p == 0) {
tmp /= p;
if (tmp % p == 0)
return 0;
mu = -mu;
}
}
if (tmp > 1)
mu = -mu;
return mu;
}
static byte[] buildValidK(long kMax) {
byte[] valid = new byte[(int) kMax + 1];
java.util.Arrays.fill(valid, (byte) 1);
if (kMax >= 0)
valid[0] = 0;
boolean[] isPrime = new boolean[(int) kMax + 1];
java.util.Arrays.fill(isPrime, true);
if (kMax >= 0)
isPrime[0] = false;
if (kMax >= 1)
isPrime[1] = false;
for (long i = 2; i * i <= kMax; ++i) {
if (isPrime[(int) i]) {
for (long j = i * i; j <= kMax; j += i)
isPrime[(int) j] = false;
}
}
for (long p = 2; p <= kMax; ++p) {
if (!isPrime[(int) p])
continue;
if (p % 4 == 3) {
for (long j = p; j <= kMax; j += p)
valid[(int) j] = 0;
}
}
return valid;
}
static long countPart2(long n, long B) {
BigInteger bn = BigInteger.valueOf(n);
BigInteger bB = BigInteger.valueOf(B);
BigInteger n2 = bn.multiply(bn).add(BigInteger.ONE);
BigInteger B2 = bB.multiply(bB);
long kMax = n2.divide(B2).longValue();
if (kMax <= 0)
return 0;
byte[] valid = buildValidK(kMax);
int threads = Runtime.getRuntime().availableProcessors();
if (threads <= 0)
threads = 4;
long[] sums = new long[threads];
Thread[] workers = new Thread[threads];
for (int t = 0; t < threads; ++t) {
final int tId = t;
final int thr = threads;
workers[t] = new Thread(() -> {
long local = 0;
for (long k = 1 + tId; k <= kMax; k += thr) {
if (valid[(int) k] == 0)
continue;
long s = (long) Math.sqrt(k);
if (s * s == k)
continue;
long[] xy = new long[2];
if (!solveNegativePell(k, n, xy))
continue;
long x0 = xy[0];
long y0 = xy[1];
BigInteger XMul = BigInteger.valueOf(x0).multiply(BigInteger.valueOf(x0))
.add(BigInteger.valueOf(k).multiply(BigInteger.valueOf(y0))
.multiply(BigInteger.valueOf(y0)));
BigInteger YMul = BigInteger.valueOf(2).multiply(BigInteger.valueOf(x0))
.multiply(BigInteger.valueOf(y0));
BigInteger YMulK = YMul.multiply(BigInteger.valueOf(k));
BigInteger nBig = BigInteger.valueOf(n);
long x = x0;
long y = y0;
while (x <= n) {
if (y > B) {
int mu = mobiusSquarefree(y);
if (mu != 0)
local += mu;
}
if (XMul.compareTo(nBig) > 0 || YMulK.compareTo(nBig) > 0)
break;
BigInteger xNext = BigInteger.valueOf(x).multiply(XMul)
.add(BigInteger.valueOf(y).multiply(YMulK));
BigInteger yNext = BigInteger.valueOf(x).multiply(YMul)
.add(BigInteger.valueOf(y).multiply(XMul));
if (xNext.compareTo(nBig) > 0)
break;
x = xNext.longValue();
y = yNext.longValue();
}
}
sums[tId] = local;
});
workers[t].start();
}
for (Thread th : workers) {
try {
th.join();
} catch (InterruptedException e) {
}
}
long total = 0;
for (long v : sums)
total += v;
return total;
}
static long computeC(long n) {
long B = Math.min(kBGlobal, n);
long total = n;
total += countPart1(n, B);
total += countPart2(n, B);
return total;
}
public static String solve() {
initPrimeRoots(Math.min(kBGlobal, kTarget));
initSmallPrimes(400000);
return Long.toString(computeC(kTarget));
}
public static void main(String[] args) {
System.out.println(solve());
}
}