Problem 813: XOR-Powers

View on Project Euler

Project Euler Problem 813 Solution

EulerSolve provides an optimized solution for Project Euler Problem 813, XOR-Powers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The solution views XOR-multiplication as polynomial multiplication over \(\mathbb{F}_2\). The enormous base value can be written as $$T^3+T+1,\qquad T=2^{2^{52}},$$ so the task becomes: build the polynomial $$g(x)=(x^3+x+1)^{3^8}=(x^3+x+1)^{6561}$$ inside \(\mathbb{F}_2[x]\), then evaluate \(g(T)\) modulo \(10^9+7\). This separates the carryless algebra from the final modular arithmetic and avoids manipulating the astronomical integer directly. Mathematical Approach The three implementations all use the same mathematical decomposition: first construct a polynomial with coefficients in \(\{0,1\}\), then substitute a modular value for \(x\). Step 1: Encode Bit Patterns as Polynomials If a nonnegative integer has binary expansion $$n=\sum_{k\ge 0} b_k 2^k,\qquad b_k\in\{0,1\},$$ we associate to it the polynomial $$P_n(x)=\sum_{k\ge 0} b_k x^k \in \mathbb{F}_2[x].$$ Under this encoding, adding coefficients modulo \(2\) is exactly XOR. Therefore carryless multiplication of integers is the same as ordinary polynomial multiplication in \(\mathbb{F}_2[x]\): if $$A(x)=\sum_i a_i x^i,\qquad B(x)=\sum_j b_j x^j,$$ then the product coefficients satisfy $$c_k=\left(\sum_{i+j=k} a_i b_j\right)\bmod 2=\bigoplus_{i+j=k}(a_i\land b_j).$$ That identity is the whole reason the XOR-power can be handled as a polynomial power....

Detailed mathematical approach

Problem Summary

The solution views XOR-multiplication as polynomial multiplication over \(\mathbb{F}_2\). The enormous base value can be written as

$$T^3+T+1,\qquad T=2^{2^{52}},$$

so the task becomes: build the polynomial

$$g(x)=(x^3+x+1)^{3^8}=(x^3+x+1)^{6561}$$

inside \(\mathbb{F}_2[x]\), then evaluate \(g(T)\) modulo \(10^9+7\). This separates the carryless algebra from the final modular arithmetic and avoids manipulating the astronomical integer directly.

Mathematical Approach

The three implementations all use the same mathematical decomposition: first construct a polynomial with coefficients in \(\{0,1\}\), then substitute a modular value for \(x\).

Step 1: Encode Bit Patterns as Polynomials

If a nonnegative integer has binary expansion

$$n=\sum_{k\ge 0} b_k 2^k,\qquad b_k\in\{0,1\},$$

we associate to it the polynomial

$$P_n(x)=\sum_{k\ge 0} b_k x^k \in \mathbb{F}_2[x].$$

Under this encoding, adding coefficients modulo \(2\) is exactly XOR. Therefore carryless multiplication of integers is the same as ordinary polynomial multiplication in \(\mathbb{F}_2[x]\): if

$$A(x)=\sum_i a_i x^i,\qquad B(x)=\sum_j b_j x^j,$$

then the product coefficients satisfy

$$c_k=\left(\sum_{i+j=k} a_i b_j\right)\bmod 2=\bigoplus_{i+j=k}(a_i\land b_j).$$

That identity is the whole reason the XOR-power can be handled as a polynomial power.

Step 2: Rewrite the Huge Base in Polynomial Form

Let

$$f(x)=x^3+x+1,\qquad T=2^{2^{52}}.$$

Then the target base integer is exactly

$$f(T)=T^3+T+1=2^{3\cdot 2^{52}}+2^{2^{52}}+1.$$

So instead of working with a gigantic ordinary integer, we only need the much smaller coefficient pattern \((1,1,0,1)\), i.e. the polynomial \(f(x)\). If the carryless product is denoted by \(\otimes\), then raising the integer \(f(T)\) to an XOR-power corresponds to raising the polynomial \(f(x)\) in \(\mathbb{F}_2[x]\) first and substituting \(x=T\) afterward.

Step 3: Compute the XOR-Power as a Polynomial Power

The required exponent is

$$3^8=6561,$$

so the polynomial to construct is

$$g(x)=f(x)^{6561}\in \mathbb{F}_2[x].$$

The implementations obtain this with binary exponentiation. Each multiplication is carryless, so every partial result remains a polynomial with coefficients \(0\) or \(1\). Since \(\deg f=3\), the final degree is

$$\deg g = 3\cdot 6561 = 19683.$$

That degree is large, but still tiny compared with the size of the original integer suggested by \(2^{2^{52}}\).

Step 4: Evaluate in the Modular Ring

Once \(g(x)\) has been constructed, the remaining task is numerical evaluation. Let

$$M=10^9+7,\qquad X\equiv T\equiv 2^{2^{52}}\pmod M.$$

If

$$g(x)=\sum_{k=0}^{19683} g_k x^k,\qquad g_k\in\{0,1\},$$

then the required value is

$$g(X)\equiv \sum_{k=0}^{19683} g_k X^k \pmod M.$$

This evaluation happens in \(\mathbb{Z}/M\mathbb{Z}\), not in \(\mathbb{F}_2[x]\). The switch of rings is valid because the polynomial coefficients are already fixed concrete bits. After the coefficient list is known, the same polynomial can be substituted into any target ring.

Step 5: Worked Example

A small example shows why carryless squaring behaves differently from ordinary squaring. Over \(\mathbb{F}_2\), cross terms cancel in pairs, so

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

The coefficient pattern is therefore

$$1000101_2,$$

which is decimal \(69\). In other words, the carryless square of the bit pattern \(1011_2\) is \(1000101_2\). The full problem uses exactly the same rule, only with exponent \(6561\) instead of \(2\), and with final substitution at \(X\equiv 2^{2^{52}}\pmod M\).

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First, the polynomial is stored as a large integer whose \(k\)-th bit indicates whether the coefficient of \(x^k\) is \(1\). Carryless multiplication is then performed by scanning the set bits of one operand and XORing shifted copies of the other operand, which is exactly the shift-and-add idea with XOR instead of ordinary addition.

Next, repeated squaring and conditional multiplication build \((x^3+x+1)^{6561}\). After that, ordinary modular exponentiation computes \(2^{2^{52}} \bmod 10^9+7\). Finally, the polynomial is evaluated from highest degree to lowest degree with Horner's method, using only modular multiplication by \(X\) and a possible addition of \(1\) when the current coefficient bit is present.

Complexity Analysis

Let \(E=6561\) and \(d=19683\). With the bitset-style representation used here, one carryless multiplication of degree-\(d\) polynomials costs \(O(d^2)\) bit operations in the naive worst case. Binary exponentiation uses \(O(\log E)\) such multiplications, while the final Horner evaluation costs \(O(d)\) modular operations. The memory usage is \(O(d)\) bits for the current large polynomial values. For this specific problem, those bounds are easily practical because \(d\) is only \(19683\).

Footnotes and References

  1. Project Euler, Problem 813: https://projecteuler.net/problem=813
  2. Carry-less product: Wikipedia — Carry-less product
  3. Finite field: Wikipedia — Finite field
  4. Exponentiation by squaring: Wikipedia — Exponentiation by squaring
  5. Horner's method: Wikipedia — Horner's method

Problem 813 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>

#include <boost/multiprecision/cpp_int.hpp>

using boost::multiprecision::bit_test;
using boost::multiprecision::cpp_int;
using boost::multiprecision::lsb;
using boost::multiprecision::msb;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

static constexpr u64 kMod = 1'000'000'007ULL;

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

static cpp_int clmul(cpp_int a, const cpp_int& b) {
    cpp_int res = 0;
    std::size_t shift = 0;
    while (a != 0) {
        const unsigned tz = static_cast<unsigned>(lsb(a));
        shift += tz;
        a >>= tz;
        res ^= (b << shift);
        a >>= 1;
        ++shift;
    }
    return res;
}

static cpp_int clpow(cpp_int base, u64 exp) {
    cpp_int res = 1;
    while (exp > 0) {
        if (exp & 1ULL) {
            res = clmul(res, base);
        }
        exp >>= 1ULL;
        if (exp > 0) {
            base = clmul(base, base);
        }
    }
    return res;
}

static u64 eval_poly_mod(const cpp_int& poly_bits, u64 x_mod) {
    const std::size_t deg = static_cast<std::size_t>(msb(poly_bits));
    u64 ans = 0;
    for (std::size_t i = deg + 1; i-- > 0;) {
        ans = static_cast<u64>((static_cast<u128>(ans) * x_mod) % kMod);
        if (bit_test(poly_bits, i)) {
            ++ans;
            if (ans == kMod) {
                ans = 0;
            }
        }
    }
    return ans;
}

int main() {
    const cpp_int eleven = 11;
    assert(clmul(eleven, eleven) == 69);
    assert(clpow(eleven, 2) == 69);

    const u64 three_pow_8 = 6561ULL;
    const u64 two_pow_52 = (1ULL << 52);

    const cpp_int f = (cpp_int(1) << 3) ^ (cpp_int(1) << 1) ^ cpp_int(1);
    const cpp_int g = clpow(f, three_pow_8);

    const u64 x = mod_pow(2, two_pow_52);
    std::cout << eval_poly_mod(g, x) << '\n';
    return 0;
}

Python

import sys

MOD = 1000000007

def mod_pow(base, exp):
    res = 1
    base %= MOD
    while exp > 0:
        if exp & 1:
            res = (res * base) % MOD
        base = (base * base) % MOD
        exp >>= 1
    return res

def clmul(a, b):
    res = 0
    shift = 0
    while a > 0:
        tz = (a & -a).bit_length() - 1
        shift += tz
        a >>= tz
        res ^= (b << shift)
        a >>= 1
        shift += 1
    return res

def clpow(base, exp):
    res = 1
    while exp > 0:
        if exp & 1:
            res = clmul(res, base)
        exp >>= 1
        if exp > 0:
            base = clmul(base, base)
    return res

def eval_poly_mod(poly_bits, x_mod):
    if poly_bits == 0:
        return 0
    deg = poly_bits.bit_length() - 1
    ans = 0
    for i in range(deg, -1, -1):
        ans = (ans * x_mod) % MOD
        if (poly_bits >> i) & 1:
            ans += 1
            if ans == MOD:
                ans = 0
    return ans

def solve():
    three_pow_8 = 6561
    two_pow_52 = 1 << 52
    
    f = 11  # X^3 + X + 1
    g = clpow(f, three_pow_8)
    
    x = mod_pow(2, two_pow_52)
    ans = eval_poly_mod(g, x)
    return str(ans)

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

Java

import java.math.BigInteger;

public class Euler813 {

    static final long MOD = 1000000007L;

    static long modPow(long base, long exp) {
        long res = 1;
        base %= MOD;
        while (exp > 0) {
            if ((exp & 1L) == 1L) {
                res = (res * base) % MOD;
            }
            base = (base * base) % MOD;
            exp >>= 1L;
        }
        return res;
    }

    static BigInteger clmul(BigInteger a, BigInteger b) {
        BigInteger res = BigInteger.ZERO;
        int shift = 0;
        while (a.signum() > 0) {
            int tz = a.getLowestSetBit();
            shift += tz;
            a = a.shiftRight(tz);
            res = res.xor(b.shiftLeft(shift));
            a = a.shiftRight(1);
            shift++;
        }
        return res;
    }

    static BigInteger clpow(BigInteger base, long exp) {
        BigInteger res = BigInteger.ONE;
        while (exp > 0) {
            if ((exp & 1L) == 1L) {
                res = clmul(res, base);
            }
            exp >>= 1L;
            if (exp > 0) {
                base = clmul(base, base);
            }
        }
        return res;
    }

    static long evalPolyMod(BigInteger polyBits, long xMod) {
        if (polyBits.signum() == 0)
            return 0;
        int deg = polyBits.bitLength() - 1;
        long ans = 0;
        for (int i = deg; i >= 0; i--) {
            ans = (ans * xMod) % MOD;
            if (polyBits.testBit(i)) {
                ans++;
                if (ans == MOD) {
                    ans = 0;
                }
            }
        }
        return ans;
    }

    public static String solve() {
        long threePow8 = 6561L;
        long twoPow52 = 1L << 52;

        BigInteger f = BigInteger.valueOf(11);
        BigInteger g = clpow(f, threePow8);

        long x = modPow(2, twoPow52);
        long ans = evalPolyMod(g, x);

        return Long.toString(ans);
    }

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