Problem 752: Powers of $1+\sqrt 7$

Project Euler Problem 752 Solution

EulerSolve provides an optimized solution for Project Euler Problem 752, Powers of $1+\sqrt 7$, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Write $$(1+\sqrt7)^k=a_k+b_k\sqrt7.$$ For each integer $n\ge 2$, let $g(n)$ be the least positive $k$ such that $$a_k\equiv 1 \pmod{n},\qquad b_k\equiv 0 \pmod{n};$$ if no such $k$ exists, then $g(n)=0$. The goal is to evaluate $$G(N)=\sum_{n=2}^{N} g(n)$$ for $N=10^6$. The C++, Python, and Java implementations do this by viewing $1+\sqrt7$ inside the finite ring $R_n=(\mathbb Z/n\mathbb Z)[\sqrt7]$ and computing its multiplicative order whenever that order exists. Mathematical Approach Let $$u=1+\sqrt7.$$ Then $g(n)$ is exactly the smallest positive exponent $k$ for which $u^k=1$ in $R_n$. Everything in the solution follows from understanding this order prime by prime, then rebuilding composite moduli with the Chinese remainder theorem. Step 1: Turn the Power Condition into Ring Arithmetic Represent an element $a+b\sqrt7$ by the pair $(a,b)$. Multiplication becomes $$(a,b)\times(c,d)=\bigl(ac+7bd,\ ad+bc\bigr)\pmod{n},$$ because $(a+b\sqrt7)(c+d\sqrt7)=(ac+7bd)+(ad+bc)\sqrt7$. So checking whether $(1+\sqrt7)^k\equiv 1$ is the same as checking whether repeated pair multiplication sends $(1,1)$ to $(1,0)$....

Detailed mathematical approach

Problem Summary

Write

$$(1+\sqrt7)^k=a_k+b_k\sqrt7.$$

For each integer $n\ge 2$, let $g(n)$ be the least positive $k$ such that

$$a_k\equiv 1 \pmod{n},\qquad b_k\equiv 0 \pmod{n};$$

if no such $k$ exists, then $g(n)=0$. The goal is to evaluate

$$G(N)=\sum_{n=2}^{N} g(n)$$

for $N=10^6$. The C++, Python, and Java implementations do this by viewing $1+\sqrt7$ inside the finite ring $R_n=(\mathbb Z/n\mathbb Z)[\sqrt7]$ and computing its multiplicative order whenever that order exists.

Mathematical Approach

Let

$$u=1+\sqrt7.$$

Then $g(n)$ is exactly the smallest positive exponent $k$ for which $u^k=1$ in $R_n$. Everything in the solution follows from understanding this order prime by prime, then rebuilding composite moduli with the Chinese remainder theorem.

Step 1: Turn the Power Condition into Ring Arithmetic

Represent an element $a+b\sqrt7$ by the pair $(a,b)$. Multiplication becomes

$$(a,b)\times(c,d)=\bigl(ac+7bd,\ ad+bc\bigr)\pmod{n},$$

because $(a+b\sqrt7)(c+d\sqrt7)=(ac+7bd)+(ad+bc)\sqrt7$.

So checking whether $(1+\sqrt7)^k\equiv 1$ is the same as checking whether repeated pair multiplication sends $(1,1)$ to $(1,0)$.

The norm of $a+b\sqrt7$ is

$$N(a+b\sqrt7)=a^2-7b^2.$$

For $u$ we get

$$N(u)=1^2-7\cdot 1^2=-6.$$

If $u^k=1$ modulo $n$, then norms give

$$(-6)^k\equiv 1 \pmod{n}.$$

This is impossible when $2\mid n$ or $3\mid n$, so every multiple of $2$ or $3$ has $g(n)=0$. Conversely, if $\gcd(n,6)=1$, then $-6$ is invertible modulo $n$, hence $u$ is a unit in a finite group and some positive order must exist.

Step 2: Prime Moduli and the Legendre Symbol

Let $p\ge 5$ be prime with $p\ne 7$. The key question is whether $7$ is a quadratic residue modulo $p$. By Euler's criterion,

$$\left(\frac{7}{p}\right)\equiv 7^{(p-1)/2}\pmod{p}\in\{1,-1\}.$$

If $\left(\frac{7}{p}\right)=1$, then $x^2-7$ splits into two distinct linear factors modulo $p$. Therefore

$$R_p\cong \mathbb F_p\times \mathbb F_p,$$

and every unit has order dividing $p-1$. In particular,

$$g(p)\mid (p-1).$$

If $\left(\frac{7}{p}\right)=-1$, then $x^2-7$ is irreducible modulo $p$, so $R_p$ is the field $\mathbb F_{p^2}$. Its multiplicative group has size $p^2-1$, hence

$$g(p)\mid (p^2-1)=(p-1)(p+1).$$

The implementation starts from the appropriate upper bound and strips away prime factors one by one. Whenever a candidate divisor $q$ satisfies $u^{M/q}=1$, the current order bound $M$ can be reduced by $q$.

Step 3: The Exceptional Prime $7$ and Order Lifting

The prime $7$ behaves differently because

$$x^2-7\equiv x^2 \pmod{7}.$$

If we write $\varepsilon=\sqrt7$ modulo $7$, then $\varepsilon^2=0$, so the binomial theorem collapses to

$$(1+\varepsilon)^k=1+k\varepsilon \pmod{7}.$$

This equals $1$ exactly when $k\equiv 0\pmod{7}$, therefore

$$g(7)=7.$$

For higher prime powers, let $t=g(p^e)$. Reducing modulo $p^e$ shows that $g(p^{e+1})$ must be a multiple of $t$, while the extra factor can only come from one more power of $p$. So the only possibilities are

$$g(p^{e+1})\in\{t,\ pt\}.$$

The implementation checks this directly: if $u^t\equiv 1\pmod{p^{e+1}}$, the order stays $t$; otherwise it is multiplied by $p$.

Step 4: Composite Moduli from the Chinese Remainder Theorem

Suppose

$$n=\prod_{i=1}^{r} p_i^{e_i},\qquad \gcd(n,6)=1.$$

The Chinese remainder theorem gives the ring decomposition

$$R_n\cong \prod_{i=1}^{r} R_{p_i^{e_i}}.$$

The image of $u$ in this product has one component in each prime-power ring. An exponent kills the whole product exactly when it kills every component, so

$$g(n)=\operatorname{lcm}\bigl(g(p_1^{e_1}),g(p_2^{e_2}),\dots,g(p_r^{e_r})\bigr).$$

This is why the algorithm spends most of its effort on prime powers: once those orders are known, every composite $n$ is reconstructed by factorization plus one $\operatorname{lcm}$ chain.

Step 5: Worked Example with $n=35$

Take

$$n=35=5\cdot 7.$$

For $p=5$, the Legendre symbol is $\left(\frac{7}{5}\right)=\left(\frac{2}{5}\right)=-1$, so $g(5)\mid 24$. Using pair arithmetic modulo $5$,

$$u^2=(3,2),\qquad u^3=(2,0).$$

Therefore

$$u^{12}=(2,0)^4=(1,0),\qquad u^6=(2,0)^2=(4,0)\ne (1,0),$$

$$g(5)=12.$$

For $p=7$, Step 3 already showed that

$$g(7)=7.$$

Now combine the two prime components:

$$g(35)=\operatorname{lcm}(12,7)=84.$$

This mirrors the full solution on a small modulus: compute prime orders, lift if necessary, then join them with the Chinese remainder theorem.

Step 6: From Local Orders to the Global Sum

Because every multiple of $2$ or $3$ contributes $0$, the total can be written as

$$G(N)=\sum_{\substack{2\le n\le N\\ \gcd(n,6)=1}} g(n).$$

So the problem is no longer "search powers separately for every $n$". Instead it becomes a preprocessing problem: build all prime-power orders up to $N$, factor each admissible $n$, take the least common multiple of its prime-power contributions, and accumulate the answer.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they represent $a+b\sqrt7$ as a pair and use binary exponentiation to test whether a candidate exponent sends $u$ back to $1$ modulo a chosen modulus. This makes order tests cheap enough to repeat many times.

Next they build a smallest-prime-factor sieve up to $N$. For each prime $p\ge 5$, they determine whether $7$ is a quadratic residue using Euler's criterion, choose the correct upper bound for $g(p)$, and reduce that bound by testing prime divisors. The exceptional prime $7$ is handled separately at the base level, and then every prime order is lifted through $p^2,p^3,\dots$ with the rule from Step 3.

Finally, for each $n\le N$ with $\gcd(n,6)=1$, the sieve provides the factorization of $n$. The implementation looks up the cached order for each prime power, takes their $\operatorname{lcm}$, and adds the result to the running total.

Complexity Analysis

Let the search limit be $N$. Building the smallest-prime-factor sieve costs $O(N\log\log N)$ time and $O(N)$ memory. Factoring each admissible $n$ then takes $O(\log N)$ time on average. The prime and prime-power preprocessing requires repeated binary exponentiation for a modest number of candidate divisors, so the overall running time is close to $O(N\log N)$ in practice, with $O(N)$ memory for the sieve and cached prime-power orders.

Footnotes and References

Problem page: https://projecteuler.net/problem=752
Legendre symbol: Wikipedia - Legendre symbol
Euler's criterion: Wikipedia - Euler's criterion
Finite field: Wikipedia - Finite field
Multiplicative order: Wikipedia - Multiplicative order
Chinese remainder theorem: Wikipedia - Chinese remainder theorem

Problem 752 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Pair {
    u64 a;
    u64 b;
};

Pair mul(const Pair x, const Pair y, const u64 mod) {
    return {
        (x.a * y.a + 7ULL * x.b * y.b) % mod,
        (x.a * y.b + x.b * y.a) % mod,
    };
}

Pair pow_u(const u64 exp, const u64 mod) {
    Pair result{1ULL, 0ULL};
    Pair base{1ULL, 1ULL};
    u64 e = exp;
    while (e > 0ULL) {
        if ((e & 1ULL) != 0ULL) {
            result = mul(result, base, mod);
        }
        base = mul(base, base, mod);
        e >>= 1ULL;
    }
    return result;
}

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1));
    for (int i = 0; i <= n; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }
    for (int i = 2; i * i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return spf;
}

std::vector<int> distinct_prime_factors(int n, const std::vector<int>& spf) {
    std::vector<int> factors;
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        factors.push_back(p);
        while (n % p == 0) {
            n /= p;
        }
    }
    return factors;
}

u64 order_mod_prime(const int p, const std::vector<int>& spf) {
    if (p == 7) {
        return 7ULL;
    }

    const u64 ls = [&]() -> u64 {
        u64 r = 1ULL;
        u64 b = 7ULL % static_cast<u64>(p);
        u64 e = static_cast<u64>((p - 1) / 2);
        const u64 mod = static_cast<u64>(p);
        while (e > 0ULL) {
            if ((e & 1ULL) != 0ULL) {
                r = (r * b) % mod;
            }
            b = (b * b) % mod;
            e >>= 1ULL;
        }
        return r;
    }();

    u64 exponent = 0ULL;
    std::vector<int> factors;
    if (ls == 1ULL) {
        exponent = static_cast<u64>(p - 1);
        factors = distinct_prime_factors(p - 1, spf);
    } else {
        exponent = static_cast<u64>(p - 1) * static_cast<u64>(p + 1);
        factors = distinct_prime_factors(p - 1, spf);
        std::vector<int> extra = distinct_prime_factors(p + 1, spf);
        factors.insert(factors.end(), extra.begin(), extra.end());
        std::sort(factors.begin(), factors.end());
        factors.erase(std::unique(factors.begin(), factors.end()), factors.end());
    }

    u64 ord = exponent;
    for (const int q : factors) {
        const u64 qq = static_cast<u64>(q);
        while (ord % qq == 0ULL) {
            const Pair test = pow_u(ord / qq, static_cast<u64>(p));
            if (test.a == 1ULL && test.b == 0ULL) {
                ord /= qq;
            } else {
                break;
            }
        }
    }
    return ord;
}

u64 G(const int n) {
    const int limit = n + 1;
    const std::vector<int> spf = build_spf(limit);

    std::vector<u64> ord_prime_power(static_cast<std::size_t>(n + 1), 0ULL);

    for (int p = 5; p <= n; ++p) {
        if (spf[static_cast<std::size_t>(p)] != p) {
            continue;
        }

        u64 ord = order_mod_prime(p, spf);
        int pk = p;
        ord_prime_power[static_cast<std::size_t>(pk)] = ord;

        while (pk <= n / p) {
            pk *= p;
            const Pair test = pow_u(ord, static_cast<u64>(pk));
            if (!(test.a == 1ULL && test.b == 0ULL)) {
                ord *= static_cast<u64>(p);
            }
            ord_prime_power[static_cast<std::size_t>(pk)] = ord;
        }
    }

    u64 total = 0ULL;
    for (int x = 2; x <= n; ++x) {
        if ((x % 2) == 0 || (x % 3) == 0) {
            continue;
        }

        int t = x;
        u64 gx = 1ULL;
        while (t > 1) {
            const int p = spf[static_cast<std::size_t>(t)];
            int pk = 1;
            while (t % p == 0) {
                t /= p;
                pk *= p;
            }
            gx = std::lcm(gx, ord_prime_power[static_cast<std::size_t>(pk)]);
        }
        total += gx;
    }

    return total;
}

}  // namespace

int main() {
    assert(G(100) == 28'891ULL);
    assert(G(1'000) == 13'131'583ULL);

    std::cout << G(1'000'000) << '\n';
    return 0;
}

Python

import math

def mul(x, y, mod):
    return (
        (x[0] * y[0] + 7 * x[1] * y[1]) % mod,
        (x[0] * y[1] + x[1] * y[0]) % mod
    )

def pow_u(exp, mod):
    result = (1, 0)
    base = (1, 1)
    e = exp
    while e > 0:
        if e & 1:
            result = mul(result, base, mod)
        base = mul(base, base, mod)
        e >>= 1
    return result

def build_spf(n):
    spf = list(range(n + 1))
    for i in range(2, int(math.isqrt(n)) + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf

def distinct_prime_factors(n, spf):
    factors = []
    while n > 1:
        p = spf[n]
        factors.append(p)
        while n % p == 0:
            n //= p
    return factors

def order_mod_prime(p, spf):
    if p == 7:
        return 7
        
    def ls():
        r = 1
        b = 7 % p
        e = (p - 1) // 2
        while e > 0:
            if e & 1:
                r = (r * b) % p
            b = (b * b) % p
            e >>= 1
        return r
        
    ls_val = ls()
    
    if ls_val == 1:
        exponent = p - 1
        factors = distinct_prime_factors(p - 1, spf)
    else:
        exponent = (p - 1) * (p + 1)
        factors = distinct_prime_factors(p - 1, spf)
        extra = distinct_prime_factors(p + 1, spf)
        factors.extend(extra)
        factors = sorted(list(set(factors)))
        
    ord_val = exponent
    for q in factors:
        while ord_val % q == 0:
            test = pow_u(ord_val // q, p)
            if test[0] == 1 and test[1] == 0:
                ord_val //= q
            else:
                break
                
    return ord_val

def G(n):
    limit = n + 1
    spf = build_spf(limit)
    ord_prime_power = [0] * (n + 1)
    
    for p in range(5, n + 1):
        if spf[p] != p:
            continue
            
        ord_val = order_mod_prime(p, spf)
        pk = p
        ord_prime_power[pk] = ord_val
        
        while pk <= n // p:
            pk *= p
            test = pow_u(ord_val, pk)
            if not (test[0] == 1 and test[1] == 0):
                ord_val *= p
            ord_prime_power[pk] = ord_val
            
    total = 0
    for x in range(2, n + 1):
        if x % 2 == 0 or x % 3 == 0:
            continue
            
        t = x
        gx = 1
        while t > 1:
            p = spf[t]
            pk = 1
            while t % p == 0:
                t //= p
                pk *= p
            gx = (gx * ord_prime_power[pk]) // math.gcd(gx, ord_prime_power[pk])
            
        total += gx
        
    return total

def solve():
    return str(G(1000000))

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler752 {
    static class Pair {
        long a, b;

        Pair(long a, long b) {
            this.a = a;
            this.b = b;
        }
    }

    static Pair mul(Pair x, Pair y, long mod) {
        return new Pair(
                (x.a * y.a + 7L * x.b * y.b) % mod,
                (x.a * y.b + x.b * y.a) % mod);
    }

    static Pair powU(long exp, long mod) {
        Pair result = new Pair(1L, 0L);
        Pair base = new Pair(1L, 1L);
        long e = exp;
        while (e > 0L) {
            if ((e & 1L) != 0L) {
                result = mul(result, base, mod);
            }
            base = mul(base, base, mod);
            e >>= 1L;
        }
        return result;
    }

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; ++i) {
            spf[i] = i;
        }
        for (int i = 2; i * i <= n; ++i) {
            if (spf[i] != i)
                continue;
            for (int j = i * i; j <= n; j += i) {
                if (spf[j] == j) {
                    spf[j] = i;
                }
            }
        }
        return spf;
    }

    static List<Integer> distinctPrimeFactors(int n, int[] spf) {
        List<Integer> factors = new ArrayList<>();
        while (n > 1) {
            int p = spf[n];
            factors.add(p);
            while (n % p == 0) {
                n /= p;
            }
        }
        return factors;
    }

    static long orderModPrime(int p, int[] spf) {
        if (p == 7) {
            return 7L;
        }

        long r = 1L;
        long b = 7L % p;
        long e = (p - 1) / 2;
        long pMod = p;
        while (e > 0L) {
            if ((e & 1L) != 0L) {
                r = (r * b) % pMod;
            }
            b = (b * b) % pMod;
            e >>= 1L;
        }
        long ls = r;

        long exponent;
        List<Integer> factors;
        if (ls == 1L) {
            exponent = p - 1;
            factors = distinctPrimeFactors(p - 1, spf);
        } else {
            exponent = (long) (p - 1) * (p + 1);
            factors = distinctPrimeFactors(p - 1, spf);
            factors.addAll(distinctPrimeFactors(p + 1, spf));
            Collections.sort(factors);
            List<Integer> uniqueFactors = new ArrayList<>();
            for (int f : factors) {
                if (uniqueFactors.isEmpty() || uniqueFactors.get(uniqueFactors.size() - 1) != f) {
                    uniqueFactors.add(f);
                }
            }
            factors = uniqueFactors;
        }

        long ord = exponent;
        for (int q : factors) {
            long qq = (long) q;
            while (ord % qq == 0L) {
                Pair test = powU(ord / qq, p);
                if (test.a == 1L && test.b == 0L) {
                    ord /= qq;
                } else {
                    break;
                }
            }
        }
        return ord;
    }

    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 long G(int n) {
        int limit = n + 1;
        int[] spf = buildSpf(limit);
        long[] ordPrimePower = new long[n + 1];

        for (int p = 5; p <= n; ++p) {
            if (spf[p] != p)
                continue;

            long ord = orderModPrime(p, spf);
            int pk = p;
            ordPrimePower[pk] = ord;

            while (pk <= n / p) {
                pk *= p;
                Pair test = powU(ord, pk);
                if (!(test.a == 1L && test.b == 0L)) {
                    ord *= p;
                }
                ordPrimePower[pk] = ord;
            }
        }

        long total = 0L;
        for (int x = 2; x <= n; ++x) {
            if (x % 2 == 0 || x % 3 == 0)
                continue;

            int t = x;
            long gx = 1L;
            while (t > 1) {
                int p = spf[t];
                int pk = 1;
                while (t % p == 0) {
                    t /= p;
                    pk *= p;
                }
                gx = lcm(gx, ordPrimePower[pk]);
            }
            total += gx;
        }

        return total;
    }

    public static String solve() {
        return Long.toString(G(1000000));
    }

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