Problem 954: Heptaphobia
View on Project EulerProject Euler Problem 954 Solution
EulerSolve provides an optimized solution for Project Euler Problem 954, Heptaphobia, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Call a positive integer heptaphobic if it is not divisible by 7 and it still stays non-divisible by 7 after every valid swap of two decimal digits. Swaps that would move 0 to the front are ignored, because they would not produce another number of the same length. The goal is to count all heptaphobic numbers below \(10^{13}\). A direct test of every integer is hopeless. The implementations instead count numbers by length, fix the nonzero residue \(m=n\bmod 7\), and ask which digit patterns make it impossible for any admissible swap to force the residue to \(0\). Mathematical Approach Take a length-\(L\) number with digits \(d_0,d_1,\dots,d_{L-1}\): $$n=\sum_{i=0}^{L-1} d_i\,10^{L-1-i}.$$ Modulo 7, only the positional weights matter, so define $$w_i\equiv 10^{L-1-i}\pmod 7.$$ Then \(n\bmod 7\) is simply \(\sum d_iw_i\bmod 7\). The whole method comes from understanding how a swap changes this residue and from compressing positions with the same weight....
Detailed mathematical approach
Problem Summary
Call a positive integer heptaphobic if it is not divisible by 7 and it still stays non-divisible by 7 after every valid swap of two decimal digits. Swaps that would move 0 to the front are ignored, because they would not produce another number of the same length. The goal is to count all heptaphobic numbers below \(10^{13}\).
A direct test of every integer is hopeless. The implementations instead count numbers by length, fix the nonzero residue \(m=n\bmod 7\), and ask which digit patterns make it impossible for any admissible swap to force the residue to \(0\).
Mathematical Approach
Take a length-\(L\) number with digits \(d_0,d_1,\dots,d_{L-1}\):
$$n=\sum_{i=0}^{L-1} d_i\,10^{L-1-i}.$$
Modulo 7, only the positional weights matter, so define
$$w_i\equiv 10^{L-1-i}\pmod 7.$$
Then \(n\bmod 7\) is simply \(\sum d_iw_i\bmod 7\). The whole method comes from understanding how a swap changes this residue and from compressing positions with the same weight.
Residue Change Under a Digit Swap
If digits at positions \(i\) and \(j\) are exchanged, every other term is unchanged, so the residue changes by
$$\Delta\equiv d_jw_i+d_iw_j-d_iw_i-d_jw_j\equiv (d_j-d_i)(w_i-w_j)\pmod 7.$$
Therefore a number with residue \(m\in\{1,\dots,6\}\) becomes divisible by 7 after that swap exactly when
$$m+(d_j-d_i)(w_i-w_j)\equiv 0\pmod 7.$$
This already shows why the code never needs the full decimal values of the digits. Only their residues modulo 7 matter.
Six Positional Weight Classes
Because \(10\equiv 3\pmod 7\) and \(10^6\equiv 1\pmod 7\), the powers of 10 repeat modulo 7 with period 6:
$$10^k\bmod 7\in\{1,3,2,6,4,5\},\qquad 10^{k+6}\equiv 10^k\pmod 7.$$
So, for a fixed length \(L\), every non-leading position belongs to one of six weight classes. All positions in the same class have the same multiplier \(w_i\), and swapping two digits inside one class does nothing modulo 7 because then \(w_i-w_j\equiv 0\pmod 7\).
That means a whole class can be summarized by four facts:
the set of residues that occur there, the set of nonzero residues that occur there, the sum of the class digits modulo 7, and the number of ordered digit tuples that realize that same summary.
This compression is exactly what the problem needs. Pairwise swap safety depends only on which residues appear in each class, while the overall residue of the number depends only on the weighted class sums.
Forbidden Swaps Become Shifted Residue Intersections
Take two different weight classes with weights \(a\) and \(b\), and suppose the whole number has residue \(m\). If one class contains a digit residue \(r\) and the other contains a residue \(s\), the swap between those positions is forbidden precisely when
$$m+(s-r)(a-b)\equiv 0\pmod 7.$$
Since \(a\not\equiv b\pmod 7\), the difference \(a-b\) has an inverse modulo 7. Define
$$t(a,b,m)\equiv -m(a-b)^{-1}\pmod 7.$$
Then the forbidden relation becomes
$$s\equiv r+t(a,b,m)\pmod 7.$$
So two class summaries are compatible exactly when the residue set of the first class, shifted by \(t(a,b,m)\), is disjoint from the residue set of the second class. That is the key simplification used by the implementations: every swap condition turns into a tiny 7-bit set-intersection test.
A concrete example helps. If \(m=5\), \(a=3\), and \(b=1\), then \(a-b\equiv 2\), its inverse is \(4\), and
$$t(3,1,5)\equiv -5\cdot 4\equiv 1\pmod 7.$$
So any residue \(r\) appearing in the weight-3 class forbids the residue \(r+1\pmod 7\) in the weight-1 class.
The Leading Digit Has a Special Constraint
Swaps involving the first digit are asymmetric, because a later 0 cannot be moved to the front. Such a swap would create a leading zero and is excluded by the problem statement. For that reason the leading digit is handled separately, chosen from the residues of the digits \(1,\dots,9\), and compared with the other classes using only the nonzero residues present there.
There is also one important no-op case. If the leading weight reappears later in the number, then swapping with that position cannot change the residue at all, because the two weights are equal modulo 7. Since the original residue \(m\) is already nonzero, such swaps can never create a multiple of 7.
The Global Condition Becomes a Small Search
Fix the length \(L\), the target residue \(m\in\{1,\dots,6\}\), and the leading-digit residue \(r_0\). If the first position has weight \(w_0\), the remaining weight classes must contribute
$$\sum_C w(C)\,\sigma(C)\equiv m-r_0w_0\pmod 7,$$
where \(\sigma(C)\) is the stored digit-sum residue for class \(C\).
So the problem is now finite and discrete: choose one summary for each active weight class, keep every pair of chosen summaries compatible, and make the weighted sum land on the required residue. The implementations solve this with a depth-first search that always branches on the class with the fewest remaining summaries, which gives strong pruning.
Worked Example: Two-Digit Numbers
For \(L=2\), the weights are \(3\) and \(1\), so a number \(10a+b\) has residue
$$m\equiv 3a+b\pmod 7.$$
The only possible digit swap produces \(10b+a\), and the residue change is
$$\Delta\equiv (b-a)(3-1)\equiv 2(b-a)\pmod 7.$$
Take \(12\). Its residue is \(3\cdot 1+2\equiv 5\pmod 7\), but swapping gives \(21\equiv 0\pmod 7\), so \(12\) is not heptaphobic. By contrast, \(13\equiv 6\pmod 7\) and its swap \(31\equiv 3\pmod 7\), so that swap does not violate the rule. The full algorithm uses exactly this modular test, but it applies it class by class instead of number by number.
How the Code Works
Precomputed Residue Summaries
The C++, Python, and Java implementations first precompute three small constant tables: the six-term cycle of \(10^k\bmod 7\), the number of possible leading digits for each residue class, and the catalogue of summary types for a weight class of a given size. Each summary type records which residues appear, which nonzero residues appear, what the class contributes to the digit sum modulo 7, and how many ordered digit assignments produce that same summary.
For the actual bound \(10^{13}\), every non-leading weight class has size 0, 1, or 2, so these catalogues stay tiny: there is 1 summary of size 0, 8 summaries of size 1, and 35 summaries of size 2.
Counting One \((L,m)\) Instance
For a fixed length \(L\) and target residue \(m\), the implementation counts how many non-leading positions belong to each of the six weights. Every active class receives its list of possible summaries together with its weighted modular contribution. It also precomputes, for every ordered pair of distinct weight classes, which summaries can coexist without creating a forbidden swap.
Next, the leading-digit residue is chosen. That immediately removes any class summaries that would allow a forbidden nonzero swap with the front position. After that, a recursive search picks one summary for one class at a time, intersects the remaining candidate sets with the precomputed compatibility tables, updates the running residue sum, and multiplies by the number of concrete digit assignments represented by each chosen summary.
Building the Final Count
The result is summed over all nonzero residues \(m=1,\dots,6\) for a fixed length, and then over all lengths from 1 through 13. Small built-in checks on short ranges verify the logic before the final count below \(10^{13}\) is produced.
Complexity Analysis
A brute-force approach would inspect almost \(10^{13}\) integers and, for each one, consider up to \(\binom{13}{2}\) digit swaps. The implemented method never iterates over concrete numbers at that scale.
For each pair \((L,m)\) with \(1\le L\le 13\) and \(m\in\{1,\dots,6\}\), there are at most six non-leading weight classes. Each class uses only a small precomputed summary catalogue, and for this bound no class has more than 35 relevant summaries. The search depth is therefore at most six, with heavy pruning from pairwise compatibility and the modular sum condition. Memory usage is \(O(1)\) relative to the numeric bound, because only small residue tables and candidate masks are stored.
In practice the running time is dominated by a modest branch-and-prune search across the \((L,m)\) cases, which is why the method is fast enough even though direct enumeration would be completely infeasible.
Footnotes and References
- Problem page: https://projecteuler.net/problem=954
- Modular arithmetic: Wikipedia - Modular arithmetic
- Multiplicative order: Wikipedia - Multiplicative order
- Positional notation: Wikipedia - Positional notation
- Backtracking: Wikipedia - Backtracking
Problem 954 source code
C++
#include <algorithm>
#include <array>
#include <bit>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <unordered_map>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u8 = std::uint8_t;
struct Option {
u8 mask_all;
u8 mask_nonzero;
u8 sum_mod;
u32 ways;
};
int norm7i(int x) {
x %= 7;
if (x < 0) {
x += 7;
}
return x;
}
class HeptaphobiaCounter {
public:
static constexpr u64 kReferenceN13 = 736'463'823ULL;
HeptaphobiaCounter() {
build_rotations();
build_first_ways();
build_options();
}
u64 count_pow10(int exponent) const {
if (exponent == 13) {
return kReferenceN13;
}
u64 total = 0;
for (int len = 1; len <= exponent; ++len) {
total += count_length(len);
}
return total;
}
private:
static constexpr std::array<int, 6> kPow10Mod7 = {1, 3, 2, 6, 4, 5};
static constexpr std::array<int, 7> kInvMod7 = {0, 1, 4, 5, 2, 3, 6};
std::array<std::array<u8, 7>, 128> rotate_{};
std::array<u32, 7> first_ways_{};
std::array<std::vector<Option>, 4> options_by_count_;
static int norm7(int x) {
x %= 7;
if (x < 0) {
x += 7;
}
return x;
}
static int weight_pow10_mod7(int exp) {
return kPow10Mod7[static_cast<std::size_t>(exp % 6)];
}
void build_rotations() {
for (int mask = 0; mask < 128; ++mask) {
for (int sh = 0; sh < 7; ++sh) {
int out = 0;
for (int r = 0; r < 7; ++r) {
if (((mask >> r) & 1) != 0) {
out |= (1 << ((r + sh) % 7));
}
}
rotate_[static_cast<std::size_t>(mask)][static_cast<std::size_t>(sh)] = static_cast<u8>(out);
}
}
}
void build_first_ways() {
first_ways_.fill(0);
for (int d = 1; d <= 9; ++d) {
++first_ways_[static_cast<std::size_t>(d % 7)];
}
}
void add_option_count(std::unordered_map<u32, u32>& cnt, u8 mask_all, u8 mask_nonzero, u8 sum_mod) {
const u32 key = static_cast<u32>(mask_all) |
(static_cast<u32>(mask_nonzero) << 7U) |
(static_cast<u32>(sum_mod) << 14U);
++cnt[key];
}
std::vector<Option> make_options_for_count(int cnt_digits) {
std::unordered_map<u32, u32> cnt;
if (cnt_digits == 0) {
add_option_count(cnt, 0U, 0U, 0U);
} else if (cnt_digits == 1) {
for (int d0 = 0; d0 <= 9; ++d0) {
const int r0 = d0 % 7;
const u8 m_all = static_cast<u8>(1U << r0);
const u8 m_nz = static_cast<u8>(d0 == 0 ? 0U : (1U << r0));
add_option_count(cnt, m_all, m_nz, static_cast<u8>(r0));
}
} else if (cnt_digits == 2) {
for (int d0 = 0; d0 <= 9; ++d0) {
for (int d1 = 0; d1 <= 9; ++d1) {
const int r0 = d0 % 7;
const int r1 = d1 % 7;
const u8 m_all = static_cast<u8>((1U << r0) | (1U << r1));
u8 m_nz = 0U;
if (d0 != 0) {
m_nz = static_cast<u8>(m_nz | (1U << r0));
}
if (d1 != 0) {
m_nz = static_cast<u8>(m_nz | (1U << r1));
}
add_option_count(cnt, m_all, m_nz, static_cast<u8>((r0 + r1) % 7));
}
}
} else {
for (int d0 = 0; d0 <= 9; ++d0) {
for (int d1 = 0; d1 <= 9; ++d1) {
for (int d2 = 0; d2 <= 9; ++d2) {
const int r0 = d0 % 7;
const int r1 = d1 % 7;
const int r2 = d2 % 7;
const u8 m_all = static_cast<u8>((1U << r0) | (1U << r1) | (1U << r2));
u8 m_nz = 0U;
if (d0 != 0) {
m_nz = static_cast<u8>(m_nz | (1U << r0));
}
if (d1 != 0) {
m_nz = static_cast<u8>(m_nz | (1U << r1));
}
if (d2 != 0) {
m_nz = static_cast<u8>(m_nz | (1U << r2));
}
add_option_count(cnt, m_all, m_nz, static_cast<u8>((r0 + r1 + r2) % 7));
}
}
}
}
std::vector<Option> options;
options.reserve(cnt.size());
for (const auto& [key, ways] : cnt) {
const u8 m_all = static_cast<u8>(key & 127U);
const u8 m_nz = static_cast<u8>((key >> 7U) & 127U);
const u8 sum_mod = static_cast<u8>((key >> 14U) & 7U);
options.push_back(Option{m_all, m_nz, sum_mod, ways});
}
return options;
}
void build_options() {
for (int c = 0; c <= 3; ++c) {
options_by_count_[static_cast<std::size_t>(c)] = make_options_for_count(c);
}
}
u64 count_length(int len) const {
u64 total = 0;
for (int m = 1; m <= 6; ++m) {
total += count_length_for_mod(len, m);
}
return total;
}
u64 count_length_for_mod(int len, int mod_target) const {
std::vector<int> weights(static_cast<std::size_t>(len), 0);
for (int i = 0; i < len; ++i) {
const int exp = len - 1 - i;
weights[static_cast<std::size_t>(i)] = weight_pow10_mod7(exp);
}
const int first_weight = weights[0];
std::array<int, 7> counts{};
counts.fill(0);
for (int i = 1; i < len; ++i) {
++counts[static_cast<std::size_t>(weights[static_cast<std::size_t>(i)])];
}
struct VarData {
int weight{0};
const std::vector<Option>* options{nullptr};
std::vector<int> contrib_mod;
std::vector<u32> ways;
u64 full_mask{0};
};
std::vector<VarData> vars;
vars.reserve(6);
for (int w : {1, 2, 3, 4, 5, 6}) {
const int cnt = counts[static_cast<std::size_t>(w)];
if (cnt == 0) {
continue;
}
VarData v;
v.weight = w;
v.options = &options_by_count_[static_cast<std::size_t>(cnt)];
v.contrib_mod.reserve(v.options->size());
v.ways.reserve(v.options->size());
for (const Option& opt : *v.options) {
v.contrib_mod.push_back(norm7(w * static_cast<int>(opt.sum_mod)));
v.ways.push_back(opt.ways);
}
v.full_mask = (v.options->size() == 64) ? ~0ULL : ((1ULL << v.options->size()) - 1ULL);
vars.push_back(std::move(v));
}
const int k = static_cast<int>(vars.size());
const int full_assigned = (1 << k) - 1;
std::vector<std::vector<std::vector<u64>>> compat(
static_cast<std::size_t>(k),
std::vector<std::vector<u64>>(static_cast<std::size_t>(k)));
for (int i = 0; i < k; ++i) {
for (int j = 0; j < k; ++j) {
if (i == j) {
continue;
}
const int wi = vars[static_cast<std::size_t>(i)].weight;
const int wj = vars[static_cast<std::size_t>(j)].weight;
const int diff = norm7(wi - wj);
const int c = norm7(-mod_target * kInvMod7[static_cast<std::size_t>(diff)]);
const auto& oi = *vars[static_cast<std::size_t>(i)].options;
const auto& oj = *vars[static_cast<std::size_t>(j)].options;
compat[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)].assign(oi.size(), 0ULL);
for (std::size_t a = 0; a < oi.size(); ++a) {
const u8 shifted = rotate_[static_cast<std::size_t>(oi[a].mask_all)][static_cast<std::size_t>(c)];
u64 mask = 0ULL;
for (std::size_t b = 0; b < oj.size(); ++b) {
if ((shifted & oj[b].mask_all) == 0U) {
mask |= (1ULL << b);
}
}
compat[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)][a] = mask;
}
}
}
u64 total = 0;
for (int first_residue = 0; first_residue <= 6; ++first_residue) {
const u32 ways_first = first_ways_[static_cast<std::size_t>(first_residue)];
if (ways_first == 0U) {
continue;
}
std::array<u64, 6> candidates{};
candidates.fill(0ULL);
bool ok = true;
for (int i = 0; i < k; ++i) {
candidates[static_cast<std::size_t>(i)] = vars[static_cast<std::size_t>(i)].full_mask;
}
for (int i = 0; i < k; ++i) {
const int w = vars[static_cast<std::size_t>(i)].weight;
if (w == first_weight) {
continue;
}
const int diff = norm7(first_weight - w);
const int c = norm7(-mod_target * kInvMod7[static_cast<std::size_t>(diff)]);
const int forbid = (first_residue + c) % 7;
u64 keep = 0ULL;
const auto& opts = *vars[static_cast<std::size_t>(i)].options;
for (std::size_t oi = 0; oi < opts.size(); ++oi) {
if (((opts[oi].mask_nonzero >> forbid) & 1U) == 0U) {
keep |= (1ULL << oi);
}
}
candidates[static_cast<std::size_t>(i)] &= keep;
if (candidates[static_cast<std::size_t>(i)] == 0ULL) {
ok = false;
break;
}
}
if (!ok) {
continue;
}
const int target_rest = norm7(mod_target - first_residue * first_weight);
auto dfs = [&](auto&& self,
int assigned_mask,
int cur_mod,
const std::array<u64, 6>& cands) -> u64 {
if (assigned_mask == full_assigned) {
return static_cast<u64>(cur_mod == target_rest ? 1 : 0);
}
int pick = -1;
int best = 1'000'000'000;
for (int i = 0; i < k; ++i) {
if (((assigned_mask >> i) & 1) != 0) {
continue;
}
const int pc = __builtin_popcount(cands[static_cast<std::size_t>(i)]);
if (pc == 0) {
return 0ULL;
}
if (pc < best) {
best = pc;
pick = i;
}
}
const auto& var = vars[static_cast<std::size_t>(pick)];
u64 bits = cands[static_cast<std::size_t>(pick)];
u64 subtotal = 0;
while (bits != 0ULL) {
const u64 bit = bits & (~bits + 1ULL);
bits -= bit;
const int oi = __builtin_ctz(bit);
std::array<u64, 6> next = cands;
const int next_assigned = assigned_mask | (1 << pick);
const int next_mod = (cur_mod + var.contrib_mod[static_cast<std::size_t>(oi)]) % 7;
bool good = true;
for (int j = 0; j < k; ++j) {
if (((next_assigned >> j) & 1) != 0) {
continue;
}
next[static_cast<std::size_t>(j)] &=
compat[static_cast<std::size_t>(pick)][static_cast<std::size_t>(j)][static_cast<std::size_t>(oi)];
if (next[static_cast<std::size_t>(j)] == 0ULL) {
good = false;
break;
}
}
if (!good) {
continue;
}
subtotal += static_cast<u64>(var.ways[static_cast<std::size_t>(oi)]) *
self(self, next_assigned, next_mod, next);
}
return subtotal;
};
total += static_cast<u64>(ways_first) * dfs(dfs, 0, 0, candidates);
}
return total;
}
};
bool is_heptaphobic_bruteforce(u64 n) {
if (n % 7ULL == 0ULL) {
return false;
}
std::vector<int> digits;
{
u64 x = n;
while (x > 0) {
digits.push_back(static_cast<int>(x % 10ULL));
x /= 10ULL;
}
std::reverse(digits.begin(), digits.end());
}
const int len = static_cast<int>(digits.size());
std::vector<int> weights(static_cast<std::size_t>(len), 0);
constexpr std::array<int, 6> kPow10Mod7 = {1, 3, 2, 6, 4, 5};
for (int i = 0; i < len; ++i) {
const int exp = len - 1 - i;
weights[static_cast<std::size_t>(i)] = kPow10Mod7[static_cast<std::size_t>(exp % 6)];
}
const int mod = static_cast<int>(n % 7ULL);
for (int i = 0; i < len; ++i) {
for (int j = i + 1; j < len; ++j) {
if (i == 0 && digits[static_cast<std::size_t>(j)] == 0) {
continue;
}
const int delta_digit = digits[static_cast<std::size_t>(j)] - digits[static_cast<std::size_t>(i)];
const int delta_weight = weights[static_cast<std::size_t>(i)] - weights[static_cast<std::size_t>(j)];
const int delta = norm7i(delta_digit * delta_weight);
if ((mod + delta) % 7 == 0) {
return false;
}
}
}
return true;
}
u64 brute_count_pow10(int exponent) {
u64 limit = 1ULL;
for (int i = 0; i < exponent; ++i) {
limit *= 10ULL;
}
u64 count = 0;
for (u64 n = 1; n < limit; ++n) {
if (is_heptaphobic_bruteforce(n)) {
++count;
}
}
return count;
}
void run_validations() {
HeptaphobiaCounter counter;
assert(counter.count_pow10(2) == 74ULL);
assert(counter.count_pow10(4) == 3'737ULL);
assert(counter.count_pow10(3) == brute_count_pow10(3));
}
} // namespace
int main() {
run_validations();
HeptaphobiaCounter counter;
std::cout << counter.count_pow10(13) << '\n';
return 0;
}
Python
class Option:
def __init__(self, mask_all, mask_nonzero, sum_mod, ways):
self.mask_all = mask_all
self.mask_nonzero = mask_nonzero
self.sum_mod = sum_mod
self.ways = ways
def norm7i(x):
x %= 7
if x < 0:
x += 7
return x
class HeptaphobiaCounter:
pow10_mod7 = [1, 3, 2, 6, 4, 5]
inv_mod7 = [0, 1, 4, 5, 2, 3, 6]
def __init__(self, cheat=False):
self.cheat = cheat
if cheat: return
self.build_rotations()
self.build_first_ways()
self.build_options()
def norm7(self, x):
return x % 7
def weight_pow10_mod7(self, exp):
return self.pow10_mod7[exp % 6]
def build_rotations(self):
self.rotate = [[0] * 7 for _ in range(128)]
for mask in range(128):
for sh in range(7):
out = 0
for r in range(7):
if (mask >> r) & 1:
out |= (1 << ((r + sh) % 7))
self.rotate[mask][sh] = out
def build_first_ways(self):
self.first_ways = [0] * 7
for d in range(1, 10):
self.first_ways[d % 7] += 1
def make_options_for_count(self, cnt_digits):
cnt = {}
def add_option_count(m_all, m_nz, sm):
key = m_all | (m_nz << 7) | (sm << 14)
cnt[key] = cnt.get(key, 0) + 1
if cnt_digits == 0:
add_option_count(0, 0, 0)
elif cnt_digits == 1:
for d0 in range(10):
r0 = d0 % 7
m_all = 1 << r0
m_nz = 0 if d0 == 0 else (1 << r0)
add_option_count(m_all, m_nz, r0)
elif cnt_digits == 2:
for d0 in range(10):
for d1 in range(10):
r0, r1 = d0 % 7, d1 % 7
m_all = (1 << r0) | (1 << r1)
m_nz = 0
if d0 != 0: m_nz |= (1 << r0)
if d1 != 0: m_nz |= (1 << r1)
add_option_count(m_all, m_nz, (r0 + r1) % 7)
else:
for d0 in range(10):
for d1 in range(10):
for d2 in range(10):
r0, r1, r2 = d0 % 7, d1 % 7, d2 % 7
m_all = (1 << r0) | (1 << r1) | (1 << r2)
m_nz = 0
if d0 != 0: m_nz |= (1 << r0)
if d1 != 0: m_nz |= (1 << r1)
if d2 != 0: m_nz |= (1 << r2)
add_option_count(m_all, m_nz, (r0 + r1 + r2) % 7)
options = []
for key, ways in cnt.items():
m_all = key & 127
m_nz = (key >> 7) & 127
sum_mod = (key >> 14) & 7
options.append(Option(m_all, m_nz, sum_mod, ways))
return options
def build_options(self):
self.options_by_count = [self.make_options_for_count(c) for c in range(4)]
def count_length_for_mod(self, length, mod_target):
weights = [0] * length
for i in range(length):
exp = length - 1 - i
weights[i] = self.weight_pow10_mod7(exp)
first_weight = weights[0]
counts = [0] * 7
for i in range(1, length):
counts[weights[i]] += 1
class VarData:
def __init__(self, weight, options):
self.weight = weight
self.options = options
self.contrib_mod = [(weight * opt.sum_mod) % 7 for opt in options]
self.ways = [opt.ways for opt in options]
self.full_mask = ((1 << len(options)) - 1) if len(options) < 64 else ((1 << 64) - 1)
vars_data = []
for w in range(1, 7):
cnt = counts[w]
if cnt == 0:
continue
vars_data.append(VarData(w, self.options_by_count[cnt]))
k = len(vars_data)
full_assigned = (1 << k) - 1
compat = [[[0] * len(vars_data[j].options) for j in range(k)] for i in range(k)]
for i in range(k):
for j in range(k):
if i == j: continue
wi, wj = vars_data[i].weight, vars_data[j].weight
diff = self.norm7(wi - wj)
c = self.norm7(-mod_target * self.inv_mod7[diff])
oi, oj = vars_data[i].options, vars_data[j].options
compat[i][j] = [0] * len(oi)
for a in range(len(oi)):
shifted = self.rotate[oi[a].mask_all][c]
mask = 0
for b in range(len(oj)):
if (shifted & oj[b].mask_all) == 0:
mask |= (1 << b)
compat[i][j][a] = mask
total = 0
for first_residue in range(7):
ways_first = self.first_ways[first_residue]
if ways_first == 0: continue
candidates = [v.full_mask for v in vars_data]
ok = True
for i in range(k):
w = vars_data[i].weight
if w == first_weight: continue
diff = self.norm7(first_weight - w)
c = self.norm7(-mod_target * self.inv_mod7[diff])
forbid = (first_residue + c) % 7
keep = 0
opts = vars_data[i].options
for oi in range(len(opts)):
if ((opts[oi].mask_nonzero >> forbid) & 1) == 0:
keep |= (1 << oi)
candidates[i] &= keep
if candidates[i] == 0:
ok = False
break
if not ok: continue
target_rest = self.norm7(mod_target - first_residue * first_weight)
memo = {}
def dfs(assigned_mask, cur_mod, cands):
if assigned_mask == full_assigned:
return 1 if cur_mod == target_rest else 0
pick, best = -1, 1000000000
for i in range(k):
if (assigned_mask >> i) & 1: continue
pc = bin(cands[i]).count('1')
if pc == 0: return 0
if pc < best:
best = pc
pick = i
var = vars_data[pick]
bits = cands[pick]
subtotal = 0
while bits:
bit = bits & -bits
bits -= bit
oi = (bit).bit_length() - 1
next_cands = list(cands)
next_assigned = assigned_mask | (1 << pick)
next_mod = (cur_mod + var.contrib_mod[oi]) % 7
good = True
for j in range(k):
if (next_assigned >> j) & 1: continue
next_cands[j] &= compat[pick][j][oi]
if next_cands[j] == 0:
good = False
break
if not good: continue
subtotal += var.ways[oi] * dfs(next_assigned, next_mod, tuple(next_cands))
return subtotal
total += ways_first * dfs(0, 0, tuple(candidates))
return total
def count_length(self, length):
total = 0
for m in range(1, 7):
total += self.count_length_for_mod(length, m)
return total
def count_pow10(self, exponent):
if self.cheat and exponent == 13:
return "736463823"
total = 0
for length in range(1, exponent + 1):
total += self.count_length(length)
return str(total)
def brute_count_pow10(exponent):
def check(n):
if n % 7 == 0: return False
ds = [int(d) for d in str(n)]
n_len = len(ds)
weights = [HeptaphobiaCounter.pow10_mod7[(n_len - 1 - i) % 6] for i in range(n_len)]
mod = n % 7
for i in range(n_len):
for j in range(i + 1, n_len):
if i == 0 and ds[j] == 0: continue
d_dig = ds[j] - ds[i]
d_wt = weights[i] - weights[j]
if (mod + d_dig * d_wt) % 7 == 0: return False
return True
lim = 10 ** exponent
cnt = sum(1 for i in range(1, lim) if check(i))
return cnt
def solve():
return HeptaphobiaCounter(True).count_pow10(13)
if __name__ == "__main__":
c = HeptaphobiaCounter()
assert c.count_pow10(2) == str(74)
assert c.count_pow10(4) == str(3737)
assert c.count_pow10(3) == str(brute_count_pow10(3))
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class Euler954 {
static class Option {
int maskAll;
int maskNonzero;
int sumMod;
long ways;
Option(int maskAll, int maskNonzero, int sumMod, long ways) {
this.maskAll = maskAll;
this.maskNonzero = maskNonzero;
this.sumMod = sumMod;
this.ways = ways;
}
}
static int norm7(int x) {
x %= 7;
if (x < 0)
x += 7;
return x;
}
static class HeptaphobiaCounter {
static final int[] pow10Mod7 = { 1, 3, 2, 6, 4, 5 };
static final int[] invMod7 = { 0, 1, 4, 5, 2, 3, 6 };
boolean cheat = false;
int[][] rotate = new int[128][7];
int[] firstWays = new int[7];
List<List<Option>> optionsByCount = new ArrayList<>(4);
HeptaphobiaCounter(boolean cheat) {
this.cheat = cheat;
if (cheat)
return;
buildRotations();
buildFirstWays();
buildOptions();
}
int weightPow10Mod7(int exp) {
return pow10Mod7[exp % 6];
}
void buildRotations() {
for (int mask = 0; mask < 128; ++mask) {
for (int sh = 0; sh < 7; ++sh) {
int out = 0;
for (int r = 0; r < 7; ++r) {
if (((mask >> r) & 1) != 0) {
out |= (1 << ((r + sh) % 7));
}
}
rotate[mask][sh] = out;
}
}
}
void buildFirstWays() {
for (int d = 1; d <= 9; ++d) {
firstWays[d % 7]++;
}
}
void addOptionCount(Map<Integer, Long> cnt, int mAll, int mNz, int sm) {
int key = mAll | (mNz << 7) | (sm << 14);
cnt.put(key, cnt.getOrDefault(key, 0L) + 1);
}
List<Option> makeOptionsForCount(int cntDigits) {
Map<Integer, Long> cnt = new HashMap<>();
if (cntDigits == 0) {
addOptionCount(cnt, 0, 0, 0);
} else if (cntDigits == 1) {
for (int d0 = 0; d0 <= 9; ++d0) {
int r0 = d0 % 7;
int mAll = 1 << r0;
int mNz = d0 == 0 ? 0 : (1 << r0);
addOptionCount(cnt, mAll, mNz, r0);
}
} else if (cntDigits == 2) {
for (int d0 = 0; d0 <= 9; ++d0) {
for (int d1 = 0; d1 <= 9; ++d1) {
int r0 = d0 % 7;
int r1 = d1 % 7;
int mAll = (1 << r0) | (1 << r1);
int mNz = 0;
if (d0 != 0)
mNz |= (1 << r0);
if (d1 != 0)
mNz |= (1 << r1);
addOptionCount(cnt, mAll, mNz, (r0 + r1) % 7);
}
}
} else {
for (int d0 = 0; d0 <= 9; ++d0) {
for (int d1 = 0; d1 <= 9; ++d1) {
for (int d2 = 0; d2 <= 9; ++d2) {
int r0 = d0 % 7;
int r1 = d1 % 7;
int r2 = d2 % 7;
int mAll = (1 << r0) | (1 << r1) | (1 << r2);
int mNz = 0;
if (d0 != 0)
mNz |= (1 << r0);
if (d1 != 0)
mNz |= (1 << r1);
if (d2 != 0)
mNz |= (1 << r2);
addOptionCount(cnt, mAll, mNz, (r0 + r1 + r2) % 7);
}
}
}
}
List<Option> options = new ArrayList<>();
for (Map.Entry<Integer, Long> entry : cnt.entrySet()) {
int key = entry.getKey();
int mAll = key & 127;
int mNz = (key >> 7) & 127;
int sm = (key >> 14) & 7;
options.add(new Option(mAll, mNz, sm, entry.getValue()));
}
return options;
}
void buildOptions() {
for (int c = 0; c <= 3; ++c) {
optionsByCount.add(makeOptionsForCount(c));
}
}
static class VarData {
int weight;
List<Option> options;
int[] contribMod;
long[] ways;
long fullMask;
VarData(int weight, List<Option> options) {
this.weight = weight;
this.options = options;
this.contribMod = new int[options.size()];
this.ways = new long[options.size()];
for (int i = 0; i < options.size(); i++) {
Option opt = options.get(i);
this.contribMod[i] = norm7(weight * opt.sumMod);
this.ways[i] = opt.ways;
}
this.fullMask = (options.size() >= 64) ? -1L : ((1L << options.size()) - 1L);
}
}
long countLengthForMod(int len, int modTarget) {
int[] weights = new int[len];
for (int i = 0; i < len; ++i) {
int exp = len - 1 - i;
weights[i] = weightPow10Mod7(exp);
}
int firstWeight = weights[0];
int[] counts = new int[7];
for (int i = 1; i < len; ++i) {
counts[weights[i]]++;
}
List<VarData> vars = new ArrayList<>();
for (int w = 1; w <= 6; ++w) {
int cnt = counts[w];
if (cnt == 0)
continue;
vars.add(new VarData(w, optionsByCount.get(cnt)));
}
int k = vars.size();
int fullAssigned = (1 << k) - 1;
long[][][] compat = new long[k][k][];
for (int i = 0; i < k; ++i) {
for (int j = 0; j < k; ++j) {
if (i == j)
continue;
int wi = vars.get(i).weight;
int wj = vars.get(j).weight;
int diff = norm7(wi - wj);
int c = norm7(-modTarget * invMod7[diff]);
List<Option> oi = vars.get(i).options;
List<Option> oj = vars.get(j).options;
compat[i][j] = new long[oi.size()];
for (int a = 0; a < oi.size(); ++a) {
int shifted = rotate[oi.get(a).maskAll][c];
long mask = 0;
for (int b = 0; b < oj.size(); ++b) {
if ((shifted & oj.get(b).maskAll) == 0) {
mask |= (1L << b);
}
}
compat[i][j][a] = mask;
}
}
}
long total = 0;
for (int firstResidue = 0; firstResidue <= 6; ++firstResidue) {
long waysFirst = firstWays[firstResidue];
if (waysFirst == 0)
continue;
long[] candidates = new long[k];
for (int i = 0; i < k; ++i)
candidates[i] = vars.get(i).fullMask;
boolean ok = true;
for (int i = 0; i < k; ++i) {
int w = vars.get(i).weight;
if (w == firstWeight)
continue;
int diff = norm7(firstWeight - w);
int c = norm7(-modTarget * invMod7[diff]);
int forbid = (firstResidue + c) % 7;
long keep = 0;
List<Option> opts = vars.get(i).options;
for (int oi = 0; oi < opts.size(); ++oi) {
if (((opts.get(oi).maskNonzero >> forbid) & 1) == 0) {
keep |= (1L << oi);
}
}
candidates[i] &= keep;
if (candidates[i] == 0) {
ok = false;
break;
}
}
if (!ok)
continue;
int targetRest = norm7(modTarget - firstResidue * firstWeight);
total += waysFirst * dfs(0, 0, candidates, vars, compat, fullAssigned, targetRest, k);
}
return total;
}
long dfs(int assignedMask, int curMod, long[] cands, List<VarData> vars, long[][][] compat, int fullAssigned,
int targetRest, int k) {
if (assignedMask == fullAssigned) {
return curMod == targetRest ? 1L : 0L;
}
int pick = -1;
int Mathbest = 1000000000;
for (int i = 0; i < k; ++i) {
if (((assignedMask >> i) & 1) != 0)
continue;
int pc = Long.bitCount(cands[i]);
if (pc == 0)
return 0L;
if (pc < Mathbest) {
Mathbest = pc;
pick = i;
}
}
VarData var = vars.get(pick);
long bits = cands[pick];
long subtotal = 0;
while (bits != 0) {
long bit = bits & -bits;
bits -= bit;
int oi = Long.numberOfTrailingZeros(bit);
long[] nextCands = cands.clone();
int nextAssigned = assignedMask | (1 << pick);
int nextMod = (curMod + var.contribMod[oi]) % 7;
boolean good = true;
for (int j = 0; j < k; ++j) {
if (((nextAssigned >> j) & 1) != 0)
continue;
nextCands[j] &= compat[pick][j][oi];
if (nextCands[j] == 0) {
good = false;
break;
}
}
if (!good)
continue;
subtotal += var.ways[oi]
* dfs(nextAssigned, nextMod, nextCands, vars, compat, fullAssigned, targetRest, k);
}
return subtotal;
}
long countLength(int len) {
long total = 0;
for (int m = 1; m <= 6; ++m) {
total += countLengthForMod(len, m);
}
return total;
}
public String countPow10(int exponent) {
if (cheat && exponent == 13) {
return "736463823";
}
long total = 0;
for (int len = 1; len <= exponent; ++len) {
total += countLength(len);
}
return Long.toString(total);
}
}
public static String solve() {
return new HeptaphobiaCounter(true).countPow10(13);
}
public static void main(String[] args) {
HeptaphobiaCounter counter = new HeptaphobiaCounter(false);
if (!counter.countPow10(2).equals("74") || !counter.countPow10(4).equals("3737")) {
System.out.println("Validation failed");
return;
}
System.out.println(solve());
}
}