Problem 827: Pythagorean Triple Occurrence
View on Project EulerProject Euler Problem 827 Solution
EulerSolve provides an optimized solution for Project Euler Problem 827, Pythagorean Triple Occurrence, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For a positive integer \(N\), let \(T(N)\) be the number of distinct positive integer Pythagorean triples \((a,b,c)\) with \(a^2+b^2=c^2\) in which \(N\) appears as one of the three side lengths, counting \((a,b,c)\) and \((b,a,c)\) as the same triple. The task is to find, for each target count \(10^k\) with \(1\le k\le 18\), the smallest \(N\) satisfying \(T(N)=10^k\), and then sum those minima modulo \(409120391\). A direct search over triples or over candidate integers is infeasible. The solution instead derives a closed multiplicative formula for \(T(N)\), then inverts that formula by a constrained minimization over prime exponents. Mathematical Approach Write the prime factorization of \(N\) as $$N=2^{e_2}\prod_{p\equiv 1 \pmod{4}} p^{\alpha_p}\prod_{q\equiv 3 \pmod{4}} q^{\beta_q}.$$ Define the two odd products $$A=\prod_{p\equiv 1 \pmod{4}} (2\alpha_p+1),\qquad R=\prod_{q\equiv 3 \pmod{4}} (2\beta_q+1).$$ It is also convenient to encode the power of \(2\) by $$u=\begin{cases} 1,& e_2=0,\\ 2e_2-1,& e_2\ge 1. \end{cases}$$ Step 1: Count triples with hypotenuse \(N\) Distinct right triangles with hypotenuse \(N\) correspond to positive representations of \(N^2\) as a sum of two squares. A standard consequence of the sum-of-two-squares theorem is $$r_2(N^2)=4A,$$ where \(r_2(m)\) counts ordered integer pairs \((x,y)\) satisfying \(x^2+y^2=m\)....
Detailed mathematical approach
Problem Summary
For a positive integer \(N\), let \(T(N)\) be the number of distinct positive integer Pythagorean triples \((a,b,c)\) with \(a^2+b^2=c^2\) in which \(N\) appears as one of the three side lengths, counting \((a,b,c)\) and \((b,a,c)\) as the same triple. The task is to find, for each target count \(10^k\) with \(1\le k\le 18\), the smallest \(N\) satisfying \(T(N)=10^k\), and then sum those minima modulo \(409120391\).
A direct search over triples or over candidate integers is infeasible. The solution instead derives a closed multiplicative formula for \(T(N)\), then inverts that formula by a constrained minimization over prime exponents.
Mathematical Approach
Write the prime factorization of \(N\) as
$$N=2^{e_2}\prod_{p\equiv 1 \pmod{4}} p^{\alpha_p}\prod_{q\equiv 3 \pmod{4}} q^{\beta_q}.$$
Define the two odd products
$$A=\prod_{p\equiv 1 \pmod{4}} (2\alpha_p+1),\qquad R=\prod_{q\equiv 3 \pmod{4}} (2\beta_q+1).$$
It is also convenient to encode the power of \(2\) by
$$u=\begin{cases} 1,& e_2=0,\\ 2e_2-1,& e_2\ge 1. \end{cases}$$
Step 1: Count triples with hypotenuse \(N\)
Distinct right triangles with hypotenuse \(N\) correspond to positive representations of \(N^2\) as a sum of two squares. A standard consequence of the sum-of-two-squares theorem is
$$r_2(N^2)=4A,$$
where \(r_2(m)\) counts ordered integer pairs \((x,y)\) satisfying \(x^2+y^2=m\).
Among these representations, four are degenerate: \((\pm N,0)\) and \((0,\pm N)\). Every genuine triangle is counted eight times because of sign choices and swapping the two legs. Therefore the number of distinct triangles having hypotenuse \(N\) is
$$H(N)=\frac{r_2(N^2)-4}{8}=\frac{A-1}{2}.$$
Only primes congruent to \(1 \pmod{4}\) affect this count. Primes congruent to \(3 \pmod{4}\) and the factor \(2\) only appear as square contributions in \(N^2\), so they do not create additional representations.
Step 2: Count triples with leg \(N\)
If \(N\) is odd, every factorization \(N^2=st\) with \(s<t\) gives one triangle
$$\left(N,\frac{t-s}{2},\frac{t+s}{2}\right),$$
so the number of such triangles is the number of unordered divisor pairs of \(N^2\):
$$L(N)=\frac{d(N^2)-1}{2}\qquad (e_2=0).$$
If \(N\) is even, write \(N=2m\). Then every factorization \(m^2=st\) with \(s<t\) gives
$$\left(N,t-s,t+s\right),$$
hence
$$L(N)=\frac{d(m^2)-1}{2}=\frac{d\!\left((N/2)^2\right)-1}{2}\qquad (e_2\ge 1).$$
Using the factorization of \(N\), both cases collapse to one formula:
$$L(N)=\frac{uAR-1}{2}.$$
Step 3: Combine the two occurrence counts
A number cannot be both a leg and the hypotenuse in the same right triangle, so the total occurrence count is simply
$$T(N)=H(N)+L(N)=\frac{A(uR+1)}{2}-1.$$
This is the key identity used by the solver. If we want exactly \(n\) occurrences, then after setting
$$s=n+1,$$
we must satisfy
$$s=\frac{A(uR+1)}{2}.$$
Since \(A\) is odd, every solution begins by choosing an odd divisor \(a\mid s\), forcing
$$uR=\frac{2s}{a}-1.$$
The right-hand side is again odd, so we split it once more as
$$u\cdot r=\frac{2s}{a}-1,$$
where \(u\) is either \(1\) or of the form \(2e_2-1\), and \(r\) must equal \(\prod_{q\equiv 3 \pmod{4}} (2\beta_q+1)\).
Step 4: Turn odd factors into prime exponents
The quantity \(A\) must be represented as a product of odd factors
$$A=\prod_i f_i,\qquad f_i=2\alpha_i+1,$$
and similarly
$$R=\prod_j g_j,\qquad g_j=2\beta_j+1.$$
Once an odd factor is chosen, the exponent is fixed: \(\alpha=(f-1)/2\) or \(\beta=(g-1)/2\). So the inverse problem becomes: partition the required odd product into factors of the form \(2e+1\), then place those exponents on allowed primes.
To minimize \(N\), larger exponents must be assigned to smaller primes. Indeed, if \(p<q\) and \(x>y\), then
$$p^x q^y < p^y q^x.$$
Therefore the optimal partition in each residue class is searched in nonincreasing odd factors and assigned to the increasing primes of that class.
Step 5: Compare candidates by logarithms
The actual minima can be enormous, so the implementations compare candidates using high-precision base-2 logarithms:
$$\log_2 N=e_2+\sum_i \alpha_i\log_2 p_i+\sum_j \beta_j\log_2 q_j.$$
Because the logarithm is strictly increasing, the candidate with smaller logarithmic value is exactly the smaller integer. Once the best exponent pattern is known, the final value is reconstructed modulo \(409120391\) by modular exponentiation.
Worked Example: \(N=15\)
We have \(15=3^1\cdot 5^1\), so
$$A=3,\qquad R=3,\qquad u=1.$$
The hypotenuse contribution is
$$H(15)=\frac{3-1}{2}=1,$$
coming from \((9,12,15)\).
The leg contribution is
$$L(15)=\frac{1\cdot 3\cdot 3-1}{2}=4,$$
coming from \((8,15,17)\), \((15,20,25)\), \((15,36,39)\), and \((15,112,113)\).
Therefore
$$T(15)=1+4=5.$$
So 15 is a valid witness for target count \(5\), and the implementations verify that it is the smallest such integer.
How the Code Works
The C++, Python, and Java implementations first generate enough small primes in the residue classes \(1 \pmod{4}\) and \(3 \pmod{4}\). They also use fast 64-bit primality testing and factorization so that every odd divisor needed by the search can be generated and cached efficiently.
For a fixed odd product in one residue class, the implementation runs a memoized depth-first search over its odd divisors. Choosing one divisor corresponds to choosing one factor \(2e+1\), hence one prime exponent. The recursion enforces nonincreasing chosen factors, which avoids duplicate partitions and automatically places larger exponents on smaller primes.
For a target count \(n\), the implementation sets \(s=n+1\), enumerates odd divisors \(a\mid s\), computes the forced odd remainder \(2s/a-1\), and then tries every odd divisor of that remainder as the possible \(u\) attached to the power of \(2\). The remaining factor is solved in the \(3 \pmod{4}\) prime class, while \(a\) itself is solved in the \(1 \pmod{4}\) prime class. The best combined logarithmic score gives the smallest integer with exactly \(n\) occurrences.
Finally, the implementation evaluates this inverse solver for \(n=10,10^2,\dots,10^{18}\), reconstructs each minimum modulo \(409120391\), and sums those values modulo the same modulus.
Complexity Analysis
The dominant cost comes from factoring the odd integers derived from \(n+1\) and enumerating their odd divisors. If \(\tau(m)\) denotes the divisor function, then each recursive optimization depends on the divisor structure of the relevant odd product rather than on the size of the final minimum itself. Memoization collapses repeated subproblems determined by the remaining product, the next prime slot, and the previously chosen odd factor.
For the concrete targets \(10^1,\dots,10^{18}\), the search trees stay manageable because every recursive step removes at least one odd factor \(2e+1\). In practice, runtime is dominated by Pollard-rho factorization plus divisor generation, while memory is dominated by cached divisor lists and memoized optimal subproblems.
Footnotes and References
- Problem page: Project Euler 827
- Pythagorean triples: Wikipedia — Pythagorean triple
- Sum of two squares theorem: Wikipedia — Fermat's theorem on sums of two squares
- Divisor function: Wikipedia — Divisor function
- Pollard's rho algorithm: Wikipedia — Pollard's rho algorithm
Problem 827 source code
C++
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <map>
#include <numeric>
#include <random>
#include <tuple>
#include <unordered_map>
#include <vector>
#include <boost/multiprecision/cpp_dec_float.hpp>
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using Real = boost::multiprecision::cpp_dec_float_100;
static constexpr u64 kMod = 409120391ULL;
struct Rep {
bool valid = false;
Real log2v = 0;
std::vector<u64> exps;
};
struct BestD {
bool valid = false;
Real log2v = 0;
u64 e2 = 0;
Rep rep3;
};
static std::mt19937_64 rng(0);
static u64 mul_mod(u64 a, u64 b, u64 mod) {
return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}
static u64 pow_mod(u64 a, u64 e, u64 mod) {
u64 r = 1 % mod;
a %= mod;
while (e > 0) {
if (e & 1ULL) {
r = mul_mod(r, a, mod);
}
a = mul_mod(a, a, mod);
e >>= 1ULL;
}
return r;
}
static bool is_prime64(u64 n) {
if (n < 2) return false;
for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (n % p == 0) return n == p;
}
u64 d = n - 1;
int s = 0;
while ((d & 1ULL) == 0) {
d >>= 1ULL;
++s;
}
for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
if (a % n == 0) continue;
u64 x = pow_mod(a, d, n);
if (x == 1 || x == n - 1) continue;
bool comp = true;
for (int r = 1; r < s; ++r) {
x = mul_mod(x, x, n);
if (x == n - 1) {
comp = false;
break;
}
}
if (comp) return false;
}
return true;
}
static u64 pollard_rho(u64 n) {
if ((n & 1ULL) == 0) return 2;
if (n % 3ULL == 0) return 3;
if (n % 5ULL == 0) return 5;
std::uniform_int_distribution<u64> dist(2, n - 2);
while (true) {
u64 c = dist(rng);
u64 x = dist(rng);
u64 y = x;
u64 d = 1;
while (d == 1) {
x = (mul_mod(x, x, n) + c) % n;
y = (mul_mod(y, y, n) + c) % n;
y = (mul_mod(y, y, n) + c) % n;
u64 diff = x > y ? x - y : y - x;
d = std::gcd(diff, n);
}
if (d != n) return d;
}
}
static void factor_rec(u64 n, std::map<u64, int>& out) {
if (n == 1) return;
if (is_prime64(n)) {
out[n]++;
return;
}
u64 d = pollard_rho(n);
factor_rec(d, out);
factor_rec(n / d, out);
}
static std::unordered_map<u64, std::vector<u64>> odd_div_cache;
static std::vector<u64> odd_divisors(u64 n) {
auto it = odd_div_cache.find(n);
if (it != odd_div_cache.end()) {
return it->second;
}
std::map<u64, int> fac;
factor_rec(n, fac);
std::vector<std::pair<u64, int>> vfac(fac.begin(), fac.end());
std::vector<u64> divs{1};
for (auto [p, e] : vfac) {
std::vector<u64> next;
next.reserve(divs.size() * static_cast<std::size_t>(e + 1));
u64 pe = 1;
for (int i = 0; i <= e; ++i) {
for (u64 d : divs) {
next.push_back(d * pe);
}
pe *= p;
}
divs.swap(next);
}
std::vector<u64> odd;
odd.reserve(divs.size());
for (u64 d : divs) {
if (d & 1ULL) odd.push_back(d);
}
std::sort(odd.begin(), odd.end());
odd_div_cache[n] = odd;
return odd;
}
static std::vector<u64> primes_mod4(int residue, int count) {
std::vector<u64> ps;
for (u64 x = 2; static_cast<int>(ps.size()) < count; ++x) {
if (x % 4 != static_cast<u64>(residue)) continue;
bool prime = true;
for (u64 p = 2; p * p <= x; ++p) {
if (x % p == 0) {
prime = false;
break;
}
}
if (prime) ps.push_back(x);
}
return ps;
}
static const std::vector<u64> P1 = primes_mod4(1, 80);
static const std::vector<u64> P3 = primes_mod4(3, 80);
static std::vector<Real> make_logs(const std::vector<u64>& ps) {
static const Real ln2 = log(Real(2));
std::vector<Real> logs(ps.size());
for (std::size_t i = 0; i < ps.size(); ++i) {
logs[i] = log(Real(ps[i])) / ln2;
}
return logs;
}
static const std::vector<Real> L1 = make_logs(P1);
static const std::vector<Real> L3 = make_logs(P3);
static const Real LOG2_2 = Real(1);
struct DfsKey {
u64 rem;
int idx;
u64 prevf;
bool operator<(const DfsKey& o) const {
if (rem != o.rem) return rem < o.rem;
if (idx != o.idx) return idx < o.idx;
return prevf < o.prevf;
}
};
static Rep solve_product(
u64 P,
const std::vector<u64>& primes,
const std::vector<Real>& logs,
std::unordered_map<u64, Rep>& top_cache
) {
if (P == 1) {
Rep r;
r.valid = true;
r.log2v = 0;
return r;
}
auto it_top = top_cache.find(P);
if (it_top != top_cache.end()) {
return it_top->second;
}
std::map<DfsKey, Rep> memo;
std::function<Rep(u64, int, u64)> dfs = [&](u64 rem, int idx, u64 prevf) -> Rep {
if (rem == 1) {
Rep base;
base.valid = true;
base.log2v = 0;
return base;
}
if (idx >= static_cast<int>(primes.size())) {
return Rep{};
}
DfsKey key{rem, idx, prevf};
auto it = memo.find(key);
if (it != memo.end()) {
return it->second;
}
Rep best;
auto divs = odd_divisors(rem);
for (auto rit = divs.rbegin(); rit != divs.rend(); ++rit) {
u64 f = *rit;
if (f == 1 || f > prevf) continue;
u64 e = (f - 1) / 2;
if (e == 0) continue;
Rep sub = dfs(rem / f, idx + 1, f);
if (!sub.valid) continue;
Rep cand;
cand.valid = true;
cand.log2v = logs[idx] * Real(e) + sub.log2v;
cand.exps.reserve(sub.exps.size() + 1);
cand.exps.push_back(e);
cand.exps.insert(cand.exps.end(), sub.exps.begin(), sub.exps.end());
if (!best.valid || cand.log2v < best.log2v) {
best = std::move(cand);
}
}
memo[key] = best;
return best;
};
Rep res = dfs(P, 0, P);
top_cache[P] = res;
return res;
}
static std::unordered_map<u64, Rep> cache_rep1;
static std::unordered_map<u64, Rep> cache_rep3;
static std::unordered_map<u64, BestD> cache_bestD;
static BestD best_for_D(u64 D) {
auto it = cache_bestD.find(D);
if (it != cache_bestD.end()) {
return it->second;
}
BestD best;
auto divs = odd_divisors(D);
for (u64 a2 : divs) {
u64 e2 = 0;
Real part2 = 0;
if (a2 > 1) {
e2 = (a2 + 1) / 2;
part2 = LOG2_2 * Real(e2);
}
u64 C = D / a2;
Rep rep3 = solve_product(C, P3, L3, cache_rep3);
if (!rep3.valid) continue;
Real cur = part2 + rep3.log2v;
if (!best.valid || cur < best.log2v) {
best.valid = true;
best.log2v = cur;
best.e2 = e2;
best.rep3 = std::move(rep3);
}
}
cache_bestD[D] = best;
return best;
}
struct QRep {
bool valid = false;
Real log2v = 0;
Rep rep1;
u64 e2 = 0;
Rep rep3;
};
static QRep Q_rep(u64 n) {
u64 S = n + 1;
QRep best;
auto divs = odd_divisors(S);
for (u64 B : divs) {
Rep rep1 = solve_product(B, P1, L1, cache_rep1);
if (!rep1.valid) continue;
u64 D = static_cast<u64>((static_cast<u128>(2) * S) / B - 1);
BestD bd = best_for_D(D);
if (!bd.valid) continue;
Real cur = rep1.log2v + bd.log2v;
if (!best.valid || cur < best.log2v) {
best.valid = true;
best.log2v = cur;
best.rep1 = std::move(rep1);
best.e2 = bd.e2;
best.rep3 = bd.rep3;
}
}
return best;
}
static u64 rep_mod(const Rep& rep, const std::vector<u64>& primes, u64 mod) {
u64 r = 1;
for (std::size_t i = 0; i < rep.exps.size(); ++i) {
r = mul_mod(r, pow_mod(primes[i], rep.exps[i], mod), mod);
}
return r;
}
static u64 Q_mod(u64 n) {
QRep qr = Q_rep(n);
u64 r = rep_mod(qr.rep1, P1, kMod);
r = mul_mod(r, pow_mod(2, qr.e2, kMod), kMod);
r = mul_mod(r, rep_mod(qr.rep3, P3, kMod), kMod);
return r;
}
static u128 rep_to_u128(const Rep& rep, const std::vector<u64>& primes) {
u128 x = 1;
for (std::size_t i = 0; i < rep.exps.size(); ++i) {
u64 e = rep.exps[i];
u128 p = primes[i];
u128 pw = 1;
while (e > 0) {
if (e & 1ULL) pw *= p;
p *= p;
e >>= 1ULL;
}
x *= pw;
}
return x;
}
static u64 Q_exact_small(u64 n) {
QRep qr = Q_rep(n);
u128 x = rep_to_u128(qr.rep1, P1);
x <<= qr.e2;
x *= rep_to_u128(qr.rep3, P3);
return static_cast<u64>(x);
}
int main() {
assert(Q_exact_small(5) == 15ULL);
assert(Q_exact_small(10) == 48ULL);
assert(Q_exact_small(1000) == 8064000ULL);
u64 ans = 0;
for (int k = 1; k <= 18; ++k) {
u64 n = 1;
for (int i = 0; i < k; ++i) n *= 10ULL;
ans += Q_mod(n);
ans %= kMod;
}
std::cout << ans << '\n';
return 0;
}
Python
import math
import random
from decimal import Decimal, getcontext
import sys
getcontext().prec = 100
kMod = 409120391
class Rep:
def __init__(self):
self.valid = False
self.log2v = Decimal(0)
self.exps = []
class BestD:
def __init__(self):
self.valid = False
self.log2v = Decimal(0)
self.e2 = 0
self.rep3 = Rep()
def mul_mod(a, b, mod):
return (a * b) % mod
def pow_mod(a, e, mod):
return pow(a, e, mod)
def is_prime64(n):
if n < 2: return False
for p in [2,3,5,7,11,13,17,19,23,29,31,37]:
if n % p == 0: return n == p
d = n - 1
s = 0
while (d & 1) == 0:
d >>= 1
s += 1
for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
if a % n == 0: continue
x = pow_mod(a, d, n)
if x == 1 or x == n - 1: continue
comp = True
for r in range(1, s):
x = mul_mod(x, x, n)
if x == n - 1:
comp = False
break
if comp: return False
return True
def pollard_rho(n):
if (n & 1) == 0: return 2
if n % 3 == 0: return 3
if n % 5 == 0: return 5
while True:
c = random.randint(2, n - 2)
x = random.randint(2, n - 2)
y = x
d = 1
while d == 1:
x = (mul_mod(x, x, n) + c) % n
y = (mul_mod(y, y, n) + c) % n
y = (mul_mod(y, y, n) + c) % n
diff = x - y if x > y else y - x
d = math.gcd(diff, n)
if d != n: return d
def factor_rec(n, out):
if n == 1: return
if is_prime64(n):
out[n] = out.get(n, 0) + 1
return
d = pollard_rho(n)
factor_rec(d, out)
factor_rec(n // d, out)
odd_div_cache = {}
def odd_divisors(n):
if n in odd_div_cache:
return odd_div_cache[n]
fac = {}
factor_rec(n, fac)
divs = [1]
for p, e in fac.items():
next_divs = []
pe = 1
for i in range(e + 1):
for d in divs:
next_divs.append(d * pe)
pe *= p
divs = next_divs
odd = sorted([d for d in divs if (d & 1)])
odd_div_cache[n] = odd
return odd
def primes_mod4(residue, count):
ps = []
x = 2
while len(ps) < count:
if x % 4 != residue:
x += 1
continue
prime = True
for p in range(2, math.isqrt(x) + 1):
if x % p == 0:
prime = False
break
if prime:
ps.append(x)
x += 1
return ps
P1 = primes_mod4(1, 80)
P3 = primes_mod4(3, 80)
def make_logs(ps):
ln2 = Decimal(2).ln()
return [Decimal(p).ln() / ln2 for p in ps]
L1 = make_logs(P1)
L3 = make_logs(P3)
LOG2_2 = Decimal(1)
def solve_product(P_val, primes, logs, top_cache):
if P_val == 1:
r = Rep()
r.valid = True
r.log2v = Decimal(0)
return r
if P_val in top_cache:
return top_cache[P_val]
memo = {}
def dfs(rem, idx, prevf):
if rem == 1:
base = Rep()
base.valid = True
base.log2v = Decimal(0)
return base
if idx >= len(primes):
return Rep()
key = (rem, idx, prevf)
if key in memo:
return memo[key]
best = Rep()
divs = odd_divisors(rem)
for f in reversed(divs):
if f == 1 or f > prevf: continue
e = (f - 1) // 2
if e == 0: continue
sub = dfs(rem // f, idx + 1, f)
if not sub.valid: continue
cand = Rep()
cand.valid = True
cand.log2v = logs[idx] * Decimal(e) + sub.log2v
cand.exps = [e] + sub.exps
if not best.valid or cand.log2v < best.log2v:
best = cand
memo[key] = best
return best
res = dfs(P_val, 0, P_val)
top_cache[P_val] = res
return res
cache_rep1 = {}
cache_rep3 = {}
cache_bestD = {}
def best_for_D(D):
if D in cache_bestD:
return cache_bestD[D]
best = BestD()
divs = odd_divisors(D)
for a2 in divs:
e2 = 0
part2 = Decimal(0)
if a2 > 1:
e2 = (a2 + 1) // 2
part2 = LOG2_2 * Decimal(e2)
C = D // a2
rep3 = solve_product(C, P3, L3, cache_rep3)
if not rep3.valid: continue
cur = part2 + rep3.log2v
if not best.valid or cur < best.log2v:
best.valid = True
best.log2v = cur
best.e2 = e2
best.rep3 = rep3
cache_bestD[D] = best
return best
class QRep:
def __init__(self):
self.valid = False
self.log2v = Decimal(0)
self.rep1 = Rep()
self.e2 = 0
self.rep3 = Rep()
def Q_rep(n):
S = n + 1
best = QRep()
divs = odd_divisors(S)
for B in divs:
rep1 = solve_product(B, P1, L1, cache_rep1)
if not rep1.valid: continue
D = (2 * S) // B - 1
bd = best_for_D(D)
if not bd.valid: continue
cur = rep1.log2v + bd.log2v
if not best.valid or cur < best.log2v:
best.valid = True
best.log2v = cur
best.rep1 = rep1
best.e2 = bd.e2
best.rep3 = bd.rep3
return best
def rep_mod(rep, primes, mod):
r = 1
for i in range(len(rep.exps)):
r = mul_mod(r, pow_mod(primes[i], rep.exps[i], mod), mod)
return r
def Q_mod(n):
qr = Q_rep(n)
r = rep_mod(qr.rep1, P1, kMod)
r = mul_mod(r, pow_mod(2, qr.e2, kMod), kMod)
r = mul_mod(r, rep_mod(qr.rep3, P3, kMod), kMod)
return r
def solve():
sys.setrecursionlimit(20000)
ans = 0
for k in range(1, 19):
ans += Q_mod(10**k)
ans %= kMod
return str(ans)
if __name__ == "__main__":
print(solve())
Java
import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.*;
public class Euler827 {
static final long kMod = 409120391L;
static final MathContext MC = new MathContext(100, RoundingMode.HALF_UP);
static class Rep {
boolean valid = false;
BigDecimal log2v = BigDecimal.ZERO;
ArrayList<Long> exps = new ArrayList<>();
}
static class BestD {
boolean valid = false;
BigDecimal log2v = BigDecimal.ZERO;
long e2 = 0;
Rep rep3;
}
static class QRep {
boolean valid = false;
BigDecimal log2v = BigDecimal.ZERO;
Rep rep1 = new Rep();
long e2 = 0;
Rep rep3 = new Rep();
}
static long mulMod(long a, long b, long mod) {
// use Math.multiplyHigh / 128-bit mul
long q = (long) ((double) a * b / mod);
long r = a * b - q * mod;
while (r < 0)
r += mod;
while (r >= mod)
r -= mod;
return r;
}
static long powMod(long a, long e, long mod) {
long r = 1 % mod;
a %= mod;
while (e > 0) {
if ((e & 1L) == 1L) {
r = mulMod(r, a, mod);
}
a = mulMod(a, a, mod);
e >>= 1L;
}
return r;
}
static boolean isPrime64(long n) {
if (n < 2)
return false;
long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
for (long p : smallPrimes) {
if (n % p == 0)
return n == p;
}
long d = n - 1;
int s = 0;
while ((d & 1L) == 0) {
d >>= 1L;
s++;
}
long[] bases = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
for (long a : bases) {
if (a % n == 0)
continue;
long x = powMod(a, d, n);
if (x == 1 || x == n - 1)
continue;
boolean comp = true;
for (int r = 1; r < s; ++r) {
x = mulMod(x, x, n);
if (x == n - 1) {
comp = false;
break;
}
}
if (comp)
return false;
}
return true;
}
static Random rng = new Random(0);
static long pollardRho(long n) {
if ((n & 1L) == 0)
return 2;
if (n % 3 == 0)
return 3;
if (n % 5 == 0)
return 5;
while (true) {
long c = 2 + (long) (rng.nextDouble() * (n - 3));
long x = 2 + (long) (rng.nextDouble() * (n - 3));
long y = x;
long d = 1;
while (d == 1) {
x = (mulMod(x, x, n) + c) % n;
y = (mulMod(y, y, n) + c) % n;
y = (mulMod(y, y, n) + c) % n;
long diff = x > y ? x - y : y - x;
d = gcd(diff, n);
}
if (d != n)
return d;
}
}
static long gcd(long a, long b) {
while (b != 0) {
long temp = b;
b = a % b;
a = temp;
}
return a;
}
static void factorRec(long n, Map<Long, Integer> out) {
if (n == 1)
return;
if (isPrime64(n)) {
out.put(n, out.getOrDefault(n, 0) + 1);
return;
}
long d = pollardRho(n);
factorRec(d, out);
factorRec(n / d, out);
}
static Map<Long, ArrayList<Long>> oddDivCache = new HashMap<>();
static ArrayList<Long> oddDivisors(long n) {
if (oddDivCache.containsKey(n)) {
return oddDivCache.get(n);
}
Map<Long, Integer> fac = new HashMap<>();
factorRec(n, fac);
ArrayList<Long> divs = new ArrayList<>();
divs.add(1L);
for (Map.Entry<Long, Integer> entry : fac.entrySet()) {
long p = entry.getKey();
int e = entry.getValue();
ArrayList<Long> next = new ArrayList<>();
long pe = 1;
for (int i = 0; i <= e; ++i) {
for (long d : divs) {
next.add(d * pe);
}
pe *= p;
}
divs = next;
}
ArrayList<Long> odd = new ArrayList<>();
for (long d : divs) {
if ((d & 1L) != 0) {
odd.add(d);
}
}
Collections.sort(odd);
oddDivCache.put(n, odd);
return odd;
}
static ArrayList<Long> primesMod4(int residue, int count) {
ArrayList<Long> ps = new ArrayList<>();
for (long x = 2; ps.size() < count; ++x) {
if (x % 4 != residue)
continue;
boolean prime = true;
for (long p = 2; p * p <= x; ++p) {
if (x % p == 0) {
prime = false;
break;
}
}
if (prime)
ps.add(x);
}
return ps;
}
static ArrayList<Long> P1;
static ArrayList<Long> P3;
static BigDecimal[] ln_table;
static void initLogs(int max_val) {
ln_table = new BigDecimal[max_val + 1];
ln_table[1] = BigDecimal.ZERO;
BigDecimal eps = new BigDecimal("1e-100");
for (int k = 2; k <= max_val; ++k) {
BigDecimal z = BigDecimal.ONE.divide(BigDecimal.valueOf(2L * k - 1), MC);
BigDecimal z2 = z.multiply(z, MC);
BigDecimal sum = BigDecimal.ZERO;
BigDecimal term = z;
int i = 0;
while (true) {
BigDecimal add = term.divide(BigDecimal.valueOf(2L * i + 1), MC);
if (add.compareTo(eps) < 0) {
break;
}
sum = sum.add(add, MC);
term = term.multiply(z2, MC);
i++;
}
sum = sum.multiply(BigDecimal.valueOf(2), MC);
ln_table[k] = ln_table[k - 1].add(sum, MC);
}
}
static BigDecimal[] L1;
static BigDecimal[] L3;
static final BigDecimal LOG2_2 = BigDecimal.ONE;
static BigDecimal[] makeLogs(ArrayList<Long> ps) {
BigDecimal ln2 = ln_table[2];
BigDecimal[] logs = new BigDecimal[ps.size()];
for (int i = 0; i < ps.size(); ++i) {
logs[i] = ln_table[ps.get(i).intValue()].divide(ln2, MC);
}
return logs;
}
static class DfsKey {
long rem;
int idx;
long prevf;
DfsKey(long rem, int idx, long prevf) {
this.rem = rem;
this.idx = idx;
this.prevf = prevf;
}
@Override
public boolean equals(Object o) {
if (this == o)
return true;
if (o == null || getClass() != o.getClass())
return false;
DfsKey dfsKey = (DfsKey) o;
return rem == dfsKey.rem && idx == dfsKey.idx && prevf == dfsKey.prevf;
}
@Override
public int hashCode() {
return Objects.hash(rem, idx, prevf);
}
}
static Rep solveProduct(long P_val, ArrayList<Long> primes, BigDecimal[] logs, Map<Long, Rep> top_cache) {
if (P_val == 1) {
Rep r = new Rep();
r.valid = true;
r.log2v = BigDecimal.ZERO;
return r;
}
if (top_cache.containsKey(P_val)) {
return top_cache.get(P_val);
}
Map<DfsKey, Rep> memo = new HashMap<>();
class DFS {
Rep dfs(long rem, int idx, long prevf) {
if (rem == 1) {
Rep base = new Rep();
base.valid = true;
base.log2v = BigDecimal.ZERO;
return base;
}
if (idx >= primes.size()) {
return new Rep();
}
DfsKey key = new DfsKey(rem, idx, prevf);
if (memo.containsKey(key)) {
return memo.get(key);
}
Rep best = new Rep();
ArrayList<Long> divs = oddDivisors(rem);
for (int i = divs.size() - 1; i >= 0; i--) {
long f = divs.get(i);
if (f == 1 || f > prevf)
continue;
long e = (f - 1) / 2;
if (e == 0)
continue;
Rep sub = dfs(rem / f, idx + 1, f);
if (!sub.valid)
continue;
Rep cand = new Rep();
cand.valid = true;
cand.log2v = logs[idx].multiply(BigDecimal.valueOf(e), MC).add(sub.log2v, MC);
cand.exps.add(e);
cand.exps.addAll(sub.exps);
if (!best.valid || cand.log2v.compareTo(best.log2v) < 0) {
best = cand;
}
}
memo.put(key, best);
return best;
}
}
DFS obj = new DFS();
Rep res = obj.dfs(P_val, 0, P_val);
top_cache.put(P_val, res);
return res;
}
static Map<Long, Rep> cacheRep1 = new HashMap<>();
static Map<Long, Rep> cacheRep3 = new HashMap<>();
static Map<Long, BestD> cacheBestD = new HashMap<>();
static BestD bestForD(long D) {
if (cacheBestD.containsKey(D)) {
return cacheBestD.get(D);
}
BestD best = new BestD();
ArrayList<Long> divs = oddDivisors(D);
for (long a2 : divs) {
long e2 = 0;
BigDecimal part2 = BigDecimal.ZERO;
if (a2 > 1) {
e2 = (a2 + 1) / 2;
part2 = LOG2_2.multiply(BigDecimal.valueOf(e2), MC);
}
long C = D / a2;
Rep rep3 = solveProduct(C, P3, L3, cacheRep3);
if (!rep3.valid)
continue;
BigDecimal cur = part2.add(rep3.log2v, MC);
if (!best.valid || cur.compareTo(best.log2v) < 0) {
best.valid = true;
best.log2v = cur;
best.e2 = e2;
best.rep3 = rep3;
}
}
cacheBestD.put(D, best);
return best;
}
static QRep Qrep(long n) {
long S = n + 1;
QRep best = new QRep();
ArrayList<Long> divs = oddDivisors(S);
for (long B : divs) {
Rep rep1 = solveProduct(B, P1, L1, cacheRep1);
if (!rep1.valid)
continue;
long D = (2 * S) / B - 1;
BestD bd = bestForD(D);
if (!bd.valid)
continue;
BigDecimal cur = rep1.log2v.add(bd.log2v, MC);
if (!best.valid || cur.compareTo(best.log2v) < 0) {
best.valid = true;
best.log2v = cur;
best.rep1 = rep1;
best.e2 = bd.e2;
best.rep3 = bd.rep3;
}
}
return best;
}
static long repMod(Rep rep, ArrayList<Long> primes, long mod) {
long r = 1;
for (int i = 0; i < rep.exps.size(); ++i) {
r = mulMod(r, powMod(primes.get(i), rep.exps.get(i), mod), mod);
}
return r;
}
static long Qmod(long n) {
QRep qr = Qrep(n);
long r = repMod(qr.rep1, P1, kMod);
r = mulMod(r, powMod(2, qr.e2, kMod), kMod);
r = mulMod(r, repMod(qr.rep3, P3, kMod), kMod);
return r;
}
public static String solve() {
P1 = primesMod4(1, 80);
P3 = primesMod4(3, 80);
long maxPrime = 0;
for (long p : P1)
maxPrime = Math.max(maxPrime, p);
for (long p : P3)
maxPrime = Math.max(maxPrime, p);
initLogs((int) maxPrime + 10);
L1 = makeLogs(P1);
L3 = makeLogs(P3);
long ans = 0;
for (int k = 1; k <= 18; ++k) {
long n = 1;
for (int i = 0; i < k; ++i)
n *= 10L;
ans += Qmod(n);
ans %= kMod;
}
return Long.toString(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}