Problem 843: Periodic Circles

Project Euler Problem 843 Solution

EulerSolve provides an optimized solution for Project Euler Problem 843, Periodic Circles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $n \ge 3$, consider binary states on a circle of length $n$. One update step replaces every entry by the xor of its two neighbors, so the dynamics are linear over $\mathbb F_2$. A positive integer $t$ belongs to $\mathcal P(n)$ if some state has least period $t$ under this update. The global target is $$\mathcal G(N)=\bigcup_{n=3}^{N}\mathcal P(n),\qquad S(N)=\sum_{v\in\mathcal G(N)} v.$$ The key point is that the implementations never simulate all states. They turn the update rule into multiplication in a quotient ring, analyze each irreducible factor separately, and reconstruct the possible periods by least common multiples. Mathematical Approach The natural state space for a circle of size $n$ is $$R_n=\mathbb F_2[x]/(x^n+1).$$ A binary state $(a_0,\dots,a_{n-1})$ is encoded as $$A(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\in R_n.$$ Multiplication by $x$ shifts left, multiplication by $x^{n-1}=x^{-1}$ shifts right, so the update operator is simply $$T_n(A)=(x+x^{-1})A=(x+x^{n-1})A.$$ Step 1: Split Off the Power of Two Write $$n=2^k m,\qquad m\text{ odd}.$$ In characteristic $2$, the Frobenius identity gives $$x^n+1=x^{2^k m}+1=\left(x^m+1\right)^{2^k}.$$ Because $m$ is odd, the derivative of $x^m+1$ is $x^{m-1}$, so $x^m+1$ is square-free....

Detailed mathematical approach

Problem Summary

For each integer $n \ge 3$, consider binary states on a circle of length $n$. One update step replaces every entry by the xor of its two neighbors, so the dynamics are linear over $\mathbb F_2$. A positive integer $t$ belongs to $\mathcal P(n)$ if some state has least period $t$ under this update.

The global target is

$$\mathcal G(N)=\bigcup_{n=3}^{N}\mathcal P(n),\qquad S(N)=\sum_{v\in\mathcal G(N)} v.$$

The key point is that the implementations never simulate all states. They turn the update rule into multiplication in a quotient ring, analyze each irreducible factor separately, and reconstruct the possible periods by least common multiples.

Mathematical Approach

The natural state space for a circle of size $n$ is

$$R_n=\mathbb F_2[x]/(x^n+1).$$

A binary state $(a_0,\dots,a_{n-1})$ is encoded as

$$A(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\in R_n.$$

Multiplication by $x$ shifts left, multiplication by $x^{n-1}=x^{-1}$ shifts right, so the update operator is simply

$$T_n(A)=(x+x^{-1})A=(x+x^{n-1})A.$$

Step 1: Split Off the Power of Two

Write

$$n=2^k m,\qquad m\text{ odd}.$$

In characteristic $2$, the Frobenius identity gives

$$x^n+1=x^{2^k m}+1=\left(x^m+1\right)^{2^k}.$$

Because $m$ is odd, the derivative of $x^m+1$ is $x^{m-1}$, so $x^m+1$ is square-free. Hence

$$x^m+1=\prod_{j=1}^{r} g_j(x)$$

with pairwise distinct irreducible factors $g_j$. Therefore

$$x^n+1=\prod_{j=1}^{r} g_j(x)^{2^k}.$$

Step 2: Decompose the State Space

Since the factors $g_j^{2^k}$ are pairwise coprime, the Chinese remainder theorem gives

$$R_n\cong \bigoplus_{j=1}^{r}\mathbb F_2[x]/\left(g_j(x)^{2^k}\right).$$

This means the dynamics split into independent local components. If one local component has period $u_j$, then the whole state has period

$$\operatorname{lcm}(u_1,\dots,u_r).$$

So the real job is to determine the possible local periods contributed by each irreducible factor.

Step 3: Reduce to a Finite Field Element

First reduce modulo $g_j$, not modulo the higher power $g_j^{2^k}$. In the field

$$K_j=\mathbb F_2[x]/(g_j(x))$$

we have $x^m=1$, so $x^{-1}=x^{m-1}$. The local multiplier becomes

$$\lambda_j=x+x^{-1}=x+x^{m-1}\pmod{g_j}.$$

There are three cases.

If $\lambda_j=0$, the lifted operator will be nilpotent on that branch and can only contribute period $1$.

If $\lambda_j=1$, the semisimple part is trivial and only powers of $2$ can appear after lifting.

If $\lambda_j\neq 0,1$, then $\lambda_j\in K_j^\times$ has a multiplicative order $r_j$.

The implementations compute $r_j$ by finding the smallest $e_j$ such that

$$\lambda_j^{2^{e_j}}=\lambda_j.$$

This means $\lambda_j$ lies in the subfield $\mathbb F_{2^{e_j}}$, so

$$r_j\mid 2^{e_j}-1.$$

Starting from $2^{e_j}-1$, they divide out prime factors whenever the corresponding power test still gives $1$, which reduces the candidate to the true multiplicative order.

Step 4: Lift from $g_j$ to $g_j^{2^k}$

Now return to the local ring $\mathbb F_2[x]/(g_j^{2^k})$. If $\lambda_j\neq 0$, the lifted multiplier can be written as

$$\lambda_j(1+u),\qquad u^{2^k}=0,$$

where $u$ lives in the nilpotent part created by the repeated factor $g_j^{2^k}$. In characteristic $2$,

$$ (1+u)^{2^s}=1+u^{2^s}. $$

So the nilpotent correction contributes only powers of $2$. As the state moves through the $g_j$-adic filtration, every exponent $2^s$ with $0\le s\le k$ can occur when the semisimple part is nontrivial, and every exponent $2^s$ with $1\le s\le k$ can occur when the semisimple part is $1$.

Therefore the local period options are

$$ \mathcal O_j=\{1\}\cup \begin{cases} \varnothing, & \lambda_j=0,\\ \{2^s:1\le s\le k\}, & \lambda_j=1,\\ \{r_j2^s:0\le s\le k\}, & \operatorname{ord}(\lambda_j)=r_j>1. \end{cases} $$

The initial $1$ accounts for choosing a local state whose periodic part is trivial.

Step 5: Reconstruct the Period Set for a Fixed $n$

Because the local factors are independent, the full period set for this circle size is exactly

$$\mathcal P(n)=\left\{\operatorname{lcm}(u_1,\dots,u_r):u_j\in\mathcal O_j\right\}.$$

This formula is the entire algorithm in one line: factor $x^m+1$, compute one local option set per irreducible factor, then close under lcm.

Step 6: Worked Examples

For $n=5$, we have $k=0$ and

$$x^5+1=(x+1)(x^4+x^3+x^2+x+1).$$

On the factor $x+1$, we substitute $x=1$, so $\lambda=1+1=0$, contributing only $1$.

On the quartic factor, let $\alpha$ be a root. Since $\alpha^5=1$,

$$\lambda=\alpha+\alpha^4.$$

Using $1+\alpha+\alpha^2+\alpha^3+\alpha^4=0$, we get

$$\lambda^2+\lambda+1=\alpha^2+\alpha^3+\alpha+\alpha^4+1=0,$$

so $\lambda$ has multiplicative order $3$. Because $k=0$, the second factor contributes $\{1,3\}$, hence

$$\mathcal P(5)=\{1,3\}.$$

For $n=6$, we have $n=2\cdot 3$, so $k=1$ and

$$x^3+1=(x+1)(x^2+x+1).$$

The linear factor again gives $\lambda=0$. On the quadratic factor, $x^2+x+1=0$ implies $x^2=x+1$, so

$$\lambda=x+x^2=1.$$

Now the local periods are $1$ and $2$, giving

$$\mathcal P(6)=\{1,2\}.$$

Step 7: Union Over All Circle Sizes

After computing $\mathcal P(n)$ for every $3\le n\le N$, we insert every period into one global set

$$\mathcal G(N)=\bigcup_{n=3}^{N}\mathcal P(n),$$

and the final answer is

$$S(N)=\sum_{v\in\mathcal G(N)} v.$$

This is why the program cares only about distinct periods, not how many states realize each one.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first precompute ordinary prime numbers large enough to factor the integers $2^e-1$ that appear during order reduction.

They then implement polynomial arithmetic over $\mathbb F_2$: degree, remainder, quotient, greatest common divisor, multiplication modulo a polynomial, and fast modular exponentiation. With that toolkit they apply Berlekamp factorization to $x^m+1$ for each odd part $m$.

For every irreducible factor, the implementation evaluates the local multiplier $x+x^{m-1}$, detects the special cases $0$ and $1$, and otherwise computes the multiplicative order by repeated Frobenius squaring followed by divisor stripping inside $2^e-1$.

Those local orders are converted into the option sets $\mathcal O_j$. An iterative lcm-closure pass combines them into $\mathcal P(n)$. The order data are cached by odd part $m$, so the sizes $m,2m,4m,\dots$ reuse the same field-factor analysis and differ only in the power-of-two lift range.

Finally the implementations take the union of all period sets for $3\le n\le N$ and sum the distinct values. The checkpoints $\mathcal P(5)=\{1,3\}$, $\mathcal P(6)=\{1,2\}$, and $S(30)=20381$ match the mathematical derivation above.

Complexity Analysis

Let $U=\{m\le N:m\text{ odd}\}$. Each odd part is factored only once, so the dominant algebraic work is the sum over $m\in U$ of factoring $x^m+1$ over $\mathbb F_2$. In this direct implementation, the Berlekamp linear algebra on degree-$m$ polynomials is cubic in $m$ up to bit-operation details, which is easily fast enough for $N=100$.

Once the factors are known, order computation for one factor of degree $d$ uses repeated squaring in the field together with tests on divisors of $2^e-1$, where $e\le d$. The period-set construction is combinatorial rather than algebraic: if the local option sets have sizes $|\mathcal O_1|,\dots,|\mathcal O_r|$, then building $\mathcal P(n)$ costs the iterative lcm-closure work induced by those choices.

Memory usage is modest. The program stores cached order lists for odd parts already processed, temporary polynomial data for one factorization, the current period set for a given $n$, and the final global set of distinct periods.

Footnotes and References

Problem page: https://projecteuler.net/problem=843
Berlekamp's algorithm: Wikipedia — Berlekamp's algorithm
Finite field: Wikipedia — Finite field
Chinese remainder theorem: Wikipedia — Chinese remainder theorem
Multiplicative order: Wikipedia — Multiplicative order

Problem 843 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <unordered_map>
#include <unordered_set>
#include <vector>

using namespace std;

namespace {

struct Poly {
    uint64_t lo = 0;
    uint64_t hi = 0;

    bool is_zero() const {
        return lo == 0 && hi == 0;
    }
    bool is_one() const {
        return lo == 1 && hi == 0;
    }
};

Poly operator^(const Poly& a, const Poly& b) {
    return {a.lo ^ b.lo, a.hi ^ b.hi};
}

Poly& operator^=(Poly& a, const Poly& b) {
    a.lo ^= b.lo;
    a.hi ^= b.hi;
    return a;
}

bool poly_equal(const Poly& a, const Poly& b) {
    return a.lo == b.lo && a.hi == b.hi;
}

int poly_degree(const Poly& p) {
    if (p.hi != 0) {
        return 64 + 63 - __builtin_clzll(p.hi);
    }
    if (p.lo != 0) {
        return 63 - __builtin_clzll(p.lo);
    }
    return -1;
}

bool poly_get_bit(const Poly& p, int pos) {
    if (pos < 64) return (p.lo >> pos) & 1ULL;
    return (p.hi >> (pos - 64)) & 1ULL;
}

void poly_toggle_bit(Poly& p, int pos) {
    if (pos < 64) {
        p.lo ^= 1ULL << pos;
    } else {
        p.hi ^= 1ULL << (pos - 64);
    }
}

Poly poly_shift_left(const Poly& p, int k) {
    unsigned __int128 v = (static_cast<unsigned __int128>(p.hi) << 64) | p.lo;
    v <<= k;
    return {static_cast<uint64_t>(v), static_cast<uint64_t>(v >> 64)};
}

Poly poly_shift_right1(const Poly& p) {
    Poly r;
    r.lo = (p.lo >> 1) | (p.hi << 63);
    r.hi = p.hi >> 1;
    return r;
}

Poly poly_shift_left1(const Poly& p) {
    Poly r;
    r.hi = (p.hi << 1) | (p.lo >> 63);
    r.lo = p.lo << 1;
    return r;
}

Poly poly_mod(Poly a, const Poly& mod) {
    int deg_mod = poly_degree(mod);
    if (deg_mod < 0) return a;
    for (;;) {
        int deg_a = poly_degree(a);
        if (deg_a < deg_mod) break;
        int shift = deg_a - deg_mod;
        a ^= poly_shift_left(mod, shift);
    }
    return a;
}

Poly poly_div(Poly a, const Poly& mod) {
    Poly q;
    int deg_mod = poly_degree(mod);
    for (;;) {
        int deg_a = poly_degree(a);
        if (deg_a < deg_mod) break;
        int shift = deg_a - deg_mod;
        poly_toggle_bit(q, shift);
        a ^= poly_shift_left(mod, shift);
    }
    return q;
}

Poly poly_gcd(Poly a, Poly b) {
    while (!b.is_zero()) {
        Poly r = poly_mod(a, b);
        a = b;
        b = r;
    }
    return a;
}

Poly poly_mul_mod(Poly a, Poly b, const Poly& mod) {
    if (a.is_zero() || b.is_zero()) return {0, 0};
    a = poly_mod(a, mod);
    b = poly_mod(b, mod);
    int deg_mod = poly_degree(mod);
    Poly res;
    while (!b.is_zero()) {
        if (b.lo & 1ULL) res ^= a;
        b = poly_shift_right1(b);
        a = poly_shift_left1(a);
        if (poly_get_bit(a, deg_mod)) a ^= mod;
    }
    return res;
}

Poly poly_pow_mod(Poly base, uint64_t exp, const Poly& mod) {
    Poly res{1, 0};
    while (exp > 0) {
        if (exp & 1ULL) res = poly_mul_mod(res, base, mod);
        base = poly_mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return res;
}

vector<uint32_t> sieve_primes(int limit) {
    vector<bool> is_prime(limit + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; i * i <= limit; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= limit; j += i) is_prime[j] = false;
    }
    vector<uint32_t> primes;
    for (int i = 2; i <= limit; ++i) {
        if (is_prime[i]) primes.push_back(static_cast<uint32_t>(i));
    }
    return primes;
}

vector<uint64_t> factorize_uint64(uint64_t n, const vector<uint32_t>& primes) {
    vector<uint64_t> factors;
    uint64_t temp = n;
    for (uint32_t p : primes) {
        uint64_t pp = static_cast<uint64_t>(p);
        if (pp * pp > temp) break;
        if (temp % pp == 0) {
            factors.push_back(pp);
            while (temp % pp == 0) temp /= pp;
        }
    }
    if (temp > 1) factors.push_back(temp);
    return factors;
}

vector<Poly> berlekamp_factor(const Poly& f) {
    int m = poly_degree(f);
    if (m <= 1) return {f};

    vector<Poly> rows(m, Poly{0, 0});
    Poly x_poly{2, 0};
    Poly x2 = poly_mul_mod(x_poly, x_poly, f);
    Poly power{1, 0};

    // Build Q matrix columns from x^(2j) mod f, then compute Q - I.
    for (int col = 0; col < m; ++col) {
        for (int row = 0; row < m; ++row) {
            if (poly_get_bit(power, row)) {
                poly_toggle_bit(rows[row], col);
            }
        }
        power = poly_mul_mod(power, x2, f);
    }
    for (int i = 0; i < m; ++i) poly_toggle_bit(rows[i], i);

    // Row-reduce to RREF.
    vector<int> pivot_col;
    pivot_col.reserve(m);
    int rank = 0;
    for (int col = 0; col < m; ++col) {
        int sel = -1;
        for (int r = rank; r < m; ++r) {
            if (poly_get_bit(rows[r], col)) {
                sel = r;
                break;
            }
        }
        if (sel == -1) continue;
        swap(rows[rank], rows[sel]);
        pivot_col.push_back(col);
        for (int r = 0; r < m; ++r) {
            if (r != rank && poly_get_bit(rows[r], col)) {
                rows[r] ^= rows[rank];
            }
        }
        ++rank;
    }

    vector<bool> is_pivot(m, false);
    for (int i = 0; i < rank; ++i) is_pivot[pivot_col[i]] = true;

    vector<Poly> basis;
    basis.reserve(m - rank);
    for (int col = 0; col < m; ++col) {
        if (is_pivot[col]) continue;
        Poly vec{0, 0};
        poly_toggle_bit(vec, col);
        for (int i = 0; i < rank; ++i) {
            int p = pivot_col[i];
            if (poly_get_bit(rows[i], col)) poly_toggle_bit(vec, p);
        }
        basis.push_back(vec);
    }

    vector<Poly> factors;
    factors.push_back(f);
    for (const Poly& b : basis) {
        if (b.is_zero() || b.is_one()) continue;
        vector<Poly> next;
        next.reserve(factors.size() * 2);
        for (const Poly& h : factors) {
            if (poly_degree(h) <= 1) {
                next.push_back(h);
                continue;
            }
            Poly g = poly_gcd(h, b);
            if (!g.is_one() && !poly_equal(g, h)) {
                Poly q = poly_div(h, g);
                next.push_back(g);
                next.push_back(q);
                continue;
            }
            Poly b1 = b;
            poly_toggle_bit(b1, 0);
            g = poly_gcd(h, b1);
            if (!g.is_one() && !poly_equal(g, h)) {
                Poly q = poly_div(h, g);
                next.push_back(g);
                next.push_back(q);
            } else {
                next.push_back(h);
            }
        }
        factors.swap(next);
    }
    return factors;
}

uint64_t order_in_factor(const Poly& lambda, const Poly& mod, const vector<uint32_t>& primes) {
    if (lambda.is_zero()) return 0;
    if (lambda.is_one()) return 1;

    int deg = poly_degree(mod);
    Poly cur = lambda;
    int e = 0;
    do {
        cur = poly_mul_mod(cur, cur, mod);
        ++e;
    } while (!poly_equal(cur, lambda) && e <= deg);

    if (e > deg || e > 63) {
        cerr << "Validation failure: unexpected field degree " << e << '\n';
        exit(1);
    }

    uint64_t order = (e == 64) ? ~0ULL : ((1ULL << e) - 1ULL);
    uint64_t reduced = order;
    vector<uint64_t> factors = factorize_uint64(order, primes);
    for (uint64_t p : factors) {
        while (reduced % p == 0) {
            uint64_t trial = reduced / p;
            if (poly_pow_mod(lambda, trial, mod).is_one()) {
                reduced = trial;
            } else {
                break;
            }
        }
    }
    return reduced;
}

vector<uint64_t> compute_orders_for_m(int m, const vector<uint32_t>& primes) {
    Poly f{1, 0};
    poly_toggle_bit(f, m);

    vector<Poly> factors = berlekamp_factor(f);
    vector<uint64_t> orders;
    orders.reserve(factors.size());
    Poly x_poly{2, 0};
    for (const Poly& g : factors) {
        Poly x_pow = poly_pow_mod(x_poly, static_cast<uint64_t>(m - 1), g);
        Poly lambda = x_poly ^ x_pow;
        lambda = poly_mod(lambda, g);
        orders.push_back(order_in_factor(lambda, g, primes));
    }
    return orders;
}

uint64_t lcm_u64(uint64_t a, uint64_t b) {
    return a / std::gcd(a, b) * b;
}

unordered_set<uint64_t> periods_for_n(int n,
                                      const vector<uint64_t>& orders,
                                      int k) {
    unordered_set<uint64_t> periods;
    periods.insert(1);
    for (uint64_t ord : orders) {
        vector<uint64_t> opts;
        opts.reserve(static_cast<size_t>(k + 2));
        opts.push_back(1);
        if (ord == 0) {
            // lambda = 0 contributes only period 1.
        } else if (ord == 1) {
            for (int s = 1; s <= k; ++s) opts.push_back(1ULL << s);
        } else {
            for (int s = 0; s <= k; ++s) opts.push_back(ord << s);
        }

        unordered_set<uint64_t> next;
        next.reserve(periods.size() * opts.size());
        for (uint64_t p : periods) {
            for (uint64_t o : opts) {
                next.insert(lcm_u64(p, o));
            }
        }
        periods.swap(next);
    }
    return periods;
}

uint64_t compute_S(int N,
                   unordered_map<int, vector<uint64_t>>& order_cache,
                   const vector<uint32_t>& primes) {
    unordered_set<uint64_t> global_periods;
    for (int n = 3; n <= N; ++n) {
        int m = n;
        int k = 0;
        while ((m & 1) == 0) {
            m >>= 1;
            ++k;
        }
        auto it = order_cache.find(m);
        if (it == order_cache.end()) {
            auto orders = compute_orders_for_m(m, primes);
            it = order_cache.emplace(m, std::move(orders)).first;
        }
        unordered_set<uint64_t> per = periods_for_n(n, it->second, k);
        global_periods.insert(per.begin(), per.end());
    }
    uint64_t sum = 0;
    for (uint64_t v : global_periods) sum += v;
    return sum;
}

bool check_periods(int n,
                   const vector<uint64_t>& expected,
                   unordered_map<int, vector<uint64_t>>& order_cache,
                   const vector<uint32_t>& primes) {
    int m = n;
    int k = 0;
    while ((m & 1) == 0) {
        m >>= 1;
        ++k;
    }
    auto it = order_cache.find(m);
    if (it == order_cache.end()) {
        auto orders = compute_orders_for_m(m, primes);
        it = order_cache.emplace(m, std::move(orders)).first;
    }
    unordered_set<uint64_t> per = periods_for_n(n, it->second, k);
    vector<uint64_t> got(per.begin(), per.end());
    vector<uint64_t> exp = expected;
    sort(got.begin(), got.end());
    sort(exp.begin(), exp.end());
    return got == exp;
}

}  // namespace

int main() {
    const int N = 100;
    vector<uint32_t> primes = sieve_primes(5'000'000);

    unordered_map<int, vector<uint64_t>> order_cache;

    if (!check_periods(5, {1, 3}, order_cache, primes) ||
        !check_periods(6, {1, 2}, order_cache, primes)) {
        cerr << "Validation failure: period sets for n=5 or n=6 mismatch\n";
        return 1;
    }

    uint64_t s30 = compute_S(30, order_cache, primes);
    if (s30 != 20381ULL) {
        cerr << "Validation failure: S(30) mismatch\n";
        return 1;
    }

    uint64_t answer = compute_S(N, order_cache, primes);
    cout << answer << '\n';
    return 0;
}

Python

from math import gcd

def solve():
    N = 100

    def sieve(n):
        ip=[True]*(n+1); ip[0]=ip[1]=False
        for i in range(2,int(n**0.5)+1):
            if ip[i]:
                for j in range(i*i,n+1,i): ip[j]=False
        return [i for i in range(2,n+1) if ip[i]]

    primes = sieve(5000000)

    def factorize(n):
        f=[]
        for p in primes:
            if p*p>n: break
            if n%p==0:
                f.append(p)
                while n%p==0: n//=p
        if n>1: f.append(n)
        return f

    # GF(2) polynomial operations using Python ints (unlimited precision)
    def pdeg(p):
        if p==0: return -1
        return p.bit_length()-1

    def pmod(a,m):
        dm=pdeg(m)
        if dm<0: return a
        while True:
            da=pdeg(a)
            if da<dm: break
            a^=m<<(da-dm)
        return a

    def pdiv(a,m):
        q=0; dm=pdeg(m)
        while True:
            da=pdeg(a)
            if da<dm: break
            s=da-dm; q|=1<<s; a^=m<<s
        return q

    def pgcd(a,b):
        while b: a,b=b,pmod(a,b)
        return a

    def pmulmod(a,b,m):
        a=pmod(a,m); b=pmod(b,m)
        r=0
        while b:
            if b&1: r^=a
            b>>=1; a<<=1
            if pdeg(a)==pdeg(m): a^=m
        return r

    def ppowmod(base,exp,m):
        r=1
        while exp>0:
            if exp&1: r=pmulmod(r,base,m)
            base=pmulmod(base,base,m); exp>>=1
        return r

    def berlekamp(f):
        m=pdeg(f)
        if m<=1: return [f]
        # Build Q matrix (x^(2j) mod f for j=0..m-1)
        rows=[0]*m; power=1; x2=pmod(4,f)  # x^2
        # Actually x=2 in poly representation: bit 1 set = x^1
        # x^2 = 4 = bit 2
        power=1  # x^0
        for col in range(m):
            for row in range(m):
                if (power>>row)&1: rows[row]^=1<<col
            power=pmulmod(power,x2,f)
        # Q-I
        for i in range(m): rows[i]^=1<<i
        # Row reduce
        pivot=[]; rank=0
        for col in range(m):
            sel=-1
            for r in range(rank,m):
                if (rows[r]>>col)&1: sel=r; break
            if sel<0: continue
            rows[rank],rows[sel]=rows[sel],rows[rank]
            pivot.append(col)
            for r in range(m):
                if r!=rank and (rows[r]>>col)&1: rows[r]^=rows[rank]
            rank+=1
        is_pivot=set(pivot)
        basis=[]
        for col in range(m):
            if col in is_pivot: continue
            vec=1<<col
            for i in range(rank):
                if (rows[i]>>col)&1: vec|=1<<pivot[i]
            basis.append(vec)
        factors=[f]
        for b in basis:
            if b==0 or b==1: continue
            nxt=[]
            for h in factors:
                if pdeg(h)<=1: nxt.append(h); continue
                g=pgcd(h,b)
                if g!=1 and g!=h:
                    nxt.append(g); nxt.append(pdiv(h,g)); continue
                b1=b^1; g=pgcd(h,b1)
                if g!=1 and g!=h:
                    nxt.append(g); nxt.append(pdiv(h,g))
                else: nxt.append(h)
            factors=nxt
        return factors

    def order_in_factor(lam,mod):
        if lam==0: return 0
        if lam==1: return 1
        deg=pdeg(mod); cur=lam; e=0
        while True:
            cur=pmulmod(cur,cur,mod); e+=1
            if cur==lam or e>deg: break
        if e>deg or e>63: return 0
        order=(1<<e)-1; reduced=order
        for p in factorize(order):
            while reduced%p==0:
                trial=reduced//p
                if ppowmod(lam,trial,mod)==1: reduced=trial
                else: break
        return reduced

    def orders_for_m(m):
        f=1|(1<<m)  # x^m + 1
        factors=berlekamp(f)
        orders=[]
        x=2  # x polynomial
        for g in factors:
            xpow=ppowmod(x,m-1,g)
            lam=pmod(x^xpow,g)
            orders.append(order_in_factor(lam,g))
        return orders

    def lcm(a,b): return a//gcd(a,b)*b

    def periods_for_n(n,orders,k):
        periods={1}
        for ord_val in orders:
            opts=[1]
            if ord_val==0: pass
            elif ord_val==1:
                for s in range(1,k+1): opts.append(1<<s)
            else:
                for s in range(k+1): opts.append(ord_val<<s)
            nxt=set()
            for p in periods:
                for o in opts: nxt.add(lcm(p,o))
            periods=nxt
        return periods

    cache={}
    global_periods=set()
    for n in range(3,N+1):
        m=n; k=0
        while m%2==0: m>>=1; k+=1
        if m not in cache: cache[m]=orders_for_m(m)
        per=periods_for_n(n,cache[m],k)
        global_periods|=per

    return str(sum(global_periods))

if __name__=='__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;

public class Euler843 {

    static class Poly {
        long lo, hi;

        Poly() {
            lo = 0;
            hi = 0;
        }

        Poly(long lo, long hi) {
            this.lo = lo;
            this.hi = hi;
        }

        boolean isZero() {
            return lo == 0 && hi == 0;
        }

        boolean isOne() {
            return lo == 1 && hi == 0;
        }

        Poly copy() {
            return new Poly(lo, hi);
        }
    }

    static Poly polyXor(Poly a, Poly b) {
        return new Poly(a.lo ^ b.lo, a.hi ^ b.hi);
    }

    static boolean polyEqual(Poly a, Poly b) {
        return a.lo == b.lo && a.hi == b.hi;
    }

    static int polyDegree(Poly p) {
        if (p.hi != 0)
            return 127 - Long.numberOfLeadingZeros(p.hi);
        if (p.lo != 0)
            return 63 - Long.numberOfLeadingZeros(p.lo);
        return -1;
    }

    static boolean polyGetBit(Poly p, int pos) {
        if (pos < 64)
            return ((p.lo >>> pos) & 1) != 0;
        return ((p.hi >>> (pos - 64)) & 1) != 0;
    }

    static void polyToggleBit(Poly p, int pos) {
        if (pos < 64)
            p.lo ^= (1L << pos);
        else
            p.hi ^= (1L << (pos - 64));
    }

    static Poly polyShiftLeft(Poly p, int k) {
        if (k == 0)
            return p.copy();
        if (k >= 64) {
            return new Poly(0, p.lo << (k - 64));
        }
        long hi = (p.hi << k) | (p.lo >>> (64 - k));
        long lo = p.lo << k;
        return new Poly(lo, hi);
    }

    static Poly polyShiftRight1(Poly p) {
        long lo = (p.lo >>> 1) | (p.hi << 63);
        long hi = p.hi >>> 1;
        return new Poly(lo, hi);
    }

    static Poly polyShiftLeft1(Poly p) {
        long hi = (p.hi << 1) | (p.lo >>> 63);
        long lo = p.lo << 1;
        return new Poly(lo, hi);
    }

    static Poly polyMod(Poly a, Poly mod) {
        int degMod = polyDegree(mod);
        if (degMod < 0)
            return a.copy();
        Poly aCpy = a.copy();
        while (true) {
            int degA = polyDegree(aCpy);
            if (degA < degMod)
                break;
            int shift = degA - degMod;
            aCpy = polyXor(aCpy, polyShiftLeft(mod, shift));
        }
        return aCpy;
    }

    static Poly polyDiv(Poly a, Poly mod) {
        Poly q = new Poly();
        int degMod = polyDegree(mod);
        Poly aCpy = a.copy();
        while (true) {
            int degA = polyDegree(aCpy);
            if (degA < degMod)
                break;
            int shift = degA - degMod;
            polyToggleBit(q, shift);
            aCpy = polyXor(aCpy, polyShiftLeft(mod, shift));
        }
        return q;
    }

    static Poly polyGcd(Poly a, Poly b) {
        Poly aCpy = a.copy();
        Poly bCpy = b.copy();
        while (!bCpy.isZero()) {
            Poly r = polyMod(aCpy, bCpy);
            aCpy = bCpy;
            bCpy = r;
        }
        return aCpy;
    }

    static Poly polyMulMod(Poly a, Poly b, Poly mod) {
        if (a.isZero() || b.isZero())
            return new Poly();
        Poly aCpy = polyMod(a, mod);
        Poly bCpy = polyMod(b, mod);
        int degMod = polyDegree(mod);
        Poly res = new Poly();

        while (!bCpy.isZero()) {
            if ((bCpy.lo & 1) != 0)
                res = polyXor(res, aCpy);
            bCpy = polyShiftRight1(bCpy);
            aCpy = polyShiftLeft1(aCpy);
            if (polyGetBit(aCpy, degMod))
                aCpy = polyXor(aCpy, mod);
        }
        return res;
    }

    static Poly polyPowMod(Poly base, long exp, Poly mod) {
        Poly res = new Poly(1, 0);
        Poly baseCpy = base.copy();
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = polyMulMod(res, baseCpy, mod);
            baseCpy = polyMulMod(baseCpy, baseCpy, mod);
            exp >>>= 1;
        }
        return res;
    }

    static ArrayList<Integer> sievePrimes(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        java.util.Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i * i <= limit; ++i) {
            if (!isPrime[i])
                continue;
            for (int j = i * i; j <= limit; j += i)
                isPrime[j] = false;
        }
        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; ++i) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static ArrayList<Long> factorizeUint64(long n, ArrayList<Integer> primes) {
        ArrayList<Long> factors = new ArrayList<>();
        long temp = n;
        for (int p : primes) {
            long pp = p;
            if (pp * pp > temp)
                break;
            if (temp % pp == 0) {
                factors.add(pp);
                while (temp % pp == 0)
                    temp /= pp;
            }
        }
        if (temp > 1)
            factors.add(temp);
        return factors;
    }

    static ArrayList<Poly> berlekampFactor(Poly f) {
        int m = polyDegree(f);
        ArrayList<Poly> resf = new ArrayList<>();
        if (m <= 1) {
            resf.add(f.copy());
            return resf;
        }

        Poly[] rows = new Poly[m];
        for (int i = 0; i < m; i++)
            rows[i] = new Poly();

        Poly xPoly = new Poly(2, 0);
        Poly x2 = polyMulMod(xPoly, xPoly, f);
        Poly power = new Poly(1, 0);

        for (int col = 0; col < m; ++col) {
            for (int row = 0; row < m; ++row) {
                if (polyGetBit(power, row)) {
                    polyToggleBit(rows[row], col);
                }
            }
            power = polyMulMod(power, x2, f);
        }
        for (int i = 0; i < m; ++i)
            polyToggleBit(rows[i], i);

        ArrayList<Integer> pivotCol = new ArrayList<>();
        int rank = 0;
        for (int col = 0; col < m; ++col) {
            int sel = -1;
            for (int r = rank; r < m; ++r) {
                if (polyGetBit(rows[r], col)) {
                    sel = r;
                    break;
                }
            }
            if (sel == -1)
                continue;
            Poly tmp = rows[rank];
            rows[rank] = rows[sel];
            rows[sel] = tmp;
            pivotCol.add(col);
            for (int r = 0; r < m; ++r) {
                if (r != rank && polyGetBit(rows[r], col)) {
                    rows[r] = polyXor(rows[r], rows[rank]);
                }
            }
            ++rank;
        }

        boolean[] isPivot = new boolean[m];
        for (int i = 0; i < rank; ++i)
            isPivot[pivotCol.get(i)] = true;

        ArrayList<Poly> basis = new ArrayList<>();
        for (int col = 0; col < m; ++col) {
            if (isPivot[col])
                continue;
            Poly vec = new Poly();
            polyToggleBit(vec, col);
            for (int i = 0; i < rank; ++i) {
                int p = pivotCol.get(i);
                if (polyGetBit(rows[i], col))
                    polyToggleBit(vec, p);
            }
            basis.add(vec);
        }

        ArrayList<Poly> factors = new ArrayList<>();
        factors.add(f.copy());

        for (Poly b : basis) {
            if (b.isZero() || b.isOne())
                continue;
            ArrayList<Poly> nextFactors = new ArrayList<>();
            for (Poly h : factors) {
                if (polyDegree(h) <= 1) {
                    nextFactors.add(h);
                    continue;
                }
                Poly g = polyGcd(h, b);
                if (!g.isOne() && !polyEqual(g, h)) {
                    Poly q = polyDiv(h, g);
                    nextFactors.add(g);
                    nextFactors.add(q);
                    continue;
                }
                Poly b1 = b.copy();
                polyToggleBit(b1, 0);
                g = polyGcd(h, b1);
                if (!g.isOne() && !polyEqual(g, h)) {
                    Poly q = polyDiv(h, g);
                    nextFactors.add(g);
                    nextFactors.add(q);
                } else {
                    nextFactors.add(h);
                }
            }
            factors = nextFactors;
        }

        return factors;
    }

    static long orderInFactor(Poly lambda, Poly mod, ArrayList<Integer> primes) {
        if (lambda.isZero())
            return 0;
        if (lambda.isOne())
            return 1;

        int deg = polyDegree(mod);
        Poly cur = lambda.copy();
        int e = 0;
        do {
            cur = polyMulMod(cur, cur, mod);
            ++e;
        } while (!polyEqual(cur, lambda) && e <= deg);

        long order = (e >= 64) ? -1L : ((1L << e) - 1L);
        long reduced = order;
        ArrayList<Long> factors = factorizeUint64(order, primes);

        for (long p : factors) {
            while (reduced % p == 0) {
                long trial = reduced / p;
                if (polyPowMod(lambda, trial, mod).isOne()) {
                    reduced = trial;
                } else {
                    break;
                }
            }
        }
        return reduced;
    }

    static ArrayList<Long> computeOrdersForM(int m, ArrayList<Integer> primes) {
        Poly f = new Poly(1, 0);
        polyToggleBit(f, m);

        ArrayList<Poly> factors = berlekampFactor(f);
        ArrayList<Long> orders = new ArrayList<>();
        Poly xPoly = new Poly(2, 0);

        for (Poly g : factors) {
            Poly xPow = polyPowMod(xPoly, m - 1, g);
            Poly lambda = polyXor(xPoly, xPow);
            lambda = polyMod(lambda, g);
            orders.add(orderInFactor(lambda, g, primes));
        }
        return orders;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    static long lcm(long a, long b) {
        return (a / gcd(a, b)) * b;
    }

    static HashSet<Long> periodsForN(int n, ArrayList<Long> orders, int k) {
        HashSet<Long> periods = new HashSet<>();
        periods.add(1L);

        for (long ord : orders) {
            ArrayList<Long> opts = new ArrayList<>();
            opts.add(1L);
            if (ord == 0) {
            } else if (ord == 1) {
                for (int s = 1; s <= k; ++s)
                    opts.add(1L << s);
            } else {
                for (int s = 0; s <= k; ++s)
                    opts.add(ord << s);
            }

            HashSet<Long> nextPeriods = new HashSet<>();
            for (long p : periods) {
                for (long o : opts) {
                    nextPeriods.add(lcm(p, o));
                }
            }
            periods = nextPeriods;
        }
        return periods;
    }

    static long computeS(int N, ArrayList<Integer> primes) {
        HashMap<Integer, ArrayList<Long>> orderCache = new HashMap<>();
        HashSet<Long> globalPeriods = new HashSet<>();

        for (int n = 3; n <= N; ++n) {
            int m = n;
            int k = 0;
            while ((m & 1) == 0) {
                m >>= 1;
                ++k;
            }
            if (!orderCache.containsKey(m)) {
                orderCache.put(m, computeOrdersForM(m, primes));
            }
            HashSet<Long> per = periodsForN(n, orderCache.get(m), k);
            globalPeriods.addAll(per);
        }

        long sum = 0;
        for (long v : globalPeriods)
            sum += v;
        return sum;
    }

    public static String solve() {
        ArrayList<Integer> primes = sievePrimes(5000000);
        return Long.toString(computeS(100, primes));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}