Problem 383: Divisibility Comparison Between Factorials

View on Project Euler

Project Euler Problem 383 Solution

EulerSolve provides an optimized solution for Project Euler Problem 383, Divisibility Comparison Between Factorials, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Define $$T_5(N)=\#\left\{1\le i\le N: v_5((2i-1)!) \lt 2v_5(i!)\right\}.$$ We must evaluate \(T_5(10^{18})\). A direct scan over all \(i\) is impossible, so the implemented solution turns the factorial inequality into a statement about carries in base \(5\), then counts the valid numbers with a memoized digit DP. Mathematical Approach Step 1: Rewrite the factorial valuation Because \((2i)! = 2i\,(2i-1)!\), we have $$v_5((2i-1)!) = v_5((2i)!)-v_5(2i).$$ Subtracting \(2v_5(i!)\) gives $$v_5((2i-1)!)-2v_5(i!)=v_5((2i)!)-2v_5(i!)-v_5(2i).$$ The middle part is exactly the \(5\)-adic valuation of the central binomial coefficient: $$v_5((2i)!)-2v_5(i!)=v_5\left(\binom{2i}{i}\right).$$ Since \(v_5(2)=0\), we also have \(v_5(2i)=v_5(i)\). Therefore the original condition is equivalent to $$v_5\left(\binom{2i}{i}\right) \lt v_5(i).$$ Legendre's formula explains why factorial valuations are natural here, but the code uses the binomial rewrite because it exposes a carry-counting interpretation immediately. Step 2: Apply Kummer's theorem in base \(5\) Kummer's theorem states that \(v_5\left(\binom{2i}{i}\right)\) equals the number of carries when adding \(i+i\) in base \(5\). Write $$i=5^a m,\qquad a=v_5(i),\qquad 5\nmid m.$$ Multiplying by \(5^a\) merely appends \(a\) trailing zeros in base \(5\), so it does not create or destroy carries; it only shifts their positions....

Detailed mathematical approach

Problem Summary

Define

$$T_5(N)=\#\left\{1\le i\le N: v_5((2i-1)!) \lt 2v_5(i!)\right\}.$$

We must evaluate \(T_5(10^{18})\). A direct scan over all \(i\) is impossible, so the implemented solution turns the factorial inequality into a statement about carries in base \(5\), then counts the valid numbers with a memoized digit DP.

Mathematical Approach

Step 1: Rewrite the factorial valuation

Because \((2i)! = 2i\,(2i-1)!\), we have

$$v_5((2i-1)!) = v_5((2i)!)-v_5(2i).$$

Subtracting \(2v_5(i!)\) gives

$$v_5((2i-1)!)-2v_5(i!)=v_5((2i)!)-2v_5(i!)-v_5(2i).$$

The middle part is exactly the \(5\)-adic valuation of the central binomial coefficient:

$$v_5((2i)!)-2v_5(i!)=v_5\left(\binom{2i}{i}\right).$$

Since \(v_5(2)=0\), we also have \(v_5(2i)=v_5(i)\). Therefore the original condition is equivalent to

$$v_5\left(\binom{2i}{i}\right) \lt v_5(i).$$

Legendre's formula explains why factorial valuations are natural here, but the code uses the binomial rewrite because it exposes a carry-counting interpretation immediately.

Step 2: Apply Kummer's theorem in base \(5\)

Kummer's theorem states that \(v_5\left(\binom{2i}{i}\right)\) equals the number of carries when adding \(i+i\) in base \(5\). Write

$$i=5^a m,\qquad a=v_5(i),\qquad 5\nmid m.$$

Multiplying by \(5^a\) merely appends \(a\) trailing zeros in base \(5\), so it does not create or destroy carries; it only shifts their positions. Hence

$$\operatorname{carries}_5(i+i)=\operatorname{carries}_5(m+m).$$

The inequality becomes

$$\operatorname{carries}_5(m+m) \lt a,$$

and because the carry count is an integer, this is the same as

$$\operatorname{carries}_5(m+m)\le a-1.$$

Step 3: Sum over the \(5\)-adic order of \(i\)

Every positive integer \(i\) has a unique decomposition \(i=5^a m\) with \(5\nmid m\). For a fixed \(a\ge 1\), the admissible values of \(m\) satisfy

$$1\le m\le \left\lfloor\frac{N}{5^a}\right\rfloor,\qquad 5\nmid m,$$

and they contribute exactly when \(\operatorname{carries}_5(m+m)\le a-1\). Define \(G(x,k)\) to be the number of positive integers \(m\le x\), not divisible by \(5\), whose doubling in base \(5\) produces at most \(k\) carries. Then

$$\boxed{T_5(N)=\sum_{a\ge 1} G\left(\left\lfloor\frac{N}{5^a}\right\rfloor,a-1\right).}$$

The sum stops as soon as \(5^a \gt N\). For \(N=10^{18}\), only \(25\) values of \(a\) need to be examined.

Step 4: Derive the least-significant-digit DP

The programs process digits from least significant to most significant. Let \(F(x,k,c)\) denote the number of integers \(0\le n\le x\) such that, when doubling \(n\) in base \(5\), the still-unprocessed higher digits create at most \(k\) additional carries, assuming an incoming carry \(c\in\{0,1\}\) from the lower digits.

If the current digit is \(d\in\{0,1,2,3,4\}\), the outgoing carry is

$$c'=\left\lfloor\frac{2d+c}{5}\right\rfloor.$$

Writing \(n=d+5q\) gives \(q\le \left\lfloor\frac{x-d}{5}\right\rfloor\), so the recurrence is

$$F(x,k,c)=\sum_{\substack{0\le d\le 4\\ d\le x}} F\left(\left\lfloor\frac{x-d}{5}\right\rfloor,k-c',c'\right).$$

The base cases are exactly the ones in the source files:

$$F(x,k,c)=0\quad \text{if }x\lt 0\text{ or }k\lt 0,\qquad F(0,k,c)=1\quad \text{for }k\ge 0.$$

To enforce \(5\nmid m\), the least significant digit must be \(1,2,3,\) or \(4\). Therefore

$$G(x,k)=\sum_{\substack{1\le d\le 4\\ d\le x}} F\left(\left\lfloor\frac{x-d}{5}\right\rfloor,k-\left\lfloor\frac{2d}{5}\right\rfloor,\left\lfloor\frac{2d}{5}\right\rfloor\right).$$

This is precisely what the helper count_non_multiples_of_five computes.

Step 5: Small examples

If \(i=20=5\cdot 4\), then \(a=1\) and \(m=4\). In base \(5\),

$$\left(4\right)_5+\left(4\right)_5=\left(13\right)_5,$$

so there is one carry. Since \(1 \nleq 0\), this \(i\) does not contribute. By contrast, \(i=35=5\cdot 7\) has \(m=\left(12\right)_5\), and

$$\left(12\right)_5+\left(12\right)_5=\left(24\right)_5,$$

which creates no carries, so it does contribute for \(a=1\). This is exactly the local criterion used by the DP.

How the Code Works

The C++, Python, and Java files all implement the same recurrence. carry_out returns \(c'=\lfloor(2d+c)/5\rfloor\). count_with_initial_carry is the memoized function \(F(x,k,c)\), keyed by \((x,k,\text{carry\_in})\). count_non_multiples_of_five handles the first digit separately so the counted number is positive and not divisible by \(5\). The outer loop iterates over \(a=1,2,\dots\) while \(5^a\le N\), accumulating \(G(\lfloor N/5^a\rfloor,a-1)\). The C++ implementation also validates the method with checkpoints \(T_5(10^3)=68\), \(T_5(10^9)=2408210\), and a brute-force comparison at \(N=20000\).

Complexity Analysis

Let \(A=\lfloor\log_5 N\rfloor\). The outer sum has \(A\) terms, and each recursive step removes one base-\(5\) digit. Every memo state has only \(5\) outgoing transitions and is determined by a triple \((x,k,c)\) with \(c\in\{0,1\}\). Thus the cost is proportional to the number of reachable memo states rather than to \(N\) itself, and the memory usage is the same order as the memo table. For the actual target \(N=10^{18}\), the reachable state space is tiny, so the method is effectively instantaneous compared with brute force.

References

  1. Problem page: https://projecteuler.net/problem=383
  2. Kummer's theorem: Wikipedia — Kummer's theorem
  3. Legendre's formula: Wikipedia — Legendre's formula
  4. \(p\)-adic valuation: Wikipedia — \(p\)-adic valuation

Problem 383 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <unordered_map>

namespace {

using i64 = long long;

struct Options {
    i64 n = 1000000000000000000LL;
    bool run_checkpoints = true;
};

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

    i64 parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<i64>(ch - '0');
    }
    value = parsed;
    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, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1;
}

int carry_out(const int digit, const int carry_in) {
    return (2 * digit + carry_in) / 5;
}

struct Key {
    i64 x;
    int k;
    int carry_in;

    bool operator==(const Key& other) const {
        return x == other.x && k == other.k && carry_in == other.carry_in;
    }
};

struct KeyHash {
    std::size_t operator()(const Key& key) const {
        std::size_t h = std::hash<i64>{}(key.x);
        h ^= std::hash<int>{}(key.k) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        h ^= std::hash<int>{}(key.carry_in) + 0x9e3779b97f4a7c15ULL + (h << 6) + (h >> 2);
        return h;
    }
};

class Solver {
   public:
    i64 solve(const i64 n) {
        memo_.clear();
        i64 answer = 0;
        i64 power_of_five = 5;

        for (int a = 1; power_of_five <= n; ++a) {
            answer += count_non_multiples_of_five(n / power_of_five, a - 1);
            if (power_of_five > n / 5) {
                break;
            }
            power_of_five *= 5;
        }

        return answer;
    }

   private:
    std::unordered_map<Key, i64, KeyHash> memo_;

    i64 count_with_initial_carry(const i64 x, const int k, const int carry_in) {
        if (x < 0 || k < 0) {
            return 0;
        }
        if (x == 0) {
            // Only n=0 is present, and it contributes zero further carries.
            return 1;
        }

        const Key key{x, k, carry_in};
        const auto it = memo_.find(key);
        if (it != memo_.end()) {
            return it->second;
        }

        i64 ways = 0;
        for (int digit = 0; digit <= 4; ++digit) {
            if (x < digit) {
                break;
            }
            const i64 q = (x - digit) / 5;
            const int next_carry = carry_out(digit, carry_in);
            ways += count_with_initial_carry(q, k - next_carry, next_carry);
        }

        memo_[key] = ways;
        return ways;
    }

    i64 count_non_multiples_of_five(const i64 x, const int k) {
        if (x <= 0 || k < 0) {
            return 0;
        }
        i64 ways = 0;
        for (int digit = 1; digit <= 4; ++digit) {
            if (x < digit) {
                break;
            }
            const i64 q = (x - digit) / 5;
            const int next_carry = carry_out(digit, 0);
            ways += count_with_initial_carry(q, k - next_carry, next_carry);
        }
        return ways;
    }
};

int v5_factorial(i64 n) {
    int total = 0;
    while (n > 0) {
        n /= 5;
        total += static_cast<int>(n);
    }
    return total;
}

i64 brute_force_count(const i64 n) {
    i64 count = 0;
    for (i64 i = 1; i <= n; ++i) {
        const int lhs = v5_factorial(2 * i - 1);
        const int rhs = 2 * v5_factorial(i);
        if (lhs < rhs) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    Solver solver;

    if (solver.solve(1000) != 68LL) {
        std::cerr << "Checkpoint failed: T5(10^3)\n";
        return false;
    }

    if (solver.solve(1000000000LL) != 2408210LL) {
        std::cerr << "Checkpoint failed: T5(10^9)\n";
        return false;
    }

    const i64 brute_n = 20000;
    const i64 expected = brute_force_count(brute_n);
    const i64 got = solver.solve(brute_n);
    if (got != expected) {
        std::cerr << "Checkpoint failed against brute force at n=" << brute_n
                  << " (expected " << expected << ", got " << got << ")\n";
        return false;
    }

    return true;
}

}  // namespace

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

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

    Solver solver;
    std::cout << solver.solve(options.n) << '\n';
    return 0;
}

Python

def solve():
    memo = {}

    def carry_out(digit, carry_in):
        return (2 * digit + carry_in) // 5

    def count_with_initial_carry(x, k, carry_in):
        if x < 0 or k < 0:
            return 0
        if x == 0:
            return 1

        key = (x, k, carry_in)
        if key in memo:
            return memo[key]

        ways = 0
        for digit in range(5):
            if x < digit:
                break
            q = (x - digit) // 5
            next_carry = carry_out(digit, carry_in)
            ways += count_with_initial_carry(q, k - next_carry, next_carry)

        memo[key] = ways
        return ways

    def count_non_multiples_of_five(x, k):
        if x <= 0 or k < 0:
            return 0
        ways = 0
        for digit in range(1, 5):
            if x < digit:
                break
            q = (x - digit) // 5
            next_carry = carry_out(digit, 0)
            ways += count_with_initial_carry(q, k - next_carry, next_carry)
        return ways

    n = 1000000000000000000
    answer = 0
    power_of_five = 5

    a = 1
    while power_of_five <= n:
        answer += count_non_multiples_of_five(n // power_of_five, a - 1)
        if power_of_five > n // 5:
            break
        power_of_five *= 5
        a += 1

    return str(answer)

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

Java

import java.util.HashMap;
import java.util.Map;

public class Euler383 {

    private static class Key {
        long x;
        int k;
        int carryIn;

        Key(long x, int k, int carryIn) {
            this.x = x;
            this.k = k;
            this.carryIn = carryIn;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (!(o instanceof Key))
                return false;
            Key key = (Key) o;
            return x == key.x && k == key.k && carryIn == key.carryIn;
        }

        @Override
        public int hashCode() {
            int result = (int) (x ^ (x >>> 32));
            result = 31 * result + k;
            result = 31 * result + carryIn;
            return result;
        }
    }

    private static final Map<Key, Long> memo = new HashMap<>();

    private static int carryOut(int digit, int carryIn) {
        return (2 * digit + carryIn) / 5;
    }

    private static long countWithInitialCarry(long x, int k, int carryIn) {
        if (x < 0 || k < 0)
            return 0;
        if (x == 0)
            return 1;

        Key key = new Key(x, k, carryIn);
        Long cached = memo.get(key);
        if (cached != null)
            return cached;

        long ways = 0;
        for (int digit = 0; digit <= 4; ++digit) {
            if (x < digit)
                break;
            long q = (x - digit) / 5;
            int nextCarry = carryOut(digit, carryIn);
            ways += countWithInitialCarry(q, k - nextCarry, nextCarry);
        }

        memo.put(key, ways);
        return ways;
    }

    private static long countNonMultiplesOfFive(long x, int k) {
        if (x <= 0 || k < 0)
            return 0;
        long ways = 0;
        for (int digit = 1; digit <= 4; ++digit) {
            if (x < digit)
                break;
            long q = (x - digit) / 5;
            int nextCarry = carryOut(digit, 0);
            ways += countWithInitialCarry(q, k - nextCarry, nextCarry);
        }
        return ways;
    }

    public static String solve() {
        memo.clear();
        long n = 1000000000000000000L;
        long answer = 0;
        long powerOfFive = 5;

        for (int a = 1; powerOfFive <= n; ++a) {
            answer += countNonMultiplesOfFive(n / powerOfFive, a - 1);
            if (powerOfFive > n / 5)
                break;
            powerOfFive *= 5;
        }

        return String.valueOf(answer);
    }

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