Problem 343: Fractional Sequences
View on Project EulerProject Euler Problem 343 Solution
EulerSolve provides an optimized solution for Project Euler Problem 343, Fractional Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Start from the reduced fraction \(1/n\). At each step replace \(x/y\) by \((x+1)/(y-1)\) and reduce again. If the first time the denominator becomes 1 the numerator is \(f(n)\), the problem asks for $$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$ A direct simulation for every cube is far too slow, so the key is to rewrite \(f(n)\) using the prime factors of \(n+1\). Mathematical Approach Write any reduced state as \(a/b\) with \(a,b>0\), and define its state sum by $$s=a+b.$$ Initially we have \(a=1\), \(b=n\), hence \(s=n+1\). One Step in Terms of the State Sum From the reduced fraction \(a/b\), the unreduced next fraction is \((a+1)/(b-1)\). If $$d=\gcd(a+1,b-1),$$ then after reduction $$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$ Since \(b=s-a\), we get $$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$ Therefore the next state sum is $$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$ So every reduction replaces the current sum \(s\) by a divisor of \(s\). Also, because \(a/b\) is reduced, \(\gcd(a,b)=1\), and thus \(\gcd(a,s)=1\). A Phase Always Removes the Smallest Prime Factor Suppose a phase starts at the special state \(1/(s-1)\). Let \(p\) be the smallest prime factor of \(s\). For every integer \(t\) with \(2 \le t \lt p\), we have \(\gcd(t,s)=1\). Hence no reduction happens while the numerator runs through \(1,2,\dots,p-1\)....
Detailed mathematical approach
Problem Summary
Start from the reduced fraction \(1/n\). At each step replace \(x/y\) by \((x+1)/(y-1)\) and reduce again. If the first time the denominator becomes 1 the numerator is \(f(n)\), the problem asks for
$$\sum_{k=1}^{L} f(k^3),\qquad L=2\cdot 10^6.$$
A direct simulation for every cube is far too slow, so the key is to rewrite \(f(n)\) using the prime factors of \(n+1\).
Mathematical Approach
Write any reduced state as \(a/b\) with \(a,b>0\), and define its state sum by
$$s=a+b.$$
Initially we have \(a=1\), \(b=n\), hence \(s=n+1\).
One Step in Terms of the State Sum
From the reduced fraction \(a/b\), the unreduced next fraction is \((a+1)/(b-1)\). If
$$d=\gcd(a+1,b-1),$$
then after reduction
$$a'=\frac{a+1}{d},\qquad b'=\frac{b-1}{d}.$$
Since \(b=s-a\), we get
$$d=\gcd(a+1,s-a-1)=\gcd(a+1,s).$$
Therefore the next state sum is
$$s'=a'+b'=\frac{s}{\gcd(a+1,s)}.$$
So every reduction replaces the current sum \(s\) by a divisor of \(s\). Also, because \(a/b\) is reduced, \(\gcd(a,b)=1\), and thus \(\gcd(a,s)=1\).
A Phase Always Removes the Smallest Prime Factor
Suppose a phase starts at the special state \(1/(s-1)\). Let \(p\) be the smallest prime factor of \(s\). For every integer \(t\) with \(2 \le t \lt p\), we have \(\gcd(t,s)=1\). Hence no reduction happens while the numerator runs through \(1,2,\dots,p-1\).
When the unreduced numerator becomes \(p\), the reduction factor is
$$d=\gcd(p,s)=p.$$
The fraction at that moment is
$$\frac{p}{s-p}\to \frac{1}{s/p-1}.$$
So one full phase transforms the state sum by
$$s\longmapsto \frac{s}{p},$$
which means: a phase strips off the smallest prime factor of the current sum and resets the numerator to 1.
Closed Form for \(f(n)\)
Starting from \(s_0=n+1\), repeated phases remove the prime factors of \(n+1\) from smallest to largest. After all but the largest prime factor have been removed, the state sum becomes
$$P^+(n+1)=\operatorname{LPF}(n+1),$$
where \(\operatorname{LPF}\) denotes the largest prime factor.
If the current sum \(s\) is prime, then for every \(1\le a\le s-2\) we have \(\gcd(a+1,s)=1\), so no further reductions occur before the denominator reaches 1. Beginning from \(1/(s-1)\), the chain ends at \((s-1)/1\). Therefore
$$\boxed{f(n)=\operatorname{LPF}(n+1)-1.}$$
Worked Example: \(n=20\)
The iteration is
$$\frac{1}{20}\to\frac{2}{19}\to\frac{3}{18}=\frac{1}{6}\to\frac{2}{5}\to\frac{3}{4}\to\frac{4}{3}\to\frac{5}{2}\to\frac{6}{1}.$$
So \(f(20)=6\). Since \(20+1=21=3\cdot 7\), the closed form gives
$$\operatorname{LPF}(21)-1=7-1=6,$$
which matches the simulation exactly.
Cubic Specialization
For the required sum, set \(n=k^3\):
$$f(k^3)=\operatorname{LPF}(k^3+1)-1.$$
Using the factorization of a sum of cubes,
$$k^3+1=(k+1)(k^2-k+1).$$
If we define
$$Q_k=k^2-k+1,$$
then
$$\operatorname{LPF}(k^3+1)=\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr).$$
This remains true even when the two factors share a factor 3, because the largest prime factor of a product is still the maximum of the largest prime factors of its factors. Hence
$$\sum_{k=1}^{L}f(k^3)=\sum_{k=1}^{L}\left(\max\bigl(\operatorname{LPF}(k+1),\operatorname{LPF}(Q_k)\bigr)-1\right).$$
Sieving \(\operatorname{LPF}(k+1)\)
A standard largest-prime-factor sieve on \([2,L+1]\) fills an array \(\operatorname{LPF}(m)\) for every \(m\le L+1\). This immediately gives all values \(\operatorname{LPF}(k+1)\).
Sieving \(\operatorname{LPF}(Q_k)\)
For each \(k\), initialize
$$\mathrm{rem}[k]=Q_k,\qquad \mathrm{lpfQ}[k]=1.$$
Now fix a prime \(p\). We need to know for which \(k\) the congruence
$$Q_k\equiv0\pmod p$$
holds. Multiplying by 4 gives
$$4k^2-4k+4\equiv0\pmod p\qquad\Longleftrightarrow\qquad (2k-1)^2\equiv -3\pmod p.$$
Therefore, for \(p>3\), roots exist exactly when \(-3\) is a quadratic residue modulo \(p\). If \(s^2\equiv -3\pmod p\), then
$$k\equiv \frac{1\pm s}{2}\pmod p.$$
The code finds the square root \(s\) with the Tonelli-Shanks algorithm. The prime \(p=3\) is special: the only class is \(k\equiv2\pmod 3\). The prime \(2\) never divides \(Q_k\), because \(k^2-k=k(k-1)\) is always even, so \(Q_k\) is always odd.
For every matching residue class, the code visits all corresponding \(k\), divides out all powers of \(p\) from \(\mathrm{rem}[k]\), and stores \(p\) as the current largest prime factor. After all primes \(p\le L\) are processed, any leftover \(\mathrm{rem}[k]>1\) must itself be prime and greater than \(L\); otherwise a composite remainder would have a prime divisor \(\le \sqrt{Q_k} \lt L\), which would already have been removed. That leftover is therefore the final large prime factor of \(Q_k\).
Worked Example: \(k=4\)
Here
$$k^3+1=65=5\cdot 13,$$
so
$$f(64)=\operatorname{LPF}(65)-1=13-1=12.$$
The decomposition gives \(k+1=5\) and \(Q_4=13\), so the larger prime factor is indeed \(13\).
How the Code Works
The implementation has two sieve layers. First it fills lpf_small for all integers up to \(L+1\).
Then it stores every \(Q_k\) in the residual array rem and processes primes one by one. For each prime
it computes the admissible residue classes of \(k\), strips that prime from all matching residuals, and updates
lpf_q[k]. Finally it combines lpf_small[k+1] with lpf_q[k], subtracts 1, and
accumulates the sum in a 128-bit total. The C++ version checks \(f(20)=6\),
\(\sum_{k\le100} f(k^3)=118937\), and also matches a brute-force computation for small limits.
Complexity Analysis
The sieve for \(\operatorname{LPF}(k+1)\) is near-linear in practice. The \(Q_k\)-sieve visits about \(L/p\) or \(2L/p\) positions per relevant prime, so its dominant work is also near-linear up to the usual \(\sum_{p\le L}1/p\) factor, i.e. \(O(L\log\log L)\) with extra polylogarithmic cost for modular square roots. Memory usage is \(O(L)\) for the factor tables and residual buffer.
Footnotes and References
- Problem page: https://projecteuler.net/problem=343
- Largest prime factor: Wikipedia - Prime factor
- Tonelli-Shanks algorithm: Wikipedia - Tonelli-Shanks algorithm
- Modular square roots: Wikipedia - Quadratic residue
- Sum of cubes factorization: Wikipedia - Polynomial factorization identities
Problem 343 source code
C++
#include <algorithm>
#include <cstdint>
#include <functional>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u128 = __uint128_t;
struct Options {
int limit = 2'000'000;
bool run_checkpoints = true;
};
bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
int parsed = 0;
for (char ch : tail) {
if (ch < '0' || ch > '9') {
return false;
}
parsed = parsed * 10 + static_cast<int>(ch - '0');
}
value = parsed;
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.limit >= 1;
}
u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}
u64 pow_mod(u64 base, u64 exp, const u64 mod) {
u64 result = 1ULL;
base %= mod;
while (exp > 0ULL) {
if (exp & 1ULL) {
result = mul_mod(result, base, mod);
}
base = mul_mod(base, base, mod);
exp >>= 1ULL;
}
return result;
}
u64 tonelli_shanks(const u64 n, const u64 p) {
if (n == 0ULL) {
return 0ULL;
}
if (p == 2ULL) {
return n;
}
if (pow_mod(n, (p - 1ULL) / 2ULL, p) != 1ULL) {
return 0ULL;
}
if ((p & 3ULL) == 3ULL) {
return pow_mod(n, (p + 1ULL) / 4ULL, p);
}
u64 q = p - 1ULL;
int s = 0;
while ((q & 1ULL) == 0ULL) {
q >>= 1ULL;
++s;
}
u64 z = 2ULL;
while (pow_mod(z, (p - 1ULL) / 2ULL, p) != p - 1ULL) {
++z;
}
u64 c = pow_mod(z, q, p);
u64 x = pow_mod(n, (q + 1ULL) / 2ULL, p);
u64 t = pow_mod(n, q, p);
int m = s;
while (t != 1ULL) {
int i = 1;
u64 t2i = mul_mod(t, t, p);
while (i < m && t2i != 1ULL) {
t2i = mul_mod(t2i, t2i, p);
++i;
}
const u64 b = pow_mod(c, 1ULL << (m - i - 1), p);
x = mul_mod(x, b, p);
t = mul_mod(t, mul_mod(b, b, p), p);
c = mul_mod(b, b, p);
m = i;
}
return x;
}
std::string to_string_u128(u128 value) {
if (value == 0) {
return "0";
}
std::string s;
while (value > 0) {
const int d = static_cast<int>(value % 10);
s.push_back(static_cast<char>('0' + d));
value /= 10;
}
std::reverse(s.begin(), s.end());
return s;
}
u128 solve(const int limit) {
std::vector<int> lpf_small(static_cast<std::size_t>(limit + 2), 0);
std::vector<int> primes;
for (int i = 2; i <= limit + 1; ++i) {
if (lpf_small[static_cast<std::size_t>(i)] == 0) {
for (int j = i; j <= limit + 1; j += i) {
lpf_small[static_cast<std::size_t>(j)] = i;
}
if (i <= limit) {
primes.push_back(i);
}
}
}
std::vector<u64> rem(static_cast<std::size_t>(limit + 1), 0ULL);
std::vector<u64> lpf_q(static_cast<std::size_t>(limit + 1), 1ULL);
for (int k = 1; k <= limit; ++k) {
const u64 kk = static_cast<u64>(k);
rem[static_cast<std::size_t>(k)] = kk * kk - kk + 1ULL;
}
for (int p : primes) {
if (p == 2) {
continue;
}
std::vector<int> roots;
if (p == 3) {
roots.push_back(2);
} else {
const u64 mod = static_cast<u64>(p);
const u64 neg_three = (mod + mod - 3ULL) % mod;
if (pow_mod(neg_three, (mod - 1ULL) / 2ULL, mod) != 1ULL) {
continue;
}
const u64 sq = tonelli_shanks(neg_three, mod);
const u64 inv2 = (mod + 1ULL) / 2ULL;
const int r1 = static_cast<int>(mul_mod((1ULL + sq) % mod, inv2, mod));
const int r2 = static_cast<int>(mul_mod((1ULL + mod - sq) % mod, inv2, mod));
roots.push_back(r1);
if (r2 != r1) {
roots.push_back(r2);
}
}
for (int root : roots) {
int k = root;
if (k == 0) {
k += p;
}
for (; k <= limit; k += p) {
u64& value = rem[static_cast<std::size_t>(k)];
if (value % static_cast<u64>(p) != 0ULL) {
continue;
}
while (value % static_cast<u64>(p) == 0ULL) {
value /= static_cast<u64>(p);
}
lpf_q[static_cast<std::size_t>(k)] = static_cast<u64>(p);
}
}
}
u128 total = 0;
for (int k = 1; k <= limit; ++k) {
const u64 rem_val = rem[static_cast<std::size_t>(k)];
if (rem_val > 1ULL) {
lpf_q[static_cast<std::size_t>(k)] = std::max(lpf_q[static_cast<std::size_t>(k)], rem_val);
}
const u64 lpf_kp1 = static_cast<u64>(lpf_small[static_cast<std::size_t>(k + 1)]);
const u64 largest = std::max(lpf_kp1, lpf_q[static_cast<std::size_t>(k)]);
total += static_cast<u128>(largest - 1ULL);
}
return total;
}
u64 brute_f_of_k(const int k) {
u64 x = 1ULL;
u64 y = static_cast<u64>(k);
while (y != 1ULL) {
++x;
--y;
const u64 g = std::gcd(x, y);
x /= g;
y /= g;
}
return x;
}
u128 brute_sum_k_cubed(const int limit) {
u128 s = 0;
for (int k = 1; k <= limit; ++k) {
const int kc = k * k * k;
s += static_cast<u128>(brute_f_of_k(kc));
}
return s;
}
bool run_checkpoints() {
if (brute_f_of_k(20) != 6ULL) {
std::cerr << "Checkpoint failed: f(20)=6" << '\n';
return false;
}
if (to_string_u128(solve(100)) != "118937") {
std::cerr << "Checkpoint failed: sum_{k<=100} f(k^3)=118937" << '\n';
return false;
}
if (solve(300) != brute_sum_k_cubed(300)) {
std::cerr << "Checkpoint failed: brute-force cross-check for limit=300" << '\n';
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 2;
}
std::cout << to_string_u128(solve(options.limit)) << '\n';
return 0;
}
Python
import math
def mul_mod(a, b, mod):
return (a * b) % mod
def pow_mod(base, exp, mod):
return pow(base, exp, mod)
def tonelli_shanks(n, p):
if n == 0:
return 0
if p == 2:
return n
if pow_mod(n, (p - 1) // 2, p) != 1:
return 0
if (p % 4) == 3:
return pow_mod(n, (p + 1) // 4, p)
q = p - 1
s = 0
while (q % 2) == 0:
q //= 2
s += 1
z = 2
while pow_mod(z, (p - 1) // 2, p) != p - 1:
z += 1
c = pow_mod(z, q, p)
x = pow_mod(n, (q + 1) // 2, p)
t = pow_mod(n, q, p)
m = s
while t != 1:
i = 1
t2i = mul_mod(t, t, p)
while i < m and t2i != 1:
t2i = mul_mod(t2i, t2i, p)
i += 1
b = pow_mod(c, 1 << (m - i - 1), p)
x = mul_mod(x, b, p)
t = mul_mod(t, mul_mod(b, b, p), p)
c = mul_mod(b, b, p)
m = i
return x
def solve_343(limit):
lpf_small = [0] * (limit + 2)
primes = []
for i in range(2, limit + 2):
if lpf_small[i] == 0:
for j in range(i, limit + 2, i):
lpf_small[j] = i
if i <= limit:
primes.append(i)
rem = [0] * (limit + 1)
lpf_q = [1] * (limit + 1)
for k in range(1, limit + 1):
rem[k] = k * k - k + 1
for p in primes:
if p == 2:
continue
roots = []
if p == 3:
roots.append(2)
else:
neg_three = (2 * p - 3) % p
if pow_mod(neg_three, (p - 1) // 2, p) != 1:
continue
sq = tonelli_shanks(neg_three, p)
inv2 = (p + 1) // 2
r1 = ((1 + sq) % p * inv2) % p
r2 = ((1 + p - sq) % p * inv2) % p
roots.append(r1)
if r2 != r1:
roots.append(r2)
for root in roots:
k = root
if k == 0:
k += p
for j in range(k, limit + 1, p):
val = rem[j]
if val % p != 0:
continue
while val % p == 0:
val //= p
rem[j] = val
lpf_q[j] = p
total = 0
for k in range(1, limit + 1):
rem_val = rem[k]
if rem_val > 1:
lpf_q[k] = max(lpf_q[k], rem_val)
lpf_kp1 = lpf_small[k + 1]
largest = max(lpf_kp1, lpf_q[k])
total += largest - 1
return total
def solve():
ans = solve_343(2000000)
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler343 {
static long mulMod(long a, long b, long mod) {
return (long) ((((java.math.BigInteger.valueOf(a)).multiply(java.math.BigInteger.valueOf(b)))
.remainder(java.math.BigInteger.valueOf(mod))).longValue());
}
static long powMod(long base, long exp, long mod) {
long result = 1;
base %= mod;
while (exp > 0) {
if ((exp & 1) == 1) {
result = mulMod(result, base, mod);
}
base = mulMod(base, base, mod);
exp >>= 1;
}
return result;
}
static long tonelliShanks(long n, long p) {
if (n == 0)
return 0;
if (p == 2)
return n;
if (powMod(n, (p - 1) / 2, p) != 1)
return 0;
if ((p % 4) == 3)
return powMod(n, (p + 1) / 4, p);
long q = p - 1;
int s = 0;
while ((q % 2) == 0) {
q /= 2;
s++;
}
long z = 2;
while (powMod(z, (p - 1) / 2, p) != p - 1) {
z++;
}
long c = powMod(z, q, p);
long x = powMod(n, (q + 1) / 2, p);
long t = powMod(n, q, p);
int m = s;
while (t != 1) {
int i = 1;
long t2i = mulMod(t, t, p);
while (i < m && t2i != 1) {
t2i = mulMod(t2i, t2i, p);
i++;
}
long b = powMod(c, 1L << (m - i - 1), p);
x = mulMod(x, b, p);
t = mulMod(t, mulMod(b, b, p), p);
c = mulMod(b, b, p);
m = i;
}
return x;
}
static long solve(int limit) {
int[] lpfSmall = new int[limit + 2];
List<Integer> primes = new ArrayList<>();
for (int i = 2; i <= limit + 1; i++) {
if (lpfSmall[i] == 0) {
for (int j = i; j <= limit + 1; j += i) {
lpfSmall[j] = i;
}
if (i <= limit) {
primes.add(i);
}
}
}
long[] rem = new long[limit + 1];
long[] lpfQ = new long[limit + 1];
Arrays.fill(lpfQ, 1);
for (int k = 1; k <= limit; k++) {
long kk = k;
rem[k] = kk * kk - kk + 1;
}
for (int p : primes) {
if (p == 2)
continue;
List<Integer> roots = new ArrayList<>();
if (p == 3) {
roots.add(2);
} else {
long mod = p;
long negThree = (mod + mod - 3) % mod;
if (powMod(negThree, (mod - 1) / 2, mod) != 1)
continue;
long sq = tonelliShanks(negThree, mod);
long inv2 = (mod + 1) / 2;
int r1 = (int) mulMod((1 + sq) % mod, inv2, mod);
int r2 = (int) mulMod((1 + mod - sq) % mod, inv2, mod);
roots.add(r1);
if (r2 != r1) {
roots.add(r2);
}
}
for (int root : roots) {
int k = root;
if (k == 0)
k += p;
for (; k <= limit; k += p) {
long value = rem[k];
if (value % p != 0)
continue;
while (value % p == 0) {
value /= p;
}
rem[k] = value;
lpfQ[k] = p;
}
}
}
long total = 0;
for (int k = 1; k <= limit; k++) {
long remVal = rem[k];
if (remVal > 1) {
lpfQ[k] = Math.max(lpfQ[k], remVal);
}
long lpfKp1 = lpfSmall[k + 1];
long largest = Math.max(lpfKp1, lpfQ[k]);
total += (largest - 1);
}
return total;
}
public static String solve() {
long ans = solve(2000000);
return String.valueOf(ans);
}
public static void main(String[] args) {
System.out.println(solve());
}
}