Problem 886: Coprime Permutations
View on Project EulerProject Euler Problem 886 Solution
EulerSolve provides an optimized solution for Project Euler Problem 886, Coprime Permutations, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For even \(n\), let \(P(n)\) be the number of permutations of \(2,3,\dots,n\) in which every adjacent pair is coprime. The target computation is \(P(34)\) modulo \(83456729\). A direct search grows far too quickly, so the implementations convert the problem into an exact count of Hamiltonian cycles in a balanced bipartite graph and then evaluate that count with a determinant/permanent inclusion-exclusion identity. Mathematical Approach Write \(n=2m\). Introduce the odd side $$U=\{1,3,5,\dots,2m-1\}$$ and the even side $$V=\{2,4,6,\dots,2m\}.$$ Define the \(m\times m\) bipartite adjacency matrix \(B\) by $$B_{u,v}=\begin{cases} 1,& \gcd(u,v)=1,\\ 0,& \text{otherwise}. \end{cases}$$ The key point is that the original permutations of \(2,\dots,2m\) can be re-expressed as alternating structures in this augmented graph. Step 1: Parity Forces an Alternating Path Among the numbers \(2,3,\dots,2m\), there are \(m\) even numbers and \(m-1\) odd numbers. Two even numbers can never be adjacent in a valid permutation, because their gcd is at least \(2\). Therefore every valid permutation must alternate parity. Since the evens are more numerous by one, the parity pattern is forced to be $$E-O-E-O-\dots-O-E.$$ So the problem is not counting arbitrary coprime permutations; it is counting alternating even-odd paths whose endpoints are both even....
Detailed mathematical approach
Problem Summary
For even \(n\), let \(P(n)\) be the number of permutations of \(2,3,\dots,n\) in which every adjacent pair is coprime. The target computation is \(P(34)\) modulo \(83456729\). A direct search grows far too quickly, so the implementations convert the problem into an exact count of Hamiltonian cycles in a balanced bipartite graph and then evaluate that count with a determinant/permanent inclusion-exclusion identity.
Mathematical Approach
Write \(n=2m\). Introduce the odd side
$$U=\{1,3,5,\dots,2m-1\}$$
and the even side
$$V=\{2,4,6,\dots,2m\}.$$
Define the \(m\times m\) bipartite adjacency matrix \(B\) by
$$B_{u,v}=\begin{cases} 1,& \gcd(u,v)=1,\\ 0,& \text{otherwise}. \end{cases}$$
The key point is that the original permutations of \(2,\dots,2m\) can be re-expressed as alternating structures in this augmented graph.
Step 1: Parity Forces an Alternating Path
Among the numbers \(2,3,\dots,2m\), there are \(m\) even numbers and \(m-1\) odd numbers. Two even numbers can never be adjacent in a valid permutation, because their gcd is at least \(2\). Therefore every valid permutation must alternate parity. Since the evens are more numerous by one, the parity pattern is forced to be
$$E-O-E-O-\dots-O-E.$$
So the problem is not counting arbitrary coprime permutations; it is counting alternating even-odd paths whose endpoints are both even.
Step 2: Add \(1\) and Turn the Path into a Cycle
The number \(1\) is coprime to every even number, so adding it on the odd side balances the bipartition. If
$$e_0,o_1,e_1,o_2,\dots,o_{m-1},e_{m-1}$$
is a valid permutation of \(2,\dots,2m\), then inserting \(1\) between the two endpoints produces the alternating cycle
$$1,e_0,o_1,e_1,\dots,o_{m-1},e_{m-1},1.$$
Conversely, deleting \(1\) from such a cycle recovers a valid permutation. Thus \(P(2m)\) is exactly the number of Hamiltonian cycles in the balanced bipartite graph \((U,V)\), read from the distinguished vertex \(1\).
Step 3: Hamiltonian Cycles Become Ordered Pairs of Perfect Matchings
Any alternating cycle in a bipartite graph can be colored with two alternating edge colors. That splits the cycle into two perfect matchings, say \(M_1\) and \(M_2\). Conversely, the union of two perfect matchings is a spanning \(2\)-regular bipartite subgraph, so it is a disjoint union of even cycles.
A single Hamiltonian cycle is therefore the special case in which the ordered pair \((M_1,M_2)\) produces exactly one cycle rather than several disconnected cycles. The permanent counts perfect matchings without signs, while the determinant supplies the signs needed to cancel the unwanted multi-cycle decompositions.
Step 4: The Inclusion-Exclusion Identity Used by the Implementations
Let the distinguished odd vertex be \(1\). For subsets \(I\subseteq U\setminus\{1\}\) and \(J\subseteq V\) with \(|I|=|J|=k\), denote by \(B_{I,J}\) the selected \(k\times k\) submatrix and by \(B_{\bar I,\bar J}\) the complementary \((m-k)\times(m-k)\) submatrix, where \(\bar I=U\setminus I\) still contains the distinguished vertex \(1\). The exact identity evaluated by the code is
$$P(2m)=\sum_{k=0}^{m-1}\ \sum_{\substack{I\subseteq U\setminus\{1\}\\J\subseteq V\\|I|=|J|=k}} (-1)^k \det(B_{I,J})^2\operatorname{perm}(B_{\bar I,\bar J})^2 \pmod{83456729}.$$
The squared permanent counts ordered pairs of perfect matchings on the complementary block. The squared determinant creates the sign-reversing correction that removes disconnected cycle covers and leaves only Hamiltonian cycles through the distinguished row.
Step 5: Compress Equal Neighborhoods
If two selectable odd vertices have identical neighborhoods, then choosing both of them inside a determinant block would create two equal rows, so the determinant would vanish. The same is true for equal columns. That means the determinant part only needs one representative from each row class and each column class, with a multiplicity factor recording how many interchangeable choices exist.
The permanent part is different: repeated columns do not vanish there. Instead, the implementations use a grouped form of Ryser's formula. If a column class has size \(s\) and exactly \(t\) columns from that class are selected in a Ryser subset, then there are
$$\binom{s}{t}$$
ways to make that choice, and the row sums depend only on \(t\), not on which particular columns were taken.
Worked Example: \(n=4\)
Here \(m=2\), so
$$U=\{1,3\},\qquad V=\{2,4\}.$$
Every odd-even pair is coprime, hence
$$B=\begin{pmatrix} 1 & 1\\ 1 & 1 \end{pmatrix}.$$
After anchoring the count at \(1\), only one non-distinguished odd row remains selectable. The formula has two types of contribution:
$$k=0:\qquad \operatorname{perm}(B)^2=2^2=4,$$
and
$$k=1:\qquad -\sum_{|J|=1}\det(B_{\{3\},J})^2\operatorname{perm}(B_{\{1\},V\setminus J})^2=-(1+1)=-2.$$
Therefore
$$P(4)=4-2=2,$$
corresponding to the two valid permutations \((2,3,4)\) and \((4,3,2)\). This is the smallest checkpoint used by the implementations.
How the Code Works
The C++, Python, and Java implementations first build the \(m\times m\) coprimality matrix between \(U\) and \(V\), together with binomial coefficients up to \(m\). They then classify rows by equal support masks, remove one copy of the distinguished row \(1\), and classify columns in the same way by applying the same mask logic to the transpose.
For every \(k\), the implementation enumerates all choices of \(k\) row classes and \(k\) column classes. The determinant block is formed from the chosen representatives; if its determinant is \(0\), that branch is discarded immediately. Otherwise the code squares the determinant, applies the sign \((-1)^k\), and multiplies by the class multiplicities.
The complementary block is then sent to a grouped permanent routine. The determinant is computed with fraction-free elimination modulo \(83456729\). The permanent is computed with Ryser's formula, but columns with identical support masks are grouped: a small prefix of classes is handled by explicit bitmask subsets, and the remaining large classes are handled recursively by choosing only how many columns are selected from each class. Each choice contributes a binomial factor \(\binom{s}{t}\).
Finally the code multiplies
$$\text{multiplicity}\times (-1)^k \times \det^2 \times \operatorname{perm}^2$$
for each branch and accumulates the result modulo \(83456729\). The implementations include small checkpoints such as \(P(4)=2\) and \(P(10)=576\), and the C++ version also cross-checks tiny cases by brute force.
Complexity Analysis
Let \(m=n/2\). If \(R\) and \(C\) are the numbers of distinct selectable row and column classes, the outer enumeration examines
$$\sum_{k=0}^{m-1}\binom{R}{k}\binom{C}{k}$$
class pairs. A determinant on a \(k\times k\) block costs \(O(k^3)\). For the complementary block of size \(r=m-k\), naive Ryser would cost \(O(r2^r)\), but grouping identical columns reduces the subset space to roughly
$$2^b\prod_i (s_i+1),$$
where \(b<10\) is the explicitly enumerated prefix size and the \(s_i\) are the remaining column-class sizes. The permanent stage then spends about \(O\!\left(r^2 2^b\prod_i (s_i+1)\right)\) arithmetic operations. Memory usage is polynomial in the current block size, plus the explicit subset table of size \(2^b\). For the target case \(n=34\), this compression is what makes the exact count practical.
Footnotes and References
- Project Euler problem page: https://projecteuler.net/problem=886
- Permanent of a matrix: Wikipedia - Permanent (mathematics)
- Determinant: Wikipedia - Determinant
- Hamiltonian path and cycle: Wikipedia - Hamiltonian path
- Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
- Ryser formula: Wikipedia - Ryser formula
Problem 886 source code
C++
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <functional>
#include <map>
#include <numeric>
#include <vector>
namespace {
using i64 = long long;
using u64 = std::uint64_t;
constexpr int MOD = 83'456'729;
int add_mod(int a, int b) {
int s = a + b;
if (s >= MOD) s -= MOD;
return s;
}
int sub_mod(int a, int b) {
int s = a - b;
if (s < 0) s += MOD;
return s;
}
int mul_mod(i64 a, i64 b) {
return static_cast<int>((static_cast<__int128>(a) * b) % MOD);
}
int mod_pow(int a, i64 e) {
i64 r = 1;
i64 x = a % MOD;
if (x < 0) x += MOD;
while (e > 0) {
if (e & 1LL) r = (r * x) % MOD;
x = (x * x) % MOD;
e >>= 1LL;
}
return static_cast<int>(r);
}
int mod_inv(int a) {
return mod_pow(a, MOD - 2);
}
std::vector<std::vector<int>> build_binom(int n) {
std::vector<std::vector<int>> C(static_cast<std::size_t>(n + 1));
C[0] = {1};
for (int i = 1; i <= n; ++i) {
C[static_cast<std::size_t>(i)].assign(static_cast<std::size_t>(i + 1), 1);
for (int j = 1; j < i; ++j) {
C[static_cast<std::size_t>(i)][static_cast<std::size_t>(j)] =
C[static_cast<std::size_t>(i - 1)][static_cast<std::size_t>(j - 1)] +
C[static_cast<std::size_t>(i - 1)][static_cast<std::size_t>(j)];
}
}
return C;
}
int determinant_mod(std::vector<std::vector<int>> A) {
const int n = static_cast<int>(A.size());
if (n == 0) return 1;
if (n == 1) return A[0][0] % MOD;
int result_inverse = 1;
for (int i = 0; i < n - 2; ++i) {
int row_to_use = -1;
for (int r = i; r < n; ++r) {
if (A[r][i] % MOD != 0) {
row_to_use = r;
break;
}
}
if (row_to_use < 0) return 0;
if (i != row_to_use) {
std::swap(A[row_to_use], A[i]);
result_inverse = (result_inverse == 0 ? 0 : MOD - result_inverse);
}
const int Aii = A[i][i] % MOD;
result_inverse = mul_mod(result_inverse, mod_pow(Aii, n - i - 2));
for (int r = i + 1; r < n; ++r) {
const int x = A[r][i] % MOD;
A[r][i] = 0;
for (int c = i + 1; c < n; ++c) {
int v = sub_mod(mul_mod(Aii, A[r][c]), mul_mod(x, A[i][c]));
A[r][c] = v;
}
}
}
int last = sub_mod(mul_mod(A[n - 2][n - 2], A[n - 1][n - 1]),
mul_mod(A[n - 2][n - 1], A[n - 1][n - 2]));
return mul_mod(last, mod_inv(result_inverse));
}
int permanent_mod_with_column_classes(const std::vector<std::vector<int>>& A,
const std::vector<std::vector<int>>& binom) {
const int n = static_cast<int>(A.size());
if (n == 0) return 1;
std::map<u64, std::vector<int>> column_classes;
for (int col = 0; col < n; ++col) {
u64 mask = 0;
for (int r = 0; r < n; ++r) {
if (A[r][col] != 0) mask |= (1ULL << r);
}
column_classes[mask].push_back(col);
}
std::vector<std::vector<int>> classes;
for (const auto& [k, inds] : column_classes) {
(void)k;
if (static_cast<int>(inds.size()) >= 3) {
classes.push_back(inds);
} else {
for (int idx : inds) {
classes.push_back({idx});
}
}
}
std::sort(classes.begin(), classes.end(), [](const auto& a, const auto& b) {
return a.size() < b.size();
});
std::vector<int> new_col_order;
new_col_order.reserve(n);
for (const auto& c : classes) {
for (int idx : c) new_col_order.push_back(idx);
}
std::vector<std::vector<int>> B(n, std::vector<int>(n));
for (int r = 0; r < n; ++r) {
for (int c = 0; c < n; ++c) {
B[r][c] = A[r][new_col_order[c]];
}
}
std::vector<int> class_sizes;
class_sizes.reserve(classes.size());
for (const auto& c : classes) class_sizes.push_back(static_cast<int>(c.size()));
const int default_block_exponent = 10;
int smallest_class_count = 0;
int smallest_class_index_count = 0;
while (smallest_class_count < static_cast<int>(class_sizes.size()) &&
smallest_class_index_count + class_sizes[smallest_class_count] < default_block_exponent) {
smallest_class_index_count += class_sizes[smallest_class_count];
++smallest_class_count;
}
int block_exponent = smallest_class_index_count;
int block_size = 1 << block_exponent;
std::vector<int> high_class_sizes;
std::vector<int> high_class_offsets;
int off = block_exponent;
for (int i = smallest_class_count; i < static_cast<int>(class_sizes.size()); ++i) {
high_class_sizes.push_back(class_sizes[i]);
high_class_offsets.push_back(off);
off += class_sizes[i];
}
std::vector<int> result_block(static_cast<std::size_t>(block_size), 0);
std::vector<int> high_mask(n, 0);
std::vector<int> popcnt(static_cast<std::size_t>(block_size), 0);
for (int s = 1; s < block_size; ++s) popcnt[static_cast<std::size_t>(s)] = popcnt[s >> 1] + (s & 1);
std::function<void(int, int, int)> rec = [&](int cls_idx, int total_ways, int high_ones) {
if (cls_idx == static_cast<int>(high_class_sizes.size())) {
for (int state = 0; state < block_size; ++state) {
int term = total_ways;
for (int r = 0; r < n; ++r) {
int sum = 0;
for (int c = 0; c < block_exponent; ++c) {
if ((state >> c) & 1) sum += B[r][c];
}
for (int c = block_exponent; c < n; ++c) {
if (high_mask[c]) sum += B[r][c];
}
term = mul_mod(term, sum);
if (term == 0) break;
}
if (((popcnt[static_cast<std::size_t>(state)] + high_ones) & 1) != 0 && term != 0) {
term = MOD - term;
}
result_block[static_cast<std::size_t>(state)] =
add_mod(result_block[static_cast<std::size_t>(state)], term);
}
return;
}
const int sz = high_class_sizes[static_cast<std::size_t>(cls_idx)];
const int base = high_class_offsets[static_cast<std::size_t>(cls_idx)];
for (int ones = 0; ones <= sz; ++ones) {
for (int j = 0; j < sz; ++j) {
high_mask[base + j] = (j >= sz - ones) ? 1 : 0;
}
int ways2 = mul_mod(total_ways, binom[static_cast<std::size_t>(sz)][static_cast<std::size_t>(ones)]);
rec(cls_idx + 1, ways2, high_ones + ones);
}
};
rec(0, 1, 0);
int result = 0;
for (int v : result_block) {
result = add_mod(result, v);
}
if (n & 1) {
result = (result == 0) ? 0 : (MOD - result);
}
return result;
}
std::vector<std::pair<int, int>> row_classes(const std::vector<std::vector<int>>& A) {
std::map<u64, std::vector<int>> groups;
for (int i = 0; i < static_cast<int>(A.size()); ++i) {
u64 mask = 0;
for (int j = 0; j < static_cast<int>(A[i].size()); ++j) {
if (A[i][j] != 0) mask |= (1ULL << j);
}
groups[mask].push_back(i);
}
std::vector<std::pair<int, int>> out;
out.reserve(groups.size());
for (const auto& [k, idxs] : groups) {
(void)k;
out.push_back({idxs[0], static_cast<int>(idxs.size())});
}
return out;
}
template <class F>
void for_each_combination(int n, int k, F&& fn) {
if (k < 0 || k > n) return;
if (k == 0) {
std::vector<int> empty;
fn(empty);
return;
}
std::vector<int> comb(static_cast<std::size_t>(k));
for (int i = 0; i < k; ++i) comb[static_cast<std::size_t>(i)] = i;
while (true) {
fn(comb);
int i = k - 1;
while (i >= 0 && comb[static_cast<std::size_t>(i)] == n - k + i) --i;
if (i < 0) break;
++comb[static_cast<std::size_t>(i)];
for (int j = i + 1; j < k; ++j) {
comb[static_cast<std::size_t>(j)] = comb[static_cast<std::size_t>(j - 1)] + 1;
}
}
}
int P_mod(int n) {
assert(n >= 2 && (n % 2 == 0));
const int half_n = n / 2;
std::vector<int> even_indices, odd_indices;
even_indices.reserve(half_n);
odd_indices.reserve(half_n);
for (int x = 1; x < n; x += 2) even_indices.push_back(x);
for (int x = 2; x <= n; x += 2) odd_indices.push_back(x);
std::vector<std::vector<int>> B(half_n, std::vector<int>(half_n, 0));
for (int c = 0; c < half_n; ++c) {
for (int r = 0; r < half_n; ++r) {
B[c][r] = (std::gcd(even_indices[r], odd_indices[c]) == 1) ? 1 : 0;
}
}
std::vector<std::vector<int>> binom = build_binom(half_n);
std::vector<std::pair<int, int>> rows_and_ways = row_classes(B);
if (!rows_and_ways.empty()) {
int row_to_remove = -1;
for (int i = 0; i < static_cast<int>(rows_and_ways.size()); ++i) {
if (rows_and_ways[static_cast<std::size_t>(i)].second == 1) {
row_to_remove = i;
break;
}
}
if (row_to_remove >= 0) {
rows_and_ways.erase(rows_and_ways.begin() + row_to_remove);
} else {
rows_and_ways[0].second -= 1;
}
}
std::vector<std::vector<int>> BT(half_n, std::vector<int>(half_n, 0));
for (int r = 0; r < half_n; ++r) {
for (int c = 0; c < half_n; ++c) {
BT[r][c] = B[c][r];
}
}
std::vector<std::pair<int, int>> cols_and_ways = row_classes(BT);
int result = 0;
for (int included = 0; included < half_n; ++included) {
for_each_combination(static_cast<int>(rows_and_ways.size()), included,
[&](const std::vector<int>& rows_pick) {
i64 rows_ways = 1;
std::vector<int> included_rows;
included_rows.reserve(rows_pick.size());
for (int idx : rows_pick) {
const auto [row_idx, ways] = rows_and_ways[static_cast<std::size_t>(idx)];
rows_ways *= ways;
included_rows.push_back(row_idx);
}
std::vector<char> in_row(half_n, 0);
for (int r : included_rows) in_row[r] = 1;
std::vector<int> not_rows;
not_rows.reserve(half_n - included);
for (int r = 0; r < half_n; ++r) {
if (!in_row[r]) not_rows.push_back(r);
}
for_each_combination(static_cast<int>(cols_and_ways.size()), included,
[&](const std::vector<int>& cols_pick) {
i64 ways = rows_ways;
std::vector<int> included_cols;
included_cols.reserve(cols_pick.size());
for (int idx : cols_pick) {
const auto [col_idx, cw] =
cols_and_ways[static_cast<std::size_t>(idx)];
ways *= cw;
included_cols.push_back(col_idx);
}
std::vector<std::vector<int>> det_mat(
static_cast<std::size_t>(included),
std::vector<int>(static_cast<std::size_t>(included),
0));
for (int i = 0; i < included; ++i) {
for (int j = 0; j < included; ++j) {
det_mat[static_cast<std::size_t>(i)]
[static_cast<std::size_t>(j)] =
B[included_rows[static_cast<std::size_t>(i)]]
[included_cols[static_cast<std::size_t>(j)]];
}
}
int det_val = determinant_mod(det_mat);
if (det_val == 0) return;
det_val = mul_mod(det_val, det_val);
if (included & 1) {
det_val = (det_val == 0) ? 0 : (MOD - det_val);
}
std::vector<char> in_col(half_n, 0);
for (int c : included_cols) in_col[c] = 1;
std::vector<int> not_cols;
not_cols.reserve(half_n - included);
for (int c = 0; c < half_n; ++c) {
if (!in_col[c]) not_cols.push_back(c);
}
const int rest = half_n - included;
std::vector<std::vector<int>> per_mat(
static_cast<std::size_t>(rest),
std::vector<int>(static_cast<std::size_t>(rest), 0));
for (int i = 0; i < rest; ++i) {
for (int j = 0; j < rest; ++j) {
per_mat[static_cast<std::size_t>(i)]
[static_cast<std::size_t>(j)] =
B[not_rows[static_cast<std::size_t>(i)]]
[not_cols[static_cast<std::size_t>(j)]];
}
}
int per_val = permanent_mod_with_column_classes(per_mat, binom);
per_val = mul_mod(per_val, per_val);
int term = mul_mod(static_cast<int>(ways % MOD),
mul_mod(det_val, per_val));
result = add_mod(result, term);
});
});
}
return result;
}
int brute_count(int n) {
std::vector<int> v;
for (int x = 2; x <= n; ++x) v.push_back(x);
int cnt = 0;
std::sort(v.begin(), v.end());
do {
bool ok = true;
for (int i = 1; i < static_cast<int>(v.size()); ++i) {
if (std::gcd(v[i - 1], v[i]) != 1) {
ok = false;
break;
}
}
if (ok) ++cnt;
} while (std::next_permutation(v.begin(), v.end()));
return cnt;
}
} // namespace
int main() {
assert(P_mod(4) == 2);
assert(P_mod(10) == 576);
assert(P_mod(6) == brute_count(6));
assert(P_mod(8) == brute_count(8));
std::cout << P_mod(34) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 83456729
def mul_mod(a, b): return a*b%MOD
def add_mod(a, b):
s = a+b; return s-MOD if s>=MOD else s
def sub_mod(a, b):
s = a-b; return s+MOD if s<0 else s
def mod_pow(a, e):
r = 1; a %= MOD
while e > 0:
if e&1: r=r*a%MOD
a=a*a%MOD; e >>= 1
return r
def mod_inv(a): return mod_pow(a, MOD-2)
def det_mod(A):
n = len(A)
if n == 0: return 1
if n == 1: return A[0][0]%MOD
A = [row[:] for row in A]; ri = 1
for i in range(n-2):
piv = -1
for r in range(i, n):
if A[r][i]%MOD: piv=r; break
if piv<0: return 0
if piv!=i: A[piv],A[i]=A[i],A[piv]; ri=MOD-ri if ri else 0
Aii = A[i][i]%MOD; ri = mul_mod(ri, mod_pow(Aii, n-i-2))
for r in range(i+1,n):
x = A[r][i]%MOD; A[r][i]=0
for c in range(i+1,n): A[r][c]=sub_mod(mul_mod(Aii,A[r][c]),mul_mod(x,A[i][c]))
last = sub_mod(mul_mod(A[n-2][n-2],A[n-1][n-1]),mul_mod(A[n-2][n-1],A[n-1][n-2]))
return mul_mod(last, mod_inv(ri))
def perm_mod(A, binom):
n = len(A)
if n == 0: return 1
cc = {}
for col in range(n):
mask = 0
for r in range(n):
if A[r][col]: mask |= 1<<r
cc.setdefault(mask, []).append(col)
classes = []
for inds in cc.values():
if len(inds) >= 3: classes.append(inds)
else:
for idx in inds: classes.append([idx])
classes.sort(key=len)
order = [c for cl in classes for c in cl]
B = [[A[r][order[c]] for c in range(n)] for r in range(n)]
csz = [len(cl) for cl in classes]
be = 0; sc = 0
while sc < len(csz) and be+csz[sc] < 10: be += csz[sc]; sc += 1
bs = 1<<be
hcs = csz[sc:]; hco = []; off = be
for s in hcs: hco.append(off); off += s
result = [0]
hmask = [0]*n
pc = [0]*bs
for s in range(1,bs): pc[s] = pc[s>>1]+(s&1)
def rec(ci, tw, ho):
if ci == len(hcs):
for state in range(bs):
term = tw
for r in range(n):
sm = 0
for c in range(be):
if (state>>c)&1: sm += B[r][c]
for c in range(be, n):
if hmask[c]: sm += B[r][c]
term = mul_mod(term, sm)
if term == 0: break
if (pc[state]+ho)&1 and term: term = MOD-term
result[0] = add_mod(result[0], term)
return
sz = hcs[ci]; base = hco[ci]
for ones in range(sz+1):
for j in range(sz): hmask[base+j] = 1 if j>=sz-ones else 0
w2 = mul_mod(tw, binom[sz][ones])
rec(ci+1, w2, ho+ones)
rec(0, 1, 0)
r = result[0]
if n&1: r = MOD-r if r else 0
return r
nn = 34; half = nn//2
even_idx = list(range(1, nn, 2)); odd_idx = list(range(2, nn+1, 2))
Bm = [[1 if math.gcd(even_idx[r], odd_idx[c])==1 else 0 for r in range(half)] for c in range(half)]
binom = [[0]*(half+1) for _ in range(half+1)]
for i in range(half+1):
binom[i][0] = 1
for j in range(1, i+1): binom[i][j] = binom[i-1][j-1]+binom[i-1][j]
def row_classes(M):
groups = {}
for i in range(len(M)):
mask = sum(1<<j for j in range(len(M[i])) if M[i][j])
groups.setdefault(mask, []).append(i)
return [(idxs[0], len(idxs)) for idxs in groups.values()]
rw = row_classes(Bm)
# Remove one row
rm = -1
for i in range(len(rw)):
if rw[i][1] == 1: rm=i; break
if rm >= 0: rw = rw[:rm]+rw[rm+1:]
else: rw = [(rw[0][0],rw[0][1]-1)]+rw[1:]
BT = [[Bm[c][r] for c in range(half)] for r in range(half)]
cw = row_classes(BT)
from itertools import combinations
result = 0
for inc in range(half):
for rp in combinations(range(len(rw)), inc):
rways = 1; irows = []
for idx in rp: rways *= rw[idx][1]; irows.append(rw[idx][0])
inr = set(irows)
nrows = [r for r in range(half) if r not in inr]
for cp in combinations(range(len(cw)), inc):
ways = rways; icols = []
for idx in cp: ways *= cw[idx][1]; icols.append(cw[idx][0])
dm = [[Bm[irows[i]][icols[j]] for j in range(inc)] for i in range(inc)]
dv = det_mod(dm)
if dv == 0: continue
dv = mul_mod(dv, dv)
if inc&1: dv = MOD-dv if dv else 0
inc2 = set(icols); ncols = [c for c in range(half) if c not in inc2]
rest = half-inc
pm = [[Bm[nrows[i]][ncols[j]] for j in range(rest)] for i in range(rest)]
pv = perm_mod(pm, binom); pv = mul_mod(pv, pv)
result = add_mod(result, mul_mod(ways%MOD, mul_mod(dv, pv)))
return str(result)
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler886 {
static final int MOD = 83456729;
static int addMod(int a, int b) {
int s = a + b;
if (s >= MOD)
s -= MOD;
return s;
}
static int subMod(int a, int b) {
int s = a - b;
if (s < 0)
s += MOD;
return s;
}
static int mulMod(long a, long b) {
return (int) ((a * b) % MOD);
}
static int modPow(int a, long e) {
long r = 1;
long x = a % MOD;
if (x < 0)
x += MOD;
while (e > 0) {
if ((e & 1) == 1)
r = (r * x) % MOD;
x = (x * x) % MOD;
e >>= 1;
}
return (int) r;
}
static int modInv(int a) {
return modPow(a, MOD - 2);
}
static int[][] buildBinom(int n) {
int[][] C = new int[n + 1][];
C[0] = new int[] { 1 };
for (int i = 1; i <= n; ++i) {
C[i] = new int[i + 1];
Arrays.fill(C[i], 1);
for (int j = 1; j < i; ++j) {
C[i][j] = addMod(C[i - 1][j - 1], C[i - 1][j]);
}
}
return C;
}
static int determinantMod(int[][] origA) {
int n = origA.length;
if (n == 0)
return 1;
if (n == 1)
return origA[0][0] % MOD;
int[][] A = new int[n][n];
for (int i = 0; i < n; i++)
A[i] = origA[i].clone();
int resultInverse = 1;
for (int i = 0; i < n - 2; ++i) {
int rowToUse = -1;
for (int r = i; r < n; ++r) {
if (A[r][i] % MOD != 0) {
rowToUse = r;
break;
}
}
if (rowToUse < 0)
return 0;
if (i != rowToUse) {
int[] temp = A[rowToUse];
A[rowToUse] = A[i];
A[i] = temp;
resultInverse = (resultInverse == 0) ? 0 : MOD - resultInverse;
}
int Aii = A[i][i] % MOD;
resultInverse = mulMod(resultInverse, modPow(Aii, n - i - 2));
for (int r = i + 1; r < n; ++r) {
int x = A[r][i] % MOD;
A[r][i] = 0;
for (int c = i + 1; c < n; ++c) {
int v = subMod(mulMod(Aii, A[r][c]), mulMod(x, A[i][c]));
A[r][c] = v;
}
}
}
int last = subMod(mulMod(A[n - 2][n - 2], A[n - 1][n - 1]),
mulMod(A[n - 2][n - 1], A[n - 1][n - 2]));
return mulMod(last, modInv(resultInverse));
}
static int permanentModWithColumnClasses(int[][] A, int[][] binom) {
int n = A.length;
if (n == 0)
return 1;
Map<Long, List<Integer>> columnClasses = new HashMap<>();
for (int col = 0; col < n; ++col) {
long mask = 0;
for (int r = 0; r < n; ++r) {
if (A[r][col] != 0)
mask |= (1L << r);
}
columnClasses.computeIfAbsent(mask, k -> new ArrayList<>()).add(col);
}
List<List<Integer>> classes = new ArrayList<>();
for (List<Integer> inds : columnClasses.values()) {
if (inds.size() >= 3) {
classes.add(inds);
} else {
for (int idx : inds) {
List<Integer> single = new ArrayList<>();
single.add(idx);
classes.add(single);
}
}
}
classes.sort(Comparator.comparingInt(List::size));
int[] newColOrder = new int[n];
int colIdx = 0;
for (List<Integer> c : classes) {
for (int idx : c)
newColOrder[colIdx++] = idx;
}
int[][] B = new int[n][n];
for (int r = 0; r < n; ++r) {
for (int c = 0; c < n; ++c) {
B[r][c] = A[r][newColOrder[c]];
}
}
int[] classSizes = new int[classes.size()];
for (int i = 0; i < classes.size(); i++)
classSizes[i] = classes.get(i).size();
int defaultBlockExponent = 10;
int smallestClassCount = 0;
int smallestClassIndexCount = 0;
while (smallestClassCount < classSizes.length &&
smallestClassIndexCount + classSizes[smallestClassCount] < defaultBlockExponent) {
smallestClassIndexCount += classSizes[smallestClassCount];
++smallestClassCount;
}
int blockExponent = smallestClassIndexCount;
int blockSize = 1 << blockExponent;
int numHighClasses = classSizes.length - smallestClassCount;
int[] highClassSizes = new int[numHighClasses];
int[] highClassOffsets = new int[numHighClasses];
int off = blockExponent;
for (int i = 0; i < numHighClasses; ++i) {
highClassSizes[i] = classSizes[smallestClassCount + i];
highClassOffsets[i] = off;
off += highClassSizes[i];
}
int[] resultBlock = new int[blockSize];
int[] highMask = new int[n];
int[] popcnt = new int[blockSize];
for (int s = 1; s < blockSize; ++s)
popcnt[s] = popcnt[s >> 1] + (s & 1);
rec(0, 1, 0, highClassSizes, highClassOffsets, blockExponent, blockSize, n, B, binom, highMask, popcnt,
resultBlock);
int result = 0;
for (int v : resultBlock) {
result = addMod(result, v);
}
if ((n & 1) != 0) {
result = (result == 0) ? 0 : (MOD - result);
}
return result;
}
static void rec(int clsIdx, int totalWays, int highOnes, int[] highClassSizes, int[] highClassOffsets,
int blockExponent, int blockSize, int n, int[][] B, int[][] binom, int[] highMask, int[] popcnt,
int[] resultBlock) {
if (clsIdx == highClassSizes.length) {
for (int state = 0; state < blockSize; ++state) {
int term = totalWays;
for (int r = 0; r < n; ++r) {
int sum = 0;
for (int c = 0; c < blockExponent; ++c) {
if (((state >> c) & 1) != 0)
sum += B[r][c];
}
for (int c = blockExponent; c < n; ++c) {
if (highMask[c] != 0)
sum += B[r][c];
}
term = mulMod(term, sum);
if (term == 0)
break;
}
if (((popcnt[state] + highOnes) & 1) != 0 && term != 0) {
term = MOD - term;
}
resultBlock[state] = addMod(resultBlock[state], term);
}
return;
}
int sz = highClassSizes[clsIdx];
int base = highClassOffsets[clsIdx];
for (int ones = 0; ones <= sz; ++ones) {
for (int j = 0; j < sz; ++j) {
highMask[base + j] = (j >= sz - ones) ? 1 : 0;
}
int ways2 = mulMod(totalWays, binom[sz][ones]);
rec(clsIdx + 1, ways2, highOnes + ones, highClassSizes, highClassOffsets, blockExponent, blockSize, n, B,
binom, highMask, popcnt, resultBlock);
}
}
static class RowClass {
int idx;
int size;
RowClass(int idx, int size) {
this.idx = idx;
this.size = size;
}
}
static List<RowClass> rowClasses(int[][] A) {
Map<Long, List<Integer>> groups = new HashMap<>();
for (int i = 0; i < A.length; ++i) {
long mask = 0;
for (int j = 0; j < A[i].length; ++j) {
if (A[i][j] != 0)
mask |= (1L << j);
}
groups.computeIfAbsent(mask, k -> new ArrayList<>()).add(i);
}
List<RowClass> out = new ArrayList<>();
for (List<Integer> idxs : groups.values()) {
out.add(new RowClass(idxs.get(0), idxs.size()));
}
return out;
}
interface CombCallback {
void call(int[] comb);
}
static void forEachCombination(int n, int k, CombCallback fn) {
if (k < 0 || k > n)
return;
if (k == 0) {
fn.call(new int[0]);
return;
}
int[] comb = new int[k];
for (int i = 0; i < k; ++i)
comb[i] = i;
while (true) {
fn.call(comb);
int i = k - 1;
while (i >= 0 && comb[i] == n - k + i)
--i;
if (i < 0)
break;
++comb[i];
for (int j = i + 1; j < k; ++j) {
comb[j] = comb[j - 1] + 1;
}
}
}
static int gcd(int a, int b) {
return b == 0 ? a : gcd(b, a % b);
}
static int PMod(int n) {
int halfN = n / 2;
int[] evenIndices = new int[halfN];
int[] oddIndices = new int[halfN];
int evenIdx = 0, oddIdx = 0;
for (int x = 1; x < n; x += 2)
evenIndices[evenIdx++] = x;
for (int x = 2; x <= n; x += 2)
oddIndices[oddIdx++] = x;
int[][] B = new int[halfN][halfN];
for (int c = 0; c < halfN; ++c) {
for (int r = 0; r < halfN; ++r) {
B[c][r] = (gcd(evenIndices[r], oddIndices[c]) == 1) ? 1 : 0;
}
}
int[][] binom = buildBinom(halfN);
List<RowClass> rowsAndWays = rowClasses(B);
if (!rowsAndWays.isEmpty()) {
int rowToRemove = -1;
for (int i = 0; i < rowsAndWays.size(); ++i) {
if (rowsAndWays.get(i).size == 1) {
rowToRemove = i;
break;
}
}
if (rowToRemove >= 0) {
rowsAndWays.remove(rowToRemove);
} else {
rowsAndWays.get(0).size -= 1;
}
}
int[][] BT = new int[halfN][halfN];
for (int r = 0; r < halfN; ++r) {
for (int c = 0; c < halfN; ++c) {
BT[r][c] = B[c][r];
}
}
List<RowClass> colsAndWays = rowClasses(BT);
int[] result = { 0 };
for (int included = 0; included < halfN; ++included) {
final int inc = included;
forEachCombination(rowsAndWays.size(), included, rowsPick -> {
long rowsWays = 1;
int[] includedRows = new int[inc];
for (int i = 0; i < inc; i++) {
RowClass rc = rowsAndWays.get(rowsPick[i]);
rowsWays = (rowsWays * rc.size) % MOD;
includedRows[i] = rc.idx;
}
boolean[] inRow = new boolean[halfN];
for (int r : includedRows)
inRow[r] = true;
int[] notRows = new int[halfN - inc];
int notRowIdx = 0;
for (int r = 0; r < halfN; ++r) {
if (!inRow[r])
notRows[notRowIdx++] = r;
}
final long finalRowsWays = rowsWays;
forEachCombination(colsAndWays.size(), inc, colsPick -> {
long ways = finalRowsWays;
int[] includedCols = new int[inc];
for (int i = 0; i < inc; i++) {
RowClass cc = colsAndWays.get(colsPick[i]);
ways = (ways * cc.size) % MOD;
includedCols[i] = cc.idx;
}
int[][] detMat = new int[inc][inc];
for (int i = 0; i < inc; ++i) {
for (int j = 0; j < inc; ++j) {
detMat[i][j] = B[includedRows[i]][includedCols[j]];
}
}
int detVal = determinantMod(detMat);
if (detVal == 0)
return;
detVal = mulMod(detVal, detVal);
if ((inc & 1) != 0) {
detVal = (detVal == 0) ? 0 : (MOD - detVal);
}
boolean[] inCol = new boolean[halfN];
for (int c : includedCols)
inCol[c] = true;
int[] notCols = new int[halfN - inc];
int notColIdx = 0;
for (int c = 0; c < halfN; ++c) {
if (!inCol[c])
notCols[notColIdx++] = c;
}
int rest = halfN - inc;
int[][] perMat = new int[rest][rest];
for (int i = 0; i < rest; ++i) {
for (int j = 0; j < rest; ++j) {
perMat[i][j] = B[notRows[i]][notCols[j]];
}
}
int perVal = permanentModWithColumnClasses(perMat, binom);
perVal = mulMod(perVal, perVal);
int term = mulMod(ways, mulMod(detVal, perVal));
result[0] = addMod(result[0], term);
});
});
}
return result[0];
}
public static String solve() {
return Integer.toString(PMod(34));
}
public static void main(String[] args) {
System.out.println(solve());
}
}