Problem 913: Row-major vs Column-major

Project Euler Problem 913 Solution

EulerSolve provides an optimized solution for Project Euler Problem 913, Row-major vs Column-major, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each pair $2\le n\le m\le 100$, the problem considers a rectangle with side lengths $n^4$ and $m^4$. Its cells are listed once in row-major order and once in column-major order, so the reindexing from one list to the other is a permutation of all $(nm)^4$ positions. For each pair, we must find the minimum number of swaps needed to realize that permutation, and then sum those values over all admissible pairs. A direct cycle walk on the full permutation would be hopeless, because the rectangle for one pair already has $Q=(nm)^4$ cells. The implementations succeed by turning the transpose permutation into modular multiplication and then counting cycles through divisors of $Q-1$. Mathematical Approach Let the rectangle have $r=m^4$ rows and $c=n^4$ columns. Define $$Q=rc=(nm)^4,\qquad L=Q-1,\qquad \alpha=r=m^4.$$ The entire task becomes a problem about the cycle structure of one modular permutation. From Matrix Indices to a Modular Permutation Take a cell in row $i$ and column $j$, with $0\le i\lt r$ and $0\le j\lt c$....

Detailed mathematical approach

Problem Summary

For each pair $2\le n\le m\le 100$, the problem considers a rectangle with side lengths $n^4$ and $m^4$. Its cells are listed once in row-major order and once in column-major order, so the reindexing from one list to the other is a permutation of all $(nm)^4$ positions. For each pair, we must find the minimum number of swaps needed to realize that permutation, and then sum those values over all admissible pairs.

A direct cycle walk on the full permutation would be hopeless, because the rectangle for one pair already has $Q=(nm)^4$ cells. The implementations succeed by turning the transpose permutation into modular multiplication and then counting cycles through divisors of $Q-1$.

Mathematical Approach

Let the rectangle have $r=m^4$ rows and $c=n^4$ columns. Define

$$Q=rc=(nm)^4,\qquad L=Q-1,\qquad \alpha=r=m^4.$$

The entire task becomes a problem about the cycle structure of one modular permutation.

From Matrix Indices to a Modular Permutation

Take a cell in row $i$ and column $j$, with $0\le i\lt r$ and $0\le j\lt c$. Its row-major index is

$$x=ic+j,$$

while its column-major index is

$$y=jr+i.$$

Now compute $rx$:

$$rx=r(ic+j)=irc+rj=iQ+rj\equiv i+rj=y \pmod{Q-1}.$$

So for every index except the last one, the reindexing permutation satisfies

$$\pi(x)\equiv \alpha x \pmod L,\qquad 0\le x\lt L.$$

The missing index $x=L=Q-1$ is an extra fixed point. Because

$$L=(nm)^4-1\equiv -1 \pmod m,$$

we have $\gcd(\alpha,L)=1$, so multiplication by $\alpha$ genuinely permutes the residues modulo $L$.

If a permutation on $Q$ elements has $C$ cycles, then the minimum number of swaps needed to realize it is $Q-C$, since a cycle of length $t$ takes exactly $t-1$ swaps. The problem is therefore reduced to counting cycles of $\pi$.

Layers Indexed by Additive Order

For a residue $x\in \mathbb{Z}/L\mathbb{Z}$, define

$$d=\frac{L}{\gcd(x,L)}.$$

This number $d$ is the additive order of $x$ in the cyclic group $\mathbb{Z}/L\mathbb{Z}$. Multiplication by $\alpha$ preserves $\gcd(x,L)$, because $\alpha$ is invertible modulo $L$. Therefore it also preserves $d$. The permutation splits into invariant layers, one layer for each divisor $d\mid L$.

Exactly $\varphi(d)$ residues have additive order $d$: they are the generators of the unique subgroup of size $d$ inside the cyclic group $\mathbb{Z}/L\mathbb{Z}$.

Why the Cycle Length Becomes a Multiplicative Order

Write

$$g=\gcd(x,L),\qquad x=gu,\qquad L=gd,$$

so that $\gcd(u,d)=1$. Then the condition that $x$ returns to itself after $t$ steps is

$$\alpha^t x\equiv x \pmod L.$$

Substituting $x=gu$ and $L=gd$ gives

$$gd\mid g(\alpha^t-1)u\iff d\mid(\alpha^t-1)u.$$

Because $u$ is invertible modulo $d$, this is equivalent to

$$\alpha^t\equiv 1 \pmod d.$$

So every orbit in the layer indexed by $d$ has length

$$\operatorname{ord}_d(\alpha),$$

the multiplicative order of $\alpha$ modulo $d$. Since the layer contains $\varphi(d)$ elements, its number of cycles is

$$\frac{\varphi(d)}{\operatorname{ord}_d(\alpha)}.$$

The Divisor Sum for the Total Cycle Count

Summing those contributions over all divisors of $L$ gives the number of cycles inside the modular part:

$$C_L=\sum_{d\mid L}\frac{\varphi(d)}{\operatorname{ord}_d(\alpha)},$$

with the natural convention $\operatorname{ord}_1(\alpha)=1$. This already counts the fixed point $x=0$. The index $x=L=Q-1$ contributes one more fixed point outside the modular residue set, so the full cycle count is

$$C=C_L+1.$$

Therefore the swap count for the pair $(n,m)$ is

$$\boxed{\text{swaps}=Q-1-\sum_{d\mid L}\frac{\varphi(d)}{\operatorname{ord}_d(\alpha)}.}$$

Worked Example: A $4\times 3$ Rectangle

The actual problem uses fourth powers, but the same mechanism is already visible in a small ordinary rectangle with $r=4$ and $c=3$. Then

$$Q=12,\qquad L=11,\qquad \alpha=4,$$

so on $0\le x\lt 11$ we have

$$\pi(x)\equiv 4x \pmod{11}.$$

The only divisors of $11$ are $1$ and $11$, hence

$$C_{11}=1+\frac{\varphi(11)}{\operatorname{ord}_{11}(4)}=1+\frac{10}{5}=3.$$

Indeed, the residues split as

$$0,$$

$$1\to 4\to 5\to 9\to 3\to 1,$$

$$2\to 8\to 10\to 7\to 6\to 2,$$

and the twelfth position $x=11$ is the extra fixed point. So the total cycle count is $4$, and the minimum number of swaps is

$$12-4=8.$$

This is the same small check used by the implementations to confirm the formula.

Factoring $s^4-1$ Without Touching the Matrix

Set $s=nm$. Then

$$L=s^4-1=(s-1)(s+1)(s^2+1).$$

Because $s\le 10^4$, the three factors $s-1$, $s+1$, and $s^2+1$ are moderate-size integers. Factoring those three numbers is enough to recover the complete prime factorization of $L$.

For each prime power $p^k\mid L$, the order satisfies

$$\operatorname{ord}_{p^k}(\alpha)\mid \varphi(p^k)=p^{k-1}(p-1).$$

The exact order is obtained by starting from $\varphi(p^k)$ and repeatedly removing prime factors whenever the modular test

$$\alpha^{\,\varphi(p^k)/q}\equiv 1 \pmod{p^k}$$

still holds.

$$d=\prod_i p_i^{e_i},$$

then multiplicativity and the Chinese remainder theorem give

$$\varphi(d)=\prod_i \varphi(p_i^{e_i}),\qquad \operatorname{ord}_d(\alpha)=\operatorname{lcm}_i\operatorname{ord}_{p_i^{e_i}}(\alpha).$$

So a depth-first enumeration of the exponent choices $0\le e_i\le v_{p_i}(L)$ visits every divisor $d\mid L$ and accumulates its cycle contribution exactly, without ever traversing the $Q$ entries of the rectangle.

How the Code Works

The C++, Python, and Java implementations follow the same number-theoretic pipeline.

Reusable Factorization Data

A prime sieve is built once. For a given pair $(n,m)$, the implementation forms $s=nm$, factors $s-1$, $s+1$, and $s^2+1$, and merges those exponents into the factorization of $L=s^4-1$. Because different pairs can share the same product $s$, this factorization is cached and reused. Factorizations of numbers of the form $p-1$ are also cached, because they are needed repeatedly when reducing $\varphi(p^k)$ to the exact order $\operatorname{ord}_{p^k}(\alpha)$.

Per-Pair Evaluation

For each prime power $p^k\mid L$, the implementation precomputes $\varphi(p^k)$ and $\operatorname{ord}_{p^k}(\alpha)$. It then runs a divisor DFS: each recursion level chooses one exponent for one prime, multiplies the running totient contribution, updates the running order with an $\operatorname{lcm}$, and descends to the next prime. At every leaf, the divisor contribution

$$\frac{\varphi(d)}{\operatorname{ord}_d(\alpha)}$$

is added to the cycle total $C_L$. The returned swap count is therefore $Q-1-C_L$.

The outer loops sum this value over all $2\le n\le m\le 100$. The implementations also compare the divisor formula against explicit cycle decompositions on small cases; for example, the ordinary $3\times 4$ rectangle gives $8$ swaps exactly, matching the theory above.

Complexity Analysis

For one pair $(n,m)$, the algorithm never iterates over all $Q=(nm)^4$ positions. Its work is concentrated in three places: factoring $s-1$, $s+1$, and $s^2+1$; computing multiplicative orders on the relevant prime powers; and enumerating the $\tau(L)$ divisors of $L$ with a DFS over prime exponents. In practice that is far smaller than a brute-force cycle walk on the full transpose permutation.

Memory usage is modest. The main data are the prime list, the caches for repeated factorizations, and the prime-power tables for the current pair. The method is efficient precisely because it manipulates the factorization of $L=(nm)^4-1$, not the $(nm)^4$ matrix entries themselves.

Footnotes and References

Problem page: https://projecteuler.net/problem=913
Row-major and column-major order: Wikipedia - Row- and column-major order
In-place matrix transposition: Wikipedia - In-place matrix transposition
Euler's totient function: Wikipedia - Euler's totient function
Multiplicative order: Wikipedia - Multiplicative order
Chinese remainder theorem: Wikipedia - Chinese remainder theorem

Problem 913 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <map>
#include <numeric>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct PrimeData {
    int exp = 0;
    std::vector<u64> phi;
    std::vector<u64> ord;
};

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

std::vector<std::pair<u64, int>> factorize_u64(u64 x, const std::vector<int>& primes) {
    std::vector<std::pair<u64, int>> out;
    u64 n = x;
    for (int p : primes) {
        if (1ULL * p * p > n) {
            break;
        }
        if (n % static_cast<u64>(p) != 0) {
            continue;
        }
        int e = 0;
        while (n % static_cast<u64>(p) == 0) {
            n /= static_cast<u64>(p);
            ++e;
        }
        out.push_back({static_cast<u64>(p), e});
    }
    if (n > 1) {
        out.push_back({n, 1});
    }
    return out;
}

void add_factorization(std::map<u64, int>& mp, u64 x, const std::vector<int>& primes) {
    if (x <= 1) {
        return;
    }
    const auto f = factorize_u64(x, primes);
    for (const auto& [p, e] : f) {
        mp[p] += e;
    }
}

u64 mod_pow(u64 a, u64 e, u64 mod) {
    if (mod == 1) {
        return 0;
    }
    u64 base = a % mod;
    u64 exp = e;
    u64 res = 1 % mod;
    while (exp > 0) {
        if (exp & 1ULL) {
            res = static_cast<u64>((u128)res * base % mod);
        }
        base = static_cast<u64>((u128)base * base % mod);
        exp >>= 1ULL;
    }
    return res;
}

u64 ipow_u64(u64 a, int e) {
    u64 r = 1;
    for (int i = 0; i < e; ++i) {
        r *= a;
    }
    return r;
}

u64 lcm_u64(u64 a, u64 b) {
    if (a == 0 || b == 0) {
        return 0;
    }
    const u64 g = std::gcd(a, b);
    return static_cast<u64>((u128)(a / g) * b);
}

const std::vector<std::pair<u64, int>>& factor_cached(
    u64 x,
    const std::vector<int>& primes,
    std::unordered_map<u64, std::vector<std::pair<u64, int>>>& cache) {
    auto it = cache.find(x);
    if (it != cache.end()) {
        return it->second;
    }
    auto [jt, _] = cache.emplace(x, factorize_u64(x, primes));
    return jt->second;
}

u64 order_prime_power(
    u64 base,
    u64 p,
    int k,
    const std::vector<int>& primes,
    std::unordered_map<u64, std::vector<std::pair<u64, int>>>& factor_cache) {
    const u64 mod = ipow_u64(p, k);
    const u64 phi = ipow_u64(p, k - 1) * (p - 1);

    std::vector<std::pair<u64, int>> fac = factor_cached(p - 1, primes, factor_cache);
    if (k > 1) {
        bool found = false;
        for (auto& [q, e] : fac) {
            if (q == p) {
                e += k - 1;
                found = true;
                break;
            }
        }
        if (!found) {
            fac.push_back({p, k - 1});
        }
    }

    u64 ord = phi;
    for (const auto& [q, e] : fac) {
        for (int i = 0; i < e; ++i) {
            if (ord % q != 0) {
                break;
            }
            const u64 cand = ord / q;
            if (mod_pow(base, cand, mod) == 1) {
                ord = cand;
            } else {
                break;
            }
        }
    }
    return ord;
}

u64 swaps_fourth_power_dims(
    int n,
    int m,
    const std::vector<int>& primes,
    std::unordered_map<u64, std::vector<std::pair<u64, int>>>& factor_cache,
    std::unordered_map<int, std::vector<std::pair<u64, int>>>& t_factor_cache) {
    const u64 uu = static_cast<u64>(n) * static_cast<u64>(m);
    const u64 N = ipow_u64(uu, 4);
    const u64 M = N - 1;
    const u64 g = ipow_u64(static_cast<u64>(m), 4);

    std::vector<std::pair<u64, int>> factors;
    auto it_t = t_factor_cache.find(static_cast<int>(uu));
    if (it_t != t_factor_cache.end()) {
        factors = it_t->second;
    } else {
        std::map<u64, int> acc;
        add_factorization(acc, uu - 1, primes);
        add_factorization(acc, uu + 1, primes);
        add_factorization(acc, uu * uu + 1, primes);

        factors.reserve(acc.size());
        for (const auto& [p, e] : acc) {
            factors.push_back({p, e});
        }
        t_factor_cache.emplace(static_cast<int>(uu), factors);
    }

    std::vector<PrimeData> data;
    data.reserve(factors.size());

    for (const auto& [p, e] : factors) {
        PrimeData pd;
        pd.exp = e;
        pd.phi.assign(static_cast<std::size_t>(e + 1), 1);
        pd.ord.assign(static_cast<std::size_t>(e + 1), 1);

        for (int k = 1; k <= e; ++k) {
            pd.phi[static_cast<std::size_t>(k)] = ipow_u64(p, k - 1) * (p - 1);
            pd.ord[static_cast<std::size_t>(k)] =
                order_prime_power(g, p, k, primes, factor_cache);
        }

        data.push_back(std::move(pd));
    }

    u64 cycles_mod_set = 0;

    std::function<void(std::size_t, u64, u64)> dfs = [&](std::size_t idx, u64 phi_cur, u64 ord_cur) {
        if (idx == data.size()) {
            cycles_mod_set += phi_cur / ord_cur;
            return;
        }

        const PrimeData& pd = data[idx];
        for (int k = 0; k <= pd.exp; ++k) {
            const u64 phi_next = static_cast<u64>((u128)phi_cur * pd.phi[static_cast<std::size_t>(k)]);
            const u64 ord_next = lcm_u64(ord_cur, pd.ord[static_cast<std::size_t>(k)]);
            dfs(idx + 1, phi_next, ord_next);
        }
    };

    dfs(0, 1, 1);

    return N - 1 - cycles_mod_set;
}

u64 brute_swaps(u64 n, u64 m) {
    const u64 N = n * m;
    std::vector<char> vis(static_cast<std::size_t>(N), 0);

    auto next = [&](u64 x) -> u64 {
        if (x == N - 1) {
            return N - 1;
        }
        return static_cast<u64>((u128)x * m % (N - 1));
    };

    u64 cycles = 0;
    for (u64 i = 0; i < N; ++i) {
        if (vis[static_cast<std::size_t>(i)]) {
            continue;
        }
        ++cycles;
        u64 x = i;
        while (!vis[static_cast<std::size_t>(x)]) {
            vis[static_cast<std::size_t>(x)] = 1;
            x = next(x);
        }
    }

    return N - cycles;
}

void print_u128(u128 x) {
    if (x == 0) {
        std::cout << 0;
        return;
    }
    std::string s;
    u128 v = x;
    while (v > 0) {
        const int digit = static_cast<int>(v % 10);
        s.push_back(static_cast<char>('0' + digit));
        v /= 10;
    }
    std::reverse(s.begin(), s.end());
    std::cout << s;
}

void validate(
    const std::vector<int>& primes,
    std::unordered_map<u64, std::vector<std::pair<u64, int>>>& factor_cache,
    std::unordered_map<int, std::vector<std::pair<u64, int>>>& t_factor_cache) {
    assert(brute_swaps(3, 4) == 8);

    u64 sum_small = 0;
    for (int n = 2; n <= 100; ++n) {
        for (int m = n; m <= 100; ++m) {
            sum_small += brute_swaps(static_cast<u64>(n), static_cast<u64>(m));
        }
    }
    assert(sum_small == 12'578'833ULL);

    for (int n = 2; n <= 5; ++n) {
        for (int m = n; m <= 5; ++m) {
            const u64 nn = ipow_u64(static_cast<u64>(n), 4);
            const u64 mm = ipow_u64(static_cast<u64>(m), 4);
            const u64 brute = brute_swaps(nn, mm);
            const u64 fast = swaps_fourth_power_dims(n, m, primes, factor_cache, t_factor_cache);
            assert(brute == fast);
        }
    }
}

}  // namespace

int main() {
    const std::vector<int> primes = sieve_primes(100'000);
    std::unordered_map<u64, std::vector<std::pair<u64, int>>> factor_cache;
    std::unordered_map<int, std::vector<std::pair<u64, int>>> t_factor_cache;

    validate(primes, factor_cache, t_factor_cache);

    u128 answer = 0;
    for (int n = 2; n <= 100; ++n) {
        for (int m = n; m <= 100; ++m) {
            answer += swaps_fourth_power_dims(n, m, primes, factor_cache, t_factor_cache);
        }
    }

    print_u128(answer);
    std::cout << '\n';
    return 0;
}

Python

import math
from collections import defaultdict

def sieve_primes(limit):
    is_prime = [True] * (limit + 1)
    is_prime[0] = is_prime[1] = False
    for p in range(2, int(math.sqrt(limit)) + 1):
        if is_prime[p]:
            for q in range(p * p, limit + 1, p):
                is_prime[q] = False
    return [i for i, prime in enumerate(is_prime) if prime]

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

def add_factorization(acc, n, primes):
    if n <= 1:
        return
    for p, e in factorize(n, primes):
        acc[p] += e

def order_prime_power(base, p, k, primes, factor_cache):
    mod = p ** k
    phi = (p ** (k - 1)) * (p - 1)
    
    if p - 1 not in factor_cache:
        factor_cache[p - 1] = dict(factorize(p - 1, primes))
    fac = dict(factor_cache[p - 1])
    
    if k > 1:
        fac[p] = fac.get(p, 0) + k - 1
        
    ord_val = phi
    for q, e in fac.items():
        for _ in range(e):
            if ord_val % q != 0:
                break
            cand = ord_val // q
            if pow(base, cand, mod) == 1:
                ord_val = cand
            else:
                break
    return ord_val

def swaps_fourth_power_dims(n, m, primes, factor_cache, t_factor_cache):
    uu = n * m
    N = uu ** 4
    g = m ** 4
    
    if uu not in t_factor_cache:
        acc = defaultdict(int)
        add_factorization(acc, uu - 1, primes)
        add_factorization(acc, uu + 1, primes)
        add_factorization(acc, uu * uu + 1, primes)
        t_factor_cache[uu] = acc
    factors = t_factor_cache[uu].items()
    
    data = []
    for p, e in factors:
        phi = [1] * (e + 1)
        ord_arr = [1] * (e + 1)
        for k in range(1, e + 1):
            phi[k] = (p ** (k - 1)) * (p - 1)
            ord_arr[k] = order_prime_power(g, p, k, primes, factor_cache)
        data.append({'exp': e, 'phi': phi, 'ord': ord_arr})
        
    cycles_mod_set = [0]
    
    def dfs(idx, phi_cur, ord_cur):
        if idx == len(data):
            cycles_mod_set[0] += phi_cur // ord_cur
            return
        pd = data[idx]
        for k in range(pd['exp'] + 1):
            phi_next = phi_cur * pd['phi'][k]
            ord_next = math.lcm(ord_cur, pd['ord'][k])
            dfs(idx + 1, phi_next, ord_next)
            
    dfs(0, 1, 1)
    return N - 1 - cycles_mod_set[0]

def solve():
    primes = sieve_primes(100000)
    factor_cache = {}
    t_factor_cache = {}
    
    ans = 0
    for n in range(2, 101):
        for m in range(n, 101):
            ans += swaps_fourth_power_dims(n, m, primes, factor_cache, t_factor_cache)
            
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler913 {

    static class PrimeData {
        int exp;
        long[] phi;
        long[] ord;
    }

    static List<Integer> sievePrimes(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++)
            isPrime[i] = true;
        for (int p = 2; (long) p * p <= limit; ++p) {
            if (isPrime[p]) {
                for (int q = p * p; q <= limit; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; ++i) {
            if (isPrime[i]) {
                primes.add(i);
            }
        }
        return primes;
    }

    static class Pair {
        long p;
        int e;

        Pair(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Pair> factorize(long x, List<Integer> primes) {
        List<Pair> factors = new ArrayList<>();
        long n = x;
        for (int p : primes) {
            if ((long) p * p > n)
                break;
            if (n % p == 0) {
                int e = 0;
                while (n % p == 0) {
                    n /= p;
                    e++;
                }
                factors.add(new Pair(p, e));
            }
        }
        if (n > 1) {
            factors.add(new Pair(n, 1));
        }
        return factors;
    }

    static void addFactorization(Map<Long, Integer> acc, long n, List<Integer> primes) {
        if (n <= 1)
            return;
        List<Pair> factors = factorize(n, primes);
        for (Pair pair : factors) {
            acc.put(pair.p, acc.getOrDefault(pair.p, 0) + pair.e);
        }
    }

    static long orderPrimePower(long base, long p, int k, List<Integer> primes,
            Map<Long, Map<Long, Integer>> factorCache) {
        long mod = 1;
        for (int i = 0; i < k; i++)
            mod *= p;

        long phi = 1;
        for (int i = 0; i < k - 1; i++)
            phi *= p;
        phi *= (p - 1);

        if (!factorCache.containsKey(p - 1)) {
            Map<Long, Integer> f = new HashMap<>();
            for (Pair pair : factorize(p - 1, primes)) {
                f.put(pair.p, pair.e);
            }
            factorCache.put(p - 1, f);
        }

        Map<Long, Integer> fac = new HashMap<>(factorCache.get(p - 1));
        if (k > 1) {
            fac.put(p, fac.getOrDefault(p, 0) + k - 1);
        }

        long ordVal = phi;
        BigInteger B = BigInteger.valueOf(base);
        BigInteger M = BigInteger.valueOf(mod);

        for (Map.Entry<Long, Integer> entry : fac.entrySet()) {
            long q = entry.getKey();
            int e = entry.getValue();
            for (int i = 0; i < e; ++i) {
                if (ordVal % q != 0)
                    break;
                long cand = ordVal / q;
                if (B.modPow(BigInteger.valueOf(cand), M).equals(BigInteger.ONE)) {
                    ordVal = cand;
                } else {
                    break;
                }
            }
        }
        return ordVal;
    }

    static long gcd(long a, long b) {
        return b == 0 ? a : gcd(b, a % b);
    }

    static long lcm(long a, long b) {
        if (a == 0 || b == 0)
            return 0;
        return (a / gcd(a, b)) * b;
    }

    static long cyclesModSet = 0;

    static void dfs(int idx, long phiCur, long ordCur, List<PrimeData> data) {
        if (idx == data.size()) {
            cyclesModSet += phiCur / ordCur;
            return;
        }
        PrimeData pd = data.get(idx);
        for (int k = 0; k <= pd.exp; ++k) {
            long phiNext = phiCur * pd.phi[k];
            long ordNext = lcm(ordCur, pd.ord[k]);
            dfs(idx + 1, phiNext, ordNext, data);
        }
    }

    static long swapsFourthPowerDims(int n, int m, List<Integer> primes, Map<Long, Map<Long, Integer>> factorCache,
            Map<Long, Map<Long, Integer>> tFactorCache) {
        long uu = (long) n * m;
        long N = uu * uu * uu * uu;
        long g = (long) m * m * m * m;

        if (!tFactorCache.containsKey(uu)) {
            Map<Long, Integer> acc = new HashMap<>();
            addFactorization(acc, uu - 1, primes);
            addFactorization(acc, uu + 1, primes);
            addFactorization(acc, uu * uu + 1, primes);
            tFactorCache.put(uu, acc);
        }

        Map<Long, Integer> factors = tFactorCache.get(uu);
        List<PrimeData> data = new ArrayList<>();

        for (Map.Entry<Long, Integer> entry : factors.entrySet()) {
            long p = entry.getKey();
            int e = entry.getValue();

            PrimeData pd = new PrimeData();
            pd.exp = e;
            pd.phi = new long[e + 1];
            pd.ord = new long[e + 1];
            pd.phi[0] = 1;
            pd.ord[0] = 1;

            for (int k = 1; k <= e; ++k) {
                long po = 1;
                for (int i = 0; i < k - 1; i++)
                    po *= p;
                pd.phi[k] = po * (p - 1);
                pd.ord[k] = orderPrimePower(g, p, k, primes, factorCache);
            }
            data.add(pd);
        }

        cyclesModSet = 0;
        dfs(0, 1, 1, data);

        return N - 1 - cyclesModSet;
    }

    public static String solve() {
        List<Integer> primes = sievePrimes(100000);
        Map<Long, Map<Long, Integer>> factorCache = new HashMap<>();
        Map<Long, Map<Long, Integer>> tFactorCache = new HashMap<>();

        BigInteger ans = BigInteger.ZERO;
        for (int n = 2; n <= 100; ++n) {
            for (int m = n; m <= 100; ++m) {
                long res = swapsFourthPowerDims(n, m, primes, factorCache, tFactorCache);
                ans = ans.add(BigInteger.valueOf(res));
            }
        }
        return ans.toString();
    }

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