Problem 486: Palindrome-containing Strings
View on Project EulerProject Euler Problem 486 Solution
EulerSolve provides an optimized solution for Project Euler Problem 486, Palindrome-containing Strings, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each \(n\), let \(A(n)\) be the number of binary strings of length \(n\) that avoid every palindromic substring of length at least \(5\). Then the number of length-\(n\) strings that do contain such a palindrome is \(2^n-A(n)\), and the cumulative function is $$F_5(n)=\sum_{k=1}^{n}\left(2^k-A(k)\right).$$ The required counting function is $$D(L)=\#\left\{n\in\mathbb{Z}:5\le n\le L,\ F_5(n)\equiv 0 \pmod{87654321}\right\}.$$ With \(L=10^{18}\), direct enumeration is hopeless. The implementations instead convert the problem into a periodic congruence in one auxiliary variable. Mathematical Approach The solution first counts palindrome-avoiding strings of one exact length, then turns that count into a closed form for \(F_5(n)\), and finally solves a modular periodicity problem. Step 1: Reduce the Forbidden Patterns to Lengths 5 and 6 A binary string contains a palindrome of length at least \(5\) if and only if it contains a palindrome of length exactly \(5\) or exactly \(6\). Every odd palindrome of length at least \(5\) has a centered palindrome of length \(5\), and every even palindrome of length at least \(6\) has a centered palindrome of length \(6\). So the avoidance problem is equivalent to forbidding only these two local patterns. Step 2: Count Exact-Length Avoiders with a Finite-State DP Let \(A(n)\) be the number of valid length-\(n\) strings....
Detailed mathematical approach
Problem Summary
For each \(n\), let \(A(n)\) be the number of binary strings of length \(n\) that avoid every palindromic substring of length at least \(5\). Then the number of length-\(n\) strings that do contain such a palindrome is \(2^n-A(n)\), and the cumulative function is
$$F_5(n)=\sum_{k=1}^{n}\left(2^k-A(k)\right).$$
The required counting function is
$$D(L)=\#\left\{n\in\mathbb{Z}:5\le n\le L,\ F_5(n)\equiv 0 \pmod{87654321}\right\}.$$
With \(L=10^{18}\), direct enumeration is hopeless. The implementations instead convert the problem into a periodic congruence in one auxiliary variable.
Mathematical Approach
The solution first counts palindrome-avoiding strings of one exact length, then turns that count into a closed form for \(F_5(n)\), and finally solves a modular periodicity problem.
Step 1: Reduce the Forbidden Patterns to Lengths 5 and 6
A binary string contains a palindrome of length at least \(5\) if and only if it contains a palindrome of length exactly \(5\) or exactly \(6\).
Every odd palindrome of length at least \(5\) has a centered palindrome of length \(5\), and every even palindrome of length at least \(6\) has a centered palindrome of length \(6\). So the avoidance problem is equivalent to forbidding only these two local patterns.
Step 2: Count Exact-Length Avoiders with a Finite-State DP
Let \(A(n)\) be the number of valid length-\(n\) strings. When a new bit \(x\) is appended, only suffixes ending at that new bit can create a new forbidden palindrome, so it is enough to remember the last five bits.
If the previous suffix is \(b_1b_2b_3b_4b_5\), then appending \(x\) creates a forbidden palindrome of length \(5\) exactly when
$$b_2=x \land b_3=b_5,$$
and it creates one of length \(6\) exactly when
$$b_1=x \land b_2=b_5 \land b_3=b_4.$$
So after length \(5\), the process is a DP on at most \(32\) suffix states. The exact totals begin as
$$A(5)=24,\quad A(6)=30,\quad A(7)=32,\quad A(8)=32,$$
and from \(n=9\) onward the totals enter a stable period of length \(6\):
$$A(9+6t+u)=a_u,\qquad (a_0,a_1,a_2,a_3,a_4,a_5)=(32,34,36,34,32,32).$$
One complete six-step block therefore contributes
$$a_0+a_1+a_2+a_3+a_4+a_5=200.$$
Step 3: Sum Those Exact Counts to Obtain \(F_5(n)\)
There are \(2^n\) binary strings of length \(n\), so the number that contain a qualifying palindrome is \(2^n-A(n)\). Summing over all lengths gives
$$F_5(n)=\sum_{k=1}^{n}\left(2^k-A(k)\right)=(2^{n+1}-2)-\sum_{k=1}^{n}A(k).$$
Write
$$n=9+6t+u,\qquad u\in\{0,1,2,3,4,5\}.$$
Because the exact-length avoidance counts are periodic, their cumulative sum is linear in \(t\):
$$\sum_{k=1}^{9+6t+u}A(k)=200t+(C_u-2),$$
with
$$\left(C_0,C_1,C_2,C_3,C_4,C_5\right)=\left(182,216,252,286,318,350\right).$$
Substituting into the previous identity yields the closed form used by all three implementations:
$$F_5(9+6t+u)=2^{10+u}64^t-(200t+C_u).$$
Step 4: Turn Divisibility into a Congruence in \(t\)
The target condition is \(F_5(n)\equiv 0 \pmod{87654321}\). For fixed \(u\), the closed form becomes
$$2^{10+u}64^t \equiv 200t+C_u \pmod{87654321}.$$
The modulus factors as
$$87654321=9\cdot 1997\cdot 4877.$$
So the implementations solve the congruence modulo \(9\), \(1997\), and \(4877\) separately. For one modulus \(m\), the exponential factor \(64^t\) has period \(\operatorname{ord}_m(64)\), while the linear term \(200t+C_u\) has period \(m\). Hence the local period is
$$P_m=\operatorname{lcm}\left(m,\operatorname{ord}_m(64)\right).$$
In this problem the three local periods are
$$P_9=9,\qquad P_{1997}=1993006,\qquad P_{4877}=11890126.$$
Step 5: Merge Local Residues with the Chinese Remainder Theorem
For each fixed \(u\), the implementation enumerates one full local period for each modulus factor and records every residue \(t\) that satisfies the congruence modulo that factor.
Two congruences
$$x\equiv a \pmod{m_1},\qquad x\equiv b \pmod{m_2}$$
are compatible exactly when
$$a\equiv b \pmod{\gcd(m_1,m_2)}.$$
Compatible residues are merged step by step into global residues modulo
$$P=\operatorname{lcm}(P_9,P_{1997},P_{4877}).$$
The important point is that the code never scans all \(P\) residue classes. It scans only the three local periods above and combines the matching classes constructively by CRT.
Worked Example
Take \(n=11\). Then \(n=9+6\cdot 0+2\), so \(t=0\) and \(u=2\). The closed form gives
$$F_5(11)=2^{12}-252=3844.$$
This matches the exact avoidance counts. Up to length \(11\), the valid-string totals are
$$A(1),\dots,A(11)=2,4,8,16,24,30,32,32,32,34,36,$$
so
$$\sum_{k=1}^{11}A(k)=250,$$
and therefore
$$F_5(11)=(2^{12}-2)-250=4094-250=3844.$$
Now move one full six-step block forward to \(n=17\). The avoidance sum grows by another \(200\), so
$$F_5(17)=2^{18}-(200+252)=261692.$$
This illustrates why large \(n\) are naturally parameterized by the pair \((t,u)\).
How the Code Works
The C++, Python, and Java implementations begin with a small finite-state DP over the last up to five bits. That DP produces the exact short values, including the checkpoints \(F_5(5)=8\), \(F_5(6)=42\), and \(F_5(11)=3844\), and it also confirms the six-term avoidance cycle used by the closed form.
Next, for every \(u\in\{0,1,2,3,4,5\}\) and every modulus factor \(m\in\{9,1997,4877\}\), the implementation scans one local period \(P_m\) and records the residues satisfying
$$2^{10+u}64^t \equiv 200t+C_u \pmod m.$$
Those local residue sets are merged with CRT into sorted global residue lists. For a query limit \(L\), the code handles \(n=5,6,7,8\) directly and, for each \(u\), counts how many
$$t\le \left\lfloor\frac{L-(9+u)}{6}\right\rfloor$$
fall into the stored global residue classes. With sorted residues, the final count is computed as full periods plus a tail found by binary search.
Complexity Analysis
Let \(R_{m,u}\) be the number of valid residues for modulus factor \(m\) in class \(u\). The local scans cost
$$O\left(\sum_{u=0}^{5}\left(P_9+P_{1997}+P_{4877}\right)\right),$$
which is a one-time preprocessing step independent of \(L\). The CRT phase depends on the residue-set sizes, not on the enormous global period \(P\). After preprocessing, evaluating \(D(L)\) only requires six limit conversions and counts in sorted arrays, so the dependence on \(L\) is effectively \(O(1)\) with small logarithmic factors from binary search. Memory usage is proportional to the total number of stored residue classes.
Footnotes and References
- Problem page: https://projecteuler.net/problem=486
- Palindrome: Wikipedia — Palindrome
- Chinese remainder theorem: Wikipedia — Chinese remainder theorem
- Multiplicative order: Wikipedia — Multiplicative order
- Finite-state machine: Wikipedia — Finite-state machine
Problem 486 source code
C++
#include <algorithm>
#include <array>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <utility>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;
using u128 = unsigned __int128;
constexpr u64 kMainMod = 87'654'321ULL;
constexpr u64 kDefaultLimit = 1'000'000'000'000'000'000ULL;
constexpr u64 kCheckpointL1 = 10'000'000ULL;
constexpr u64 kCheckpointL2 = 5'000'000'000ULL;
// For n = 9 + 6t + u (u in [0,5]):
// F5(n) = 2^(10+u) * 64^t - (200t + C[u]).
constexpr std::array<u64, 6> kC = {182ULL, 216ULL, 252ULL, 286ULL, 318ULL, 350ULL};
constexpr std::array<u64, 3> kPrimePowerFactors = {9ULL, 1997ULL, 4877ULL};
struct Options {
u64 limit = kDefaultLimit;
bool allow_multithreading = true;
bool run_checkpoints = true;
unsigned requested_threads = 0;
};
struct ModData {
u64 mod = 1;
u64 period = 1;
std::array<std::vector<u32>, 6> residues_by_u;
};
bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const u64 digit = static_cast<u64>(c - '0');
if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
return false;
}
parsed = parsed * 10ULL + digit;
}
value = parsed;
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u64 parsed = 0;
if (!parse_u64_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
u64 limit = 0;
if (parse_u64_after_prefix(arg, "--limit=", limit)) {
options.limit = limit;
continue;
}
unsigned threads = 0;
if (parse_unsigned_after_prefix(arg, "--threads=", threads)) {
options.requested_threads = threads;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return true;
}
u64 mul_mod(u64 a, u64 b, u64 mod) {
return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % mod);
}
u64 pow_mod(u64 base, u64 exp, u64 mod) {
u64 result = 1 % mod;
u64 cur = base % mod;
u64 e = exp;
while (e > 0) {
if (e & 1ULL) {
result = mul_mod(result, cur, mod);
}
cur = mul_mod(cur, cur, mod);
e >>= 1ULL;
}
return result;
}
u64 gcd_u64(u64 a, u64 b) {
while (b != 0) {
const u64 t = a % b;
a = b;
b = t;
}
return a;
}
u64 lcm_u64(u64 a, u64 b) {
return (a / gcd_u64(a, b)) * b;
}
std::vector<u64> distinct_prime_factors(u64 n) {
std::vector<u64> factors;
for (u64 p = 2; p * p <= n; ++p) {
if (n % p != 0) {
continue;
}
factors.push_back(p);
while (n % p == 0) {
n /= p;
}
}
if (n > 1) {
factors.push_back(n);
}
return factors;
}
u64 euler_phi(u64 n) {
if (n == 0) {
return 0;
}
u64 result = n;
u64 x = n;
for (u64 p = 2; p * p <= x; ++p) {
if (x % p != 0) {
continue;
}
while (x % p == 0) {
x /= p;
}
result -= result / p;
}
if (x > 1) {
result -= result / x;
}
return result;
}
u64 multiplicative_order(u64 a, u64 mod) {
if (gcd_u64(a, mod) != 1ULL) {
return 0;
}
u64 order = euler_phi(mod);
const std::vector<u64> factors = distinct_prime_factors(order);
for (const u64 p : factors) {
while (order % p == 0ULL && pow_mod(a, order / p, mod) == 1ULL) {
order /= p;
}
}
return order;
}
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload) {
if (!allow_multithreading || workload < 2) {
return 1;
}
unsigned threads = requested_threads;
if (threads == 0) {
threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
}
threads = std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
return threads;
}
// We keep enough prefix values to validate all F5 sample checkpoints.
std::vector<u64> compute_F5_prefix(int max_n) {
std::vector<u64> f(static_cast<std::size_t>(max_n) + 1ULL, 0ULL);
// State is last up to 5 bits.
std::vector<std::pair<u32, u64>> states;
states.emplace_back(0U, 1ULL);
f[0] = 0ULL;
for (int n = 1; n <= max_n; ++n) {
u64 next_count[32] = {0ULL};
for (const auto& [bits, count] : states) {
for (u32 add = 0; add <= 1; ++add) {
const u32 combined = (bits << 1U) | add;
bool ok = true;
if (n >= 5) {
const u32 v5 = combined & 31U;
const u32 b0 = (v5 >> 0U) & 1U;
const u32 b1 = (v5 >> 1U) & 1U;
const u32 b3 = (v5 >> 3U) & 1U;
const u32 b4 = (v5 >> 4U) & 1U;
if (b0 == b4 && b1 == b3) {
ok = false;
}
}
if (ok && n >= 6) {
const u32 v6 = combined & 63U;
const u32 b0 = (v6 >> 0U) & 1U;
const u32 b1 = (v6 >> 1U) & 1U;
const u32 b2 = (v6 >> 2U) & 1U;
const u32 b3 = (v6 >> 3U) & 1U;
const u32 b4 = (v6 >> 4U) & 1U;
const u32 b5 = (v6 >> 5U) & 1U;
if (b0 == b5 && b1 == b4 && b2 == b3) {
ok = false;
}
}
if (!ok) {
continue;
}
const u32 next_bits = (n <= 5) ? (combined & ((1U << n) - 1U)) : (combined & 31U);
next_count[next_bits] += count;
}
}
states.clear();
states.reserve(32);
u64 total_a = 0ULL;
if (n <= 5) {
const u32 limit = 1U << n;
for (u32 mask = 0; mask < limit; ++mask) {
if (next_count[mask] != 0ULL) {
states.emplace_back(mask, next_count[mask]);
total_a += next_count[mask];
}
}
} else {
for (u32 mask = 0; mask < 32U; ++mask) {
if (next_count[mask] != 0ULL) {
states.emplace_back(mask, next_count[mask]);
total_a += next_count[mask];
}
}
}
const u64 total_strings = (n < 63) ? (1ULL << n) : 0ULL;
f[static_cast<std::size_t>(n)] = f[static_cast<std::size_t>(n - 1)] + (total_strings - total_a);
}
return f;
}
std::vector<u32> compute_residues_for_mod_u(u64 mod, u64 period, int u) {
std::vector<u32> residues;
const u64 lhs_mul = pow_mod(2ULL, static_cast<u64>(10 + u), mod);
u64 pow64 = 1ULL % mod;
u64 rhs = kC[static_cast<std::size_t>(u)] % mod;
residues.reserve(static_cast<std::size_t>(period / mod + 8ULL));
for (u64 t = 0; t < period; ++t) {
if (mul_mod(lhs_mul, pow64, mod) == rhs) {
residues.push_back(static_cast<u32>(t));
}
pow64 = mul_mod(pow64, 64ULL % mod, mod);
rhs += 200ULL;
rhs %= mod;
}
return residues;
}
u64 mod_inverse(u64 a, u64 mod) {
i64 t = 0;
i64 new_t = 1;
i64 r = static_cast<i64>(mod);
i64 new_r = static_cast<i64>(a % mod);
while (new_r != 0) {
const i64 q = r / new_r;
const i64 temp_t = t - q * new_t;
t = new_t;
new_t = temp_t;
const i64 temp_r = r - q * new_r;
r = new_r;
new_r = temp_r;
}
if (r != 1) {
return 0;
}
i64 x = t;
x %= static_cast<i64>(mod);
if (x < 0) {
x += static_cast<i64>(mod);
}
return static_cast<u64>(x);
}
struct CRTPrecompute {
u64 m1 = 1;
u64 m2 = 1;
u64 g = 1;
u64 m1_div_g = 1;
u64 m2_div_g = 1;
u64 lcm = 1;
u64 inv_m1_div_g_mod_m2_div_g = 0;
};
CRTPrecompute build_crt_precompute(u64 m1, u64 m2) {
CRTPrecompute data;
data.m1 = m1;
data.m2 = m2;
data.g = gcd_u64(m1, m2);
data.m1_div_g = m1 / data.g;
data.m2_div_g = m2 / data.g;
data.lcm = data.m1_div_g * m2;
data.inv_m1_div_g_mod_m2_div_g = mod_inverse(data.m1_div_g % data.m2_div_g, data.m2_div_g);
return data;
}
u64 crt_merge_pair(u64 a, u64 b, const CRTPrecompute& crt) {
// Assumes compatibility: (a - b) % g == 0.
const i64 diff = static_cast<i64>(b) - static_cast<i64>(a);
const i64 reduced = diff / static_cast<i64>(crt.g);
i64 k = reduced % static_cast<i64>(crt.m2_div_g);
if (k < 0) {
k += static_cast<i64>(crt.m2_div_g);
}
const u64 ku = static_cast<u64>(k);
const u64 t = mul_mod(ku, crt.inv_m1_div_g_mod_m2_div_g, crt.m2_div_g);
const u64 x = static_cast<u64>((static_cast<u128>(a) +
static_cast<u128>(crt.m1) * static_cast<u128>(t)) %
crt.lcm);
return x;
}
std::vector<u64> combine_residues_for_u(const std::vector<u32>& mod1997,
u64 period1997,
const std::vector<u32>& mod4877,
u64 period4877,
const std::vector<u32>& mod9,
u64 period9,
u64& out_global_period) {
const CRTPrecompute crt12 = build_crt_precompute(period1997, period4877);
const CRTPrecompute crt129 = build_crt_precompute(crt12.lcm, period9);
out_global_period = crt129.lcm;
std::vector<u32> mod4877_even;
std::vector<u32> mod4877_odd;
mod4877_even.reserve(mod4877.size() / 2 + 1);
mod4877_odd.reserve(mod4877.size() / 2 + 1);
for (const u32 value : mod4877) {
if ((value & 1U) == 0U) {
mod4877_even.push_back(value);
} else {
mod4877_odd.push_back(value);
}
}
std::vector<u64> combined;
combined.reserve(mod1997.size() * (mod4877.size() / 2 + 1) * std::max<std::size_t>(1, mod9.size()));
for (const u32 a32 : mod1997) {
const u64 a = static_cast<u64>(a32);
const std::vector<u32>& candidates = ((a32 & 1U) == 0U) ? mod4877_even : mod4877_odd;
for (const u32 b32 : candidates) {
const u64 b = static_cast<u64>(b32);
const u64 merged12 = crt_merge_pair(a, b, crt12);
for (const u32 c32 : mod9) {
const u64 c = static_cast<u64>(c32);
const u64 merged = crt_merge_pair(merged12, c, crt129);
combined.push_back(merged);
}
}
}
std::sort(combined.begin(), combined.end());
combined.erase(std::unique(combined.begin(), combined.end()), combined.end());
return combined;
}
u64 count_from_residues(const std::vector<u64>& residues, u64 period, u64 t_limit) {
if (residues.empty()) {
return 0ULL;
}
const u64 full = t_limit / period;
const u64 rem = t_limit % period;
const auto it = std::upper_bound(residues.begin(), residues.end(), rem);
const u64 in_tail = static_cast<u64>(std::distance(residues.begin(), it));
return full * static_cast<u64>(residues.size()) + in_tail;
}
class Euler486Solver {
public:
explicit Euler486Solver(const Options& options)
: allow_multithreading_(options.allow_multithreading),
requested_threads_(options.requested_threads) {
initialize_mod_data();
build_mod_residue_tables();
build_global_residue_tables();
}
u64 D(u64 limit) const {
if (limit < 5ULL) {
return 0ULL;
}
u64 answer = 0ULL;
const u64 small_end = std::min<u64>(limit, 8ULL);
for (u64 n = 5ULL; n <= small_end; ++n) {
if (small_f5_mod_[static_cast<std::size_t>(n)] == 0ULL) {
++answer;
}
}
if (limit < 9ULL) {
return answer;
}
for (int u = 0; u < 6; ++u) {
const u64 n0 = 9ULL + static_cast<u64>(u);
if (n0 > limit) {
continue;
}
const u64 t_limit = (limit - n0) / 6ULL;
answer += count_from_residues(global_residues_by_u_[static_cast<std::size_t>(u)],
global_period_,
t_limit);
}
return answer;
}
bool run_checkpoints() const {
bool ok = true;
// F5 sample values.
ok &= checkpoint_equal("F5(4)", small_f5_[4], 0ULL);
ok &= checkpoint_equal("F5(5)", small_f5_[5], 8ULL);
ok &= checkpoint_equal("F5(6)", small_f5_[6], 42ULL);
ok &= checkpoint_equal("F5(11)", small_f5_[11], 3'844ULL);
// D(L) sample values.
ok &= checkpoint_equal("D(10^7)", D(kCheckpointL1), 0ULL);
ok &= checkpoint_equal("D(5*10^9)", D(kCheckpointL2), 51ULL);
// Additional structural checkpoint used by the derivation.
static constexpr std::array<u64, 6> expected_a_cycle = {32ULL, 34ULL, 36ULL, 34ULL, 32ULL, 32ULL};
for (int i = 0; i < 6; ++i) {
const u64 n = 9ULL + static_cast<u64>(i);
ok &= checkpoint_equal((std::string("A(") + std::to_string(n) + ")").c_str(),
small_a_[static_cast<std::size_t>(n)],
expected_a_cycle[static_cast<std::size_t>(i)]);
}
return ok;
}
private:
bool allow_multithreading_ = true;
unsigned requested_threads_ = 0;
std::array<ModData, 3> mod_data_;
std::array<std::vector<u64>, 6> global_residues_by_u_;
u64 global_period_ = 1ULL;
std::vector<u64> small_f5_;
std::vector<u64> small_a_;
std::vector<u64> small_f5_mod_;
static bool checkpoint_equal(const char* label, u64 got, u64 expected) {
const bool ok = (got == expected);
std::cout << "[checkpoint] " << label << " = " << got;
if (!ok) {
std::cout << " (expected " << expected << ")";
}
std::cout << '\n';
return ok;
}
void initialize_mod_data() {
for (std::size_t i = 0; i < mod_data_.size(); ++i) {
mod_data_[i].mod = kPrimePowerFactors[i];
const u64 order = multiplicative_order(64ULL % mod_data_[i].mod, mod_data_[i].mod);
mod_data_[i].period = lcm_u64(mod_data_[i].mod, order);
}
small_f5_ = compute_F5_prefix(20);
small_f5_mod_.assign(small_f5_.size(), 0ULL);
for (std::size_t i = 0; i < small_f5_.size(); ++i) {
small_f5_mod_[i] = small_f5_[i] % kMainMod;
}
// Also keep A(n) values for validation checks.
small_a_.assign(21, 0ULL);
small_a_[0] = 1ULL;
for (std::size_t n = 1; n <= 20; ++n) {
const u64 total = (n < 63) ? (1ULL << n) : 0ULL;
small_a_[n] = total - (small_f5_[n] - small_f5_[n - 1]);
}
}
void build_mod_residue_tables() {
struct Task {
int mod_index = 0;
int u = 0;
};
std::vector<Task> tasks;
tasks.reserve(18);
for (int mi = 0; mi < 3; ++mi) {
for (int u = 0; u < 6; ++u) {
tasks.push_back(Task{mi, u});
}
}
const unsigned threads = choose_thread_count(allow_multithreading_,
requested_threads_,
tasks.size());
if (threads == 1U) {
for (const Task& task : tasks) {
ModData& md = mod_data_[static_cast<std::size_t>(task.mod_index)];
md.residues_by_u[static_cast<std::size_t>(task.u)] =
compute_residues_for_mod_u(md.mod, md.period, task.u);
}
return;
}
std::atomic<std::size_t> next_task(0);
std::vector<std::thread> workers;
workers.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
workers.emplace_back([&]() {
while (true) {
const std::size_t index = next_task.fetch_add(1);
if (index >= tasks.size()) {
break;
}
const Task task = tasks[index];
ModData& md = mod_data_[static_cast<std::size_t>(task.mod_index)];
md.residues_by_u[static_cast<std::size_t>(task.u)] =
compute_residues_for_mod_u(md.mod, md.period, task.u);
}
});
}
for (auto& worker : workers) {
worker.join();
}
}
void build_global_residue_tables() {
// Index mapping by modulus value.
const int idx9 = 0;
const int idx1997 = 1;
const int idx4877 = 2;
std::array<std::pair<int, std::vector<u64>>, 6> temp;
for (int u = 0; u < 6; ++u) {
temp[static_cast<std::size_t>(u)].first = u;
}
const unsigned threads = choose_thread_count(allow_multithreading_, requested_threads_, 6);
if (threads == 1U) {
for (int u = 0; u < 6; ++u) {
u64 period = 1;
temp[static_cast<std::size_t>(u)].second = combine_residues_for_u(
mod_data_[idx1997].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx1997].period,
mod_data_[idx4877].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx4877].period,
mod_data_[idx9].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx9].period,
period);
global_period_ = period;
}
} else {
std::atomic<int> next_u(0);
std::atomic<u64> period_seen(0ULL);
std::vector<std::thread> workers;
workers.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
workers.emplace_back([&]() {
while (true) {
const int u = next_u.fetch_add(1);
if (u >= 6) {
break;
}
u64 period = 1ULL;
auto combined = combine_residues_for_u(
mod_data_[idx1997].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx1997].period,
mod_data_[idx4877].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx4877].period,
mod_data_[idx9].residues_by_u[static_cast<std::size_t>(u)],
mod_data_[idx9].period,
period);
temp[static_cast<std::size_t>(u)].second = std::move(combined);
u64 expected = 0ULL;
if (period_seen.compare_exchange_strong(expected, period)) {
// first writer sets the period
} else {
if (expected != period) {
std::cerr << "Inconsistent global period across residues.\n";
std::terminate();
}
}
}
});
}
for (auto& worker : workers) {
worker.join();
}
global_period_ = period_seen.load();
}
for (const auto& item : temp) {
global_residues_by_u_[static_cast<std::size_t>(item.first)] = item.second;
}
}
};
} // namespace
int main(int argc, char** argv) {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
Options options;
if (!parse_arguments(argc, argv, options)) {
std::cerr << "Usage: ./Euler486 [--limit=N] [--threads=T] [--single-thread] "
"[--skip-checkpoints]\n";
return 1;
}
Euler486Solver solver(options);
bool ok = true;
if (options.run_checkpoints) {
ok = solver.run_checkpoints();
}
const u64 answer = solver.D(options.limit);
std::cout << "D(" << options.limit << ") = " << answer << '\n';
if (!ok) {
std::cerr << "Checkpoint validation failed.\n";
return 2;
}
return 0;
}
Python
def solve():
MOD = 87654321; LIMIT = 10**18
C = [182,216,252,286,318,350]; PPF = [9,1997,4877]
def mm(a,b,m): return a*b%m
def pm(b,e,m):
r=1; b%=m
while e>0:
if e&1: r=r*b%m
b=b*b%m; e>>=1
return r
def gcd(a,b):
while b: a,b=b,a%b
return a
def lcm(a,b): return a//gcd(a,b)*b
def euler_phi(n):
r=n; x=n
p=2
while p*p<=x:
if x%p==0:
while x%p==0: x//=p
r-=r//p
p+=1
if x>1: r-=r//x
return r
def mult_order(a, m):
if gcd(a,m)!=1: return 0
o=euler_phi(m); x=o
p=2
while p*p<=x:
if x%p==0:
while o%p==0 and pm(a,o//p,m)==1: o//=p
p+=1
while x%p==0: x//=p
if x>1:
while o%x==0 and pm(a,o//x,m)==1: o//=x
return o
def mod_inv(a,m):
t,nt,r,nr=0,1,m,a%m
while nr:
q=r//nr; t,nt=nt,t-q*nt; r,nr=nr,r-q*nr
if r!=1: return 0
return t%m
def crt_merge(a,m1,b,m2):
g=gcd(m1,m2); m1g=m1//g; m2g=m2//g; lc=m1g*m2
inv=mod_inv(m1g%m2g,m2g)
diff=b-a; dg=diff//g
k=dg%m2g; t=k*inv%m2g
x=(a+m1*t)%lc
return x,lc
# Compute mod data
mod_data = []
for pp in PPF:
o = mult_order(64%pp, pp); period = lcm(pp, o)
residues = [[] for _ in range(6)]
for u in range(6):
lhs_mul = pm(2, 10+u, pp)
p64 = 1%pp; rhs = C[u]%pp
res = []
for t in range(period):
if mm(lhs_mul, p64, pp) == rhs: res.append(t)
p64 = mm(p64, 64%pp, pp); rhs = (rhs+200)%pp
residues[u] = res
mod_data.append((pp, period, residues))
# Combine residues for each u
global_period = 1
global_res = [None]*6
for u in range(6):
r1997 = mod_data[1][2][u]; p1997 = mod_data[1][1]
r4877 = mod_data[2][2][u]; p4877 = mod_data[2][1]
r9 = mod_data[0][2][u]; p9 = mod_data[0][1]
# CRT merge 1997 and 4877
combined12 = []
even4877 = [x for x in r4877 if x%2==0]
odd4877 = [x for x in r4877 if x%2==1]
for a in r1997:
cands = even4877 if a%2==0 else odd4877
for b in cands:
g = gcd(p1997, p4877)
if (b-a)%g != 0: continue
m12, lc12 = crt_merge(a, p1997, b, p4877)
combined12.append((m12, lc12))
# merge with mod9
combined = set()
for m12, lc12 in combined12:
for c in r9:
g = gcd(lc12, p9)
if (c-m12)%g != 0: continue
m, lc = crt_merge(m12, lc12, c, p9)
combined.add(m)
global_period = lc
global_res[u] = sorted(combined)
# Small cases (n=5..8)
def compute_f5_prefix():
f = [0]*21; states = [(0, 1)]
for n in range(1, 21):
nc = [0]*32
for bits,cnt in states:
for add in range(2):
cb = (bits<<1)|add; ok=True
if n>=5:
v5=cb&31; b0=v5&1; b1=(v5>>1)&1; b3=(v5>>3)&1; b4=(v5>>4)&1
if b0==b4 and b1==b3: ok=False
if ok and n>=6:
v6=cb&63
if (v6&1)==((v6>>5)&1) and ((v6>>1)&1)==((v6>>4)&1) and ((v6>>2)&1)==((v6>>3)&1): ok=False
if not ok: continue
nb = (cb&((1<<n)-1)) if n<=5 else (cb&31)
nc[nb] += cnt
states = [(m,nc[m]) for m in range(32 if n>5 else (1<<n)) if nc[m]>0]
total_a = sum(c for _,c in states)
total_str = 1<<n if n<63 else 0
f[n] = f[n-1] + (total_str - total_a)
return f
f5 = compute_f5_prefix()
f5m = [x%MOD for x in f5]
# Count
import bisect
def count_res(residues, period, t_limit):
if not residues: return 0
full = t_limit // period; rem = t_limit % period
tail = bisect.bisect_right(residues, rem)
return full*len(residues)+tail
answer = 0
for n in range(5, min(LIMIT, 8)+1):
if f5m[n] == 0: answer += 1
if LIMIT >= 9:
for u in range(6):
n0 = 9+u
if n0 > LIMIT: continue
t_limit = (LIMIT-n0)//6
answer += count_res(global_res[u], global_period, t_limit)
return str(answer)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
public class Euler486 {
private static final long[] K_C = { 182L, 216L, 252L, 286L, 318L, 350L };
private static final long[] K_PRIME_POWER_FACTORS = { 9L, 1997L, 4877L };
private static final long K_MAIN_MOD = 87654321L;
private static long gcdU64(long a, long b) {
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return a;
}
private static long lcmU64(long a, long b) {
return (a / gcdU64(a, b)) * b;
}
private static List<Long> distinctPrimeFactors(long n) {
List<Long> factors = new ArrayList<>();
for (long p = 2; p * p <= n; ++p) {
if (n % p == 0) {
factors.add(p);
while (n % p == 0) {
n /= p;
}
}
}
if (n > 1) {
factors.add(n);
}
return factors;
}
private static long eulerPhi(long n) {
if (n == 0)
return 0;
long res = n;
long x = n;
for (long p = 2; p * p <= x; ++p) {
if (x % p == 0) {
while (x % p == 0)
x /= p;
res -= res / p;
}
}
if (x > 1)
res -= res / x;
return res;
}
private static long powMod(long base, long exp, long mod) {
long res = 1 % mod;
long cur = base % mod;
while (exp > 0) {
if ((exp & 1) != 0) {
res = (res * cur) % mod;
}
cur = (cur * cur) % mod;
exp >>= 1;
}
return res;
}
private static long multiplicativeOrder(long a, long mod) {
if (gcdU64(a, mod) != 1)
return 0;
long order = eulerPhi(mod);
List<Long> factors = distinctPrimeFactors(order);
for (long p : factors) {
while (order % p == 0 && powMod(a, order / p, mod) == 1) {
order /= p;
}
}
return order;
}
private static long[] computeF5Prefix(int maxN) {
long[] f = new long[maxN + 1];
List<long[]> states = new ArrayList<>();
states.add(new long[] { 0, 1 }); // bits, count
for (int n = 1; n <= maxN; ++n) {
long[] nextCount = new long[32];
for (long[] st : states) {
long bits = st[0];
long count = st[1];
for (long add = 0; add <= 1; ++add) {
long combined = (bits << 1) | add;
boolean ok = true;
if (n >= 5) {
long v5 = combined & 31;
long b0 = (v5 >> 0) & 1, b1 = (v5 >> 1) & 1, b3 = (v5 >> 3) & 1, b4 = (v5 >> 4) & 1;
if (b0 == b4 && b1 == b3)
ok = false;
}
if (ok && n >= 6) {
long v6 = combined & 63;
long b0 = (v6 >> 0) & 1, b1 = (v6 >> 1) & 1, b2 = (v6 >> 2) & 1;
long b3 = (v6 >> 3) & 1, b4 = (v6 >> 4) & 1, b5 = (v6 >> 5) & 1;
if (b0 == b5 && b1 == b4 && b2 == b3)
ok = false;
}
if (!ok)
continue;
int nextBits = (n <= 5) ? (int) (combined & ((1 << n) - 1)) : (int) (combined & 31);
nextCount[nextBits] += count;
}
}
states.clear();
long totalA = 0;
for (int mask = 0; mask < 32; ++mask) {
if (nextCount[mask] != 0) {
states.add(new long[] { mask, nextCount[mask] });
totalA += nextCount[mask];
}
}
long totalStrings = (n < 63) ? (1L << n) : 0L;
f[n] = f[n - 1] + (totalStrings - totalA);
}
return f;
}
private static List<Long> computeResiduesForModU(long mod, long period, int u) {
List<Long> residues = new ArrayList<>();
long lhsMul = powMod(2, 10 + u, mod);
long pow64 = 1 % mod;
long rhs = K_C[u] % mod;
for (long t = 0; t < period; ++t) {
if ((lhsMul * pow64) % mod == rhs) {
residues.add(t);
}
pow64 = (pow64 * 64) % mod;
rhs = (rhs + 200) % mod;
}
return residues;
}
private static long modInverse(long a, long mod) {
long t = 0, newT = 1;
long r = mod, newR = a % mod;
while (newR != 0) {
long q = r / newR;
long tempT = t - q * newT;
t = newT;
newT = tempT;
long tempR = r - q * newR;
r = newR;
newR = tempR;
}
if (r != 1)
return 0;
long res = t % mod;
if (res < 0)
res += mod;
return res;
}
static class CRTPrecompute {
long m1, m2, g, m1DivG, m2DivG, lcm, invM1DivGModM2DivG;
CRTPrecompute(long m1, long m2) {
this.m1 = m1;
this.m2 = m2;
this.g = gcdU64(m1, m2);
this.m1DivG = m1 / g;
this.m2DivG = m2 / g;
this.lcm = m1DivG * m2;
this.invM1DivGModM2DivG = modInverse(this.m1DivG % this.m2DivG, this.m2DivG);
}
}
private static long crtMergePair(long a, long b, CRTPrecompute crt) {
long diff = b - a;
long reduced = diff / crt.g;
long k = reduced % crt.m2DivG;
if (k < 0)
k += crt.m2DivG;
long t = (k * crt.invM1DivGModM2DivG) % crt.m2DivG;
return (a + crt.m1 * t) % crt.lcm;
}
private static class MergeResult {
List<Long> residues;
long period;
MergeResult(List<Long> residues, long period) {
this.residues = residues;
this.period = period;
}
}
private static MergeResult combineResiduesForU(List<Long> mod1997, long p1997,
List<Long> mod4877, long p4877,
List<Long> mod9, long p9) {
CRTPrecompute crt12 = new CRTPrecompute(p1997, p4877);
CRTPrecompute crt129 = new CRTPrecompute(crt12.lcm, p9);
long globalPeriod = crt129.lcm;
List<Long> mod4877Even = new ArrayList<>();
List<Long> mod4877Odd = new ArrayList<>();
for (long v : mod4877) {
if (v % 2 == 0)
mod4877Even.add(v);
else
mod4877Odd.add(v);
}
List<Long> combined = new ArrayList<>();
for (long a : mod1997) {
List<Long> candidates = (a % 2 == 0) ? mod4877Even : mod4877Odd;
for (long b : candidates) {
long merged12 = crtMergePair(a, b, crt12);
for (long c : mod9) {
long merged = crtMergePair(merged12, c, crt129);
combined.add(merged);
}
}
}
Collections.sort(combined);
List<Long> unique = new ArrayList<>();
for (Long v : combined) {
if (unique.isEmpty() || !unique.get(unique.size() - 1).equals(v)) {
unique.add(v);
}
}
return new MergeResult(unique, globalPeriod);
}
private static long countFromResidues(List<Long> residues, long period, long tLimit) {
if (residues.isEmpty())
return 0;
long full = tLimit / period;
long rem = tLimit % period;
int inTail = Collections.binarySearch(residues, rem);
if (inTail < 0) {
inTail = -(inTail + 1);
} else {
inTail++;
while (inTail < residues.size() && residues.get(inTail) == rem) {
inTail++;
}
}
return full * residues.size() + inTail;
}
public static void main(String[] args) {
long[] mods = K_PRIME_POWER_FACTORS;
long[] periods = new long[3];
List<List<Long>>[] residuesByU = new List[3];
for (int i = 0; i < 3; i++) {
long order = multiplicativeOrder(64 % mods[i], mods[i]);
periods[i] = lcmU64(mods[i], order);
residuesByU[i] = new ArrayList<>();
for (int u = 0; u < 6; u++) {
residuesByU[i].add(computeResiduesForModU(mods[i], periods[i], u));
}
}
long[] smallF5 = computeF5Prefix(20);
long[] smallF5Mod = new long[smallF5.length];
for (int i = 0; i < smallF5.length; i++) {
smallF5Mod[i] = smallF5[i] % K_MAIN_MOD;
}
List<List<Long>> globalResiduesByU = new ArrayList<>();
long globalPeriod = 1;
for (int u = 0; u < 6; u++) {
MergeResult res = combineResiduesForU(
residuesByU[1].get(u), periods[1],
residuesByU[2].get(u), periods[2],
residuesByU[0].get(u), periods[0]);
globalResiduesByU.add(res.residues);
globalPeriod = res.period;
}
long limit = 1000000000000000000L;
long answer = 0;
long smallEnd = Math.min(limit, 8);
for (int n = 5; n <= smallEnd; n++) {
if (smallF5Mod[n] == 0)
answer++;
}
if (limit >= 9) {
for (int u = 0; u < 6; u++) {
long n0 = 9 + u;
if (n0 > limit)
continue;
long tLimit = (limit - n0) / 6;
answer += countFromResidues(globalResiduesByU.get(u), globalPeriod, tLimit);
}
}
System.out.println(answer);
}
}