Problem 411: Uphill Paths

View on Project Euler

Project Euler Problem 411 Solution

EulerSolve provides an optimized solution for Project Euler Problem 411, Uphill Paths, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a given \(n\), define the stations $$P_i=\left(2^i \bmod n,\ 3^i \bmod n\right),\qquad 0 \le i \le 2n.$$ An uphill path from \((0,0)\) to \((n,n)\) uses only right and up moves, so the stations visited by such a path must appear in an order whose coordinates never decrease. Let \(S(n)\) be the maximum number of stations that can be visited by one uphill path. The program evaluates $$\sum_{k=1}^{30} S(k^5).$$ Mathematical Approach The implementation does not search over grid paths directly. Instead, it reduces the problem to a finite set of modular points and then solves an order problem on those points. 1) Reduce the Station Sequence to a Finite Prefix Write $$n = 2^{\alpha} 3^{\beta} m,\qquad \gcd(m,6)=1.$$ Consider the \(x\)-coordinate \(x_i = 2^i \bmod n\). For every \(i \ge \alpha\), the factor \(2^{\alpha}\) already divides \(2^i\), so modulo \(2^{\alpha}\) the value is fixed at \(0\). On the coprime part \(n / 2^{\alpha}\), the base \(2\) is invertible, hence the remaining behavior is purely periodic with period $$T_2 = \operatorname{ord}_{n / 2^{\alpha}}(2),$$ with the harmless convention \(T_2=1\) when \(n / 2^{\alpha}=1\). Thus the \(x\)-coordinate is eventually periodic after \(\alpha\) steps. Exactly the same argument applies to the \(y\)-coordinate \(y_i = 3^i \bmod n\)....

Detailed mathematical approach

Problem Summary

For a given \(n\), define the stations

$$P_i=\left(2^i \bmod n,\ 3^i \bmod n\right),\qquad 0 \le i \le 2n.$$

An uphill path from \((0,0)\) to \((n,n)\) uses only right and up moves, so the stations visited by such a path must appear in an order whose coordinates never decrease. Let \(S(n)\) be the maximum number of stations that can be visited by one uphill path. The program evaluates

$$\sum_{k=1}^{30} S(k^5).$$

Mathematical Approach

The implementation does not search over grid paths directly. Instead, it reduces the problem to a finite set of modular points and then solves an order problem on those points.

1) Reduce the Station Sequence to a Finite Prefix

Write

$$n = 2^{\alpha} 3^{\beta} m,\qquad \gcd(m,6)=1.$$

Consider the \(x\)-coordinate \(x_i = 2^i \bmod n\). For every \(i \ge \alpha\), the factor \(2^{\alpha}\) already divides \(2^i\), so modulo \(2^{\alpha}\) the value is fixed at \(0\). On the coprime part \(n / 2^{\alpha}\), the base \(2\) is invertible, hence the remaining behavior is purely periodic with period

$$T_2 = \operatorname{ord}_{n / 2^{\alpha}}(2),$$

with the harmless convention \(T_2=1\) when \(n / 2^{\alpha}=1\). Thus the \(x\)-coordinate is eventually periodic after \(\alpha\) steps.

Exactly the same argument applies to the \(y\)-coordinate \(y_i = 3^i \bmod n\). After \(\beta\) steps it is periodic with period

$$T_3 = \operatorname{ord}_{n / 3^{\beta}}(3),$$

again taking \(T_3=1\) when the modulus equals \(1\).

Therefore the pair sequence \((x_i,y_i)\) is periodic from

$$p = \max(\alpha,\beta)$$

onward, with common period

$$T = \operatorname{lcm}(T_2,T_3).$$

Every station with index \(i \ge p+T\) repeats one already seen in the block \(p, p+1, \dots, p+T-1\). Since the original problem only uses indices \(0 \le i \le 2n\), it is enough to generate

$$M = \min(2n+1,\ p+T)$$

points and then discard duplicates.

2) Convert Uphill Paths into a Longest Nondecreasing Subsequence

After deduplication, sort the distinct stations lexicographically by \(x\) and then by \(y\):

$$Q_1,\dots,Q_r,\qquad Q_j=(x_j,y_j),\qquad x_1 \le x_2 \le \cdots \le x_r.$$

If an uphill path visits stations \(Q_{j_1}, Q_{j_2}, \dots, Q_{j_t}\), then necessarily

$$x_{j_1} \le x_{j_2} \le \cdots \le x_{j_t},\qquad y_{j_1} \le y_{j_2} \le \cdots \le y_{j_t}.$$

Because the list is already sorted by \(x\), the first condition is automatic when we take an increasing index subsequence. So the entire problem becomes

$$S(n)=\operatorname{LNDS}(y_1,y_2,\dots,y_r),$$

where \(\operatorname{LNDS}\) denotes the length of the longest nondecreasing subsequence.

The word nondecreasing matters. Equal \(x\)-coordinates are allowed because a path can move vertically between two stations, and equal \(y\)-coordinates are allowed because it can move horizontally. That is why the implementation uses the patience-sorting variant based on the first tail strictly greater than the new \(y\)-value.

3) Worked Example: \(n=22\)

For \(n=22=2\cdot 11\), we have \(\alpha=1\) and \(\beta=0\). Hence

$$T_2=\operatorname{ord}_{11}(2)=10,\qquad T_3=\operatorname{ord}_{22}(3)=5,$$

so

$$p=1,\qquad T=\operatorname{lcm}(10,5)=10,\qquad M=\min(45,11)=11.$$

Thus all relevant stations already occur among the first \(11\) indices. After sorting and removing repeats, the \(y\)-sequence is

$$1,3,9,15,5,1,1,5,15,9,3.$$

Its longest nondecreasing subsequence has length \(5\), so

$$S(22)=5.$$

This matches the checkpoint used by the implementation. The same program also verifies \(S(123)=14\) and \(S(10000)=48\).

How the Code Works

The C++, Python, and Java implementations build a smallest-prime-factor table once up to \(30^5\). That table makes Euler's totient function and multiplicative orders fast to evaluate for every modulus \(n=k^5\).

For each \(k\in\{1,\dots,30\}\), the implementation forms \(n=k^5\), computes the preperiod length \(p\) and the common period \(T\), generates the first \(M\) stations, sorts and deduplicates them, and finally applies patience sorting with binary search to the \(y\)-coordinates. The Python entry point reuses the same compiled computation, while the Java version mirrors the same mathematics directly.

Complexity Analysis

For one modulus \(n\), let \(M=\min(2n+1,p+T)\), and let \(r\le M\) be the number of distinct stations after deduplication. Generating the modular points costs \(O(M)\). Sorting the distinct points costs \(O(r\log r)\), and the longest-nondecreasing-subsequence step also costs \(O(r\log r)\). Therefore one evaluation of \(S(n)\) runs in

$$O(r\log r)$$

time and uses \(O(r)\) memory.

Across all \(k \le 30\), the one-time prime-factor sieve up to \(30^5\) is near-linear and uses \(O(30^5)\) memory. After that precomputation, the point ordering work dominates.

References

  1. Problem page: https://projecteuler.net/problem=411
  2. Multiplicative order: Wikipedia — Multiplicative order
  3. Euler's totient function: Wikipedia — Euler's totient function
  4. Chinese remainder theorem: Wikipedia — Chinese remainder theorem
  5. Longest increasing subsequence and patience sorting: Wikipedia — Longest increasing subsequence

Problem 411 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

struct Options {
    int k_max = 30;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--k-max=", options.k_max)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.k_max >= 1;
}

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(n / 10);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (const int p : primes) {
            const i64 x = 1LL * i * p;
            if (x > n) {
                break;
            }
            spf[static_cast<std::size_t>(x)] = p;
            if (p == spf[static_cast<std::size_t>(i)]) {
                break;
            }
        }
    }
    return spf;
}

void radix_sort_u64(std::vector<u64>& values) {
    if (values.size() <= 1U) {
        return;
    }

    std::vector<u64> tmp(values.size());
    constexpr int RADIX_BITS = 16;
    constexpr int RADIX_SIZE = 1 << RADIX_BITS;
    constexpr u64 RADIX_MASK = static_cast<u64>(RADIX_SIZE - 1);
    std::array<std::uint32_t, RADIX_SIZE> count{};

    for (int pass = 0; pass < 4; ++pass) {
        const int shift = pass * RADIX_BITS;
        count.fill(0);

        for (const u64 v : values) {
            ++count[static_cast<std::size_t>((v >> shift) & RADIX_MASK)];
        }

        std::uint32_t acc = 0;
        for (int i = 0; i < RADIX_SIZE; ++i) {
            const std::uint32_t c = count[static_cast<std::size_t>(i)];
            count[static_cast<std::size_t>(i)] = acc;
            acc += c;
        }

        for (const u64 v : values) {
            tmp[count[static_cast<std::size_t>((v >> shift) & RADIX_MASK)]++] = v;
        }
        values.swap(tmp);
    }
}

u64 mod_pow_u64(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0) {
        if (e & 1ULL) {
            result = static_cast<u64>((__uint128_t)result * cur % mod);
        }
        cur = static_cast<u64>((__uint128_t)cur * cur % mod);
        e >>= 1ULL;
    }
    return result;
}

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

int euler_phi(int x, const std::vector<int>& spf) {
    if (x == 1) {
        return 1;
    }
    int phi = x;
    int n = x;
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        phi = phi / p * (p - 1);
        while (n % p == 0) {
            n /= p;
        }
    }
    return phi;
}

int multiplicative_order(int base, int mod, const std::vector<int>& spf) {
    if (mod == 1) {
        return 1;
    }
    int ord = euler_phi(mod, spf);
    std::vector<int> fac = unique_prime_factors(ord, spf);
    for (const int p : fac) {
        while (ord % p == 0) {
            const int cand = ord / p;
            if (mod_pow_u64(static_cast<u64>(base), static_cast<u64>(cand), static_cast<u64>(mod)) == 1ULL) {
                ord = cand;
            } else {
                break;
            }
        }
    }
    return ord;
}

int S_value(int n, const std::vector<int>& spf) {
    if (n == 1) {
        return 1;
    }

    int pre2 = 0;
    int tmp = n;
    while ((tmp & 1) == 0) {
        ++pre2;
        tmp >>= 1;
    }
    const int mod2 = tmp;
    const int per2 = multiplicative_order(2, mod2, spf);

    int pre3 = 0;
    tmp = n;
    while (tmp % 3 == 0) {
        ++pre3;
        tmp /= 3;
    }
    const int mod3 = tmp;
    const int per3 = multiplicative_order(3, mod3, spf);

    const int pre = std::max(pre2, pre3);
    const i64 period = std::lcm(static_cast<i64>(per2), static_cast<i64>(per3));
    const i64 needed = std::min<i64>(2LL * n + 1, pre + period);

    std::vector<u64> points;
    points.reserve(static_cast<std::size_t>(needed));

    int x = 1 % n;
    int y = 1 % n;
    for (i64 i = 0; i < needed; ++i) {
        const u64 key = (static_cast<u64>(static_cast<std::uint32_t>(x)) << 32) |
                        static_cast<u64>(static_cast<std::uint32_t>(y));
        points.push_back(key);
        x = static_cast<int>((2LL * x) % n);
        y = static_cast<int>((3LL * y) % n);
    }

    radix_sort_u64(points);
    points.erase(std::unique(points.begin(), points.end()), points.end());

    std::vector<int> tails;
    tails.reserve(256);
    for (const u64 key : points) {
        const int yy = static_cast<int>(key & 0xffffffffULL);
        auto it = std::upper_bound(tails.begin(), tails.end(), yy);
        if (it == tails.end()) {
            tails.push_back(yy);
        } else {
            *it = yy;
        }
    }

    return static_cast<int>(tails.size());
}

i64 solve_sum(int k_max, const std::vector<int>& spf) {
    i64 ans = 0;
    for (int k = 1; k <= k_max; ++k) {
        int n = 1;
        for (int t = 0; t < 5; ++t) {
            n *= k;
        }
        ans += S_value(n, spf);
    }
    return ans;
}

bool run_checkpoints(const std::vector<int>& spf) {
    if (S_value(22, spf) != 5) {
        std::cerr << "Checkpoint failed: S(22)\n";
        return false;
    }
    if (S_value(123, spf) != 14) {
        std::cerr << "Checkpoint failed: S(123)\n";
        return false;
    }
    if (S_value(10000, spf) != 48) {
        std::cerr << "Checkpoint failed: S(10000)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const int n_max = 24300000;  // 30^5
    const std::vector<int> spf = build_spf(n_max);

    if (options.run_checkpoints && !run_checkpoints(spf)) {
        return 2;
    }

    std::cout << solve_sum(options.k_max, spf) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;

public class Euler411 {
    static int[] buildSPF(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long x = (long) i * p;
                if (x > n)
                    break;
                spf[(int) x] = p;
                if (p == spf[i])
                    break;
            }
        }
        return spf;
    }

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

    static int eulerPhi(int x, int[] spf) {
        int phi = x, n = x;
        while (n > 1) {
            int p = spf[n];
            phi = phi / p * (p - 1);
            while (n % p == 0)
                n /= p;
        }
        return phi;
    }

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

    static int multOrder(int base, int mod, int[] spf) {
        if (mod == 1)
            return 1;
        int ord = eulerPhi(mod, spf);
        for (int p : uniquePF(ord, spf)) {
            while (ord % p == 0) {
                int cand = ord / p;
                if (modPow(base, cand, mod) == 1)
                    ord = cand;
                else
                    break;
            }
        }
        return ord;
    }

    static int sValue(int n, int[] spf) {
        if (n == 1)
            return 1;
        int pre2 = 0, tmp = n;
        while ((tmp & 1) == 0) {
            pre2++;
            tmp >>= 1;
        }
        int mod2 = tmp, per2 = multOrder(2, mod2, spf);
        int pre3 = 0;
        tmp = n;
        while (tmp % 3 == 0) {
            pre3++;
            tmp /= 3;
        }
        int mod3 = tmp, per3 = multOrder(3, mod3, spf);
        int pre = Math.max(pre2, pre3);
        long period = lcm(per2, per3);
        long needed = Math.min(2L * n + 1, pre + period);

        long[] points = new long[(int) needed];
        int x = 1 % n, y = 1 % n;
        for (int i = 0; i < needed; i++) {
            points[i] = ((long) x << 32) | (y & 0xFFFFFFFFL);
            x = (int) ((2L * x) % n);
            y = (int) ((3L * y) % n);
        }
        Arrays.sort(points);
        // Remove duplicates
        int uLen = 0;
        for (int i = 0; i < points.length; i++) {
            if (i == 0 || points[i] != points[i - 1])
                points[uLen++] = points[i];
        }
        // LIS on y-values
        int[] tails = new int[uLen];
        int lisLen = 0;
        for (int i = 0; i < uLen; i++) {
            int yy = (int) (points[i] & 0xFFFFFFFFL);
            int pos = upperBound(tails, lisLen, yy);
            if (pos == lisLen)
                tails[lisLen++] = yy;
            else
                tails[pos] = yy;
        }
        return lisLen;
    }

    static int upperBound(int[] arr, int len, int val) {
        int lo = 0, hi = len;
        while (lo < hi) {
            int mid = (lo + hi) >>> 1;
            if (arr[mid] <= val)
                lo = mid + 1;
            else
                hi = mid;
        }
        return lo;
    }

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

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

    public static String solve() {
        int kMax = 30;
        int nMax = 1;
        for (int t = 0; t < 5; t++)
            nMax *= kMax;
        int[] spf = buildSPF(nMax);
        long ans = 0;
        for (int k = 1; k <= kMax; k++) {
            int n = 1;
            for (int t = 0; t < 5; t++)
                n *= k;
            ans += sValue(n, spf);
        }
        return String.valueOf(ans);
    }

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