Problem 160: Factorial Trailing Digits

Project Euler Problem 160 Solution

EulerSolve provides an optimized solution for Project Euler Problem 160, Factorial Trailing Digits, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $F(n)$ denote the last five non-zero digits of $n!$. The problem asks for $F(10^{12})$, so any method that tries to form the factorial directly is hopelessly too large. The successful approach is to separate the trailing zeros from the rest of the factorial, compute the reduced factorial modulo $5^5=3125$, observe what happens modulo $2^5=32$, and then combine the two congruences into the final residue modulo $10^5$. Mathematical Approach Write $v_p(m)$ for the exponent of the prime $p$ in $m$. The whole computation is driven by the fact that trailing zeros come from matching one factor 2 with one factor 5. Separate trailing zeros from the surviving digits The number of trailing zeros in $n!$ is the number of available factors 5, because factors 2 are more abundant. Therefore $$e=v_5(n!)=\sum_{k\ge 1}\left\lfloor \frac{n}{5^k}\right\rfloor,$$ and the desired quantity is $$F(n)\equiv \frac{n!}{10^e}\pmod{100000}.$$ This identity says exactly what must be preserved: remove the $e$ matched pairs $(2,5)$, keep everything else, and then read the result modulo $10^5$. Remove all factors of 5 recursively It is convenient to strip out the factors of 5 first and delay the matching factors of 2 until later. Define $$U(n)=\prod_{m=1}^{n}\frac{m}{5^{v_5(m)}}.$$ Then $U(n)=n!/5^e$....

Detailed mathematical approach

Problem Summary

Write $N!$ in decimal, remove the zeros at its end and keep the last five of the remaining digits: these five digits are $f(N)$. The statement gives $9!=362880$, so $f(9)=36288$; $10!=3628800$, so $f(10)=36288$; and $20!=2432902008176640000$, so $f(20)=17664$. We must find $f(10^{12})$. ^[1] In arithmetic terms, if $N!$ ends in exactly $z$ zeros, then $f(N)$ is the remainder of $N!/10^{z}$ after division by $10^5=100000$, written with five digits. The number $(10^{12})!$ has about $1.16\cdot10^{13}$ digits; we need only five of them, so we work with remainders instead of the whole number.

The whole solution in one sentence: we never compute $N!$ itself, only the remainders of $N!/10^{z}$ after division by $32$ and by $3125$, whose product is $100000$; the first remainder is $0$ for every $N\ge8$, the second follows from a short recursion over the base-$5$ digits of $N$ and a table of products modulo $3125$, and the Chinese remainder theorem joins the two into the five digits.

Mathematical Approach

We write $v_p(m)$ for the exponent of the prime $p$ in $m$, the number of times $p$ divides $m$: for example $v_2(40)=3$ and $v_5(40)=1$, because $40=2^3\cdot5$. We write $a\equiv b\pmod m$ when $m$ divides $a-b$, that is, when $a$ and $b$ leave the same remainder after division by $m$.

1. Trailing zeros are pairs of a 2 and a 5

Since $10=2\cdot5$, the number $N!$ is divisible by $10^k$ exactly when it is divisible by $2^k$ and by $5^k$, so it ends in $z=\min(v_2(N!),v_5(N!))$ zeros. Among $1,2,\dots,N$ there are $\lfloor N/p\rfloor$ multiples of $p$, $\lfloor N/p^2\rfloor$ multiples of $p^2$, and so on, and a number with exactly $e$ factors $p$ is counted in the first $e$ of these counts. This gives Legendre's formula ^[3]

$$v_p(N!)=\lfloor N/p\rfloor+\lfloor N/p^2\rfloor+\lfloor N/p^3\rfloor+\cdots.$$

Every term for $p=2$ is at least the corresponding term for $p=5$, so $v_2(N!)\ge v_5(N!)$ and $z=v_5(N!)$: the factors $5$ are the scarce ones. ^[2] For $N=20$ we get $z=4+0=4$ and $v_2(20!)=10+5+2+1=18$; indeed $20!$ ends in four zeros, and $20!/10^4=243290200817664$ ends in $17664$. For $N=10^{12}$ the formula gives $z=249999999997$. So the task is to compute

$$f(N)=\frac{N!}{10^{z}}\bmod 10^5,\qquad z=v_5(N!).$$

2. Why the five digits cannot simply be carried along

A first idea is to multiply by $1,2,3,\dots$ in turn, remove the zeros after each step and keep only the last five digits. This fails as early as $N=15$. From $14!=87178291200$ we keep $82912=f(14)$, and $82912\cdot15=1243680$ suggests $f(15)=24368$. But $15!=1307674368000$, so $f(15)=74368$. The factor $15$ created a new zero, and after its removal the fifth digit from the end depends on the sixth digit from the end of $14!/10^2=871782912$, which had already been thrown away. Removing zeros is a division by $10$, and a division needs the digits that the reduction discarded.

The cure is to strip every $m\le N$ of its factors $2$ and $5$, only counting them, and let $d(N)=v_2(N!)-v_5(N!)$ be the number of factors $2$ that have no partner $5$. Then $N!/10^{z}$ is exactly $2^{d(N)}$ times the product of the stripped numbers, a formula with multiplications only, so every step may be reduced modulo $10^5$. This is correct, but for $N=10^{12}$ it needs $10^{12}$ steps. The following sections replace it by a small table, built once, and a few dozen steps.

3. Splitting 100000 into 32 and 3125

We have $10^5=2^5\cdot5^5=32\cdot3125$, and $32$ and $3125$ have no common factor. By the Chinese remainder theorem, a number $x$ with $0\le x\lt100000$ is determined by its remainders after division by $32$ and by $3125$. ^[4] For example, $17664=552\cdot32$ and $17664=5\cdot3125+2039$, so $f(20)$ is the only number below $100000$ with the remainders $0$ and $2039$.

4. The remainder modulo 32 is zero

After the $z$ pairs are removed, $N!/10^{z}$ still contains $d(N)$ factors $2$, so it is a multiple of $32$ exactly when $d(N)\ge5$. For $N=1,2,\dots,12$ the values of $d(N)$ are $0,1,1,3,2,3,3,6,6,6,6,8$; they drop at $N=5$, because the factor $5$ brings a $5$ but no $2$. Still, $d(N)\ge5$ for every $N\ge8$. Indeed $v_5(N!)\lt N/5+N/25+\dots=N/4$, while the first two terms of Legendre's formula give

$$v_2(N!)\ge\lfloor N/2\rfloor+\lfloor N/4\rfloor\ge\frac{N-1}2+\frac{N-3}4,$$

so that

$$d(N)\gt\frac{N-1}2+\frac{N-3}4-\frac N4=\frac N2-\frac54.$$

For $N\ge12$ the right-hand side is at least $4.75$, so the integer $d(N)$ is at least $5$; the list above covers $N=8,\dots,11$. For $N\le7$ the remainder is not $0$; for example $f(7)=504=15\cdot32+24$. The programs therefore compute $f(N)$ directly for $N\lt10$. For $N=10^{12}$ we have $d=749999999990$.

5. Removing the factors 5 one base-5 digit at a time

Let $U(N)=N!/5^{z}$ be the factorial with all its factors $5$ removed, and let $A(N)$ be the product of the numbers $m\le N$ that are not divisible by $5$. The multiples of $5$ up to $N$ are $5\cdot1,5\cdot2,\dots,5\cdot\lfloor N/5\rfloor$, so

$$N!=A(N)\cdot5^{\lfloor N/5\rfloor}\cdot\lfloor N/5\rfloor!.$$

The factor $A(N)$ contains no $5$, so removing all factors $5$ on both sides gives

$$U(N)=A(N)\,U(\lfloor N/5\rfloor),\qquad U(0)=1.$$

Dividing $N$ by $5$ and rounding down deletes the last base-$5$ digit of $N$, so $U(N)$ is a product with one factor for each base-$5$ digit of $N$:

$$U(N)=A(N)\,A(\lfloor N/5\rfloor)\,A(\lfloor N/25\rfloor)\cdots.$$

For $N=10$, the numbers $1,\dots,10$ without their factors $5$ are $1,2,3,4,1,6,7,8,9,2$, whose product is $U(10)=145152$${}=A(10)\,A(2)=72576\cdot2$; dividing by $2^{z}=2^2$, the partners of the two factors $5$, gives $36288=f(10)$. In base $5$, the number $10^{12}$ is written as $112341$ followed by twelve zeros, so for the target the product has $18$ factors.

6. A block of 3125 numbers contributes $-1$

Modulo $3125$ only the remainders of the factors of $A(N)$ matter, and these factors are exactly the numbers coprime to $3125=5^5$; their remainders are called units. Every unit $u$ has a partner $u'$ with $uu'\equiv1\pmod{3125}$, its inverse: if $ua\equiv ub$, then $3125$ divides $u(a-b)$ and, as $u$ has no factor $5$, also $a-b$; so $u\cdot0,u\cdot1,\dots,u\cdot3124$ leave $3125$ different remainders, and one of them is $1$. The smaller modulus $25$ shows what happens. The $20$ numbers from $1$ to $24$ that are not multiples of $5$ contain the nine pairs $2\cdot13$, $3\cdot17$, $4\cdot19$, $6\cdot21$, $7\cdot18$, $8\cdot22$, $9\cdot14$, $11\cdot16$ and $12\cdot23$, whose products $26$, $51$, $76$, $126$, $126$, $176$, $126$, $176$ and $276$ are all $\equiv1\pmod{25}$, while $1$ and $24$ are their own partners, since $24\cdot24=576\equiv1$. Hence the product of all twenty is $\equiv24\equiv-1\pmod{25}$.

The same happens modulo $3125$. A number that is its own partner satisfies $u^2\equiv1$, so $3125$ divides $(u-1)(u+1)$. The two factors differ by $2$, so at most one of them is divisible by $5$, and that one must contain all of $5^5$: $u\equiv1$ or $u\equiv-1$. Every other unit differs from its partner, so the units pair off into products $\equiv1$, and the product of all $2500$ numbers from $1$ to $3125$ that are coprime to $3125$ is $\equiv1\cdot(-1)=-1\pmod{3125}$. This is Gauss's generalisation of Wilson's theorem, which gives $-1$ for every power of an odd prime. ^[5] Any $3125$ consecutive integers contain one number from each remainder class, so their $2500$ members not divisible by $5$ also multiply to $-1$ modulo $3125$.

7. One table for every level

Write $N=3125q+r$ with $0\le r\lt3125$. The numbers up to $N$ consist of $q$ complete blocks of $3125$, each contributing $-1$ by Section 6, and the numbers $3125q+t$ with $1\le t\le r$, which leave the same remainders as $t$. Hence

$$A(N)\equiv(-1)^{q}\,P(r)\pmod{3125},$$

where $P(r)$ is the product of the numbers $t\le r$ not divisible by $5$, reduced modulo $3125$, and $P(0)=1$. One table of $P(0),\dots,P(3125)$ therefore gives $A$ at every level of the recursion; its last entry, $P(3125)=3124$, is the $-1$ of Section 6. For example, $A(3130)$ consists of one complete block and the numbers $3126,\dots,3129$, so $A(3130)\equiv-P(5)=-24$${}\equiv3101\pmod{3125}$.

8. Dividing by the matching factors 2

Section 5 removed the factors $5$ of the $z$ trailing zeros; their partners $2^{z}$ are still inside $U(N)$, and $N!/10^{z}=U(N)/2^{z}$. Since $2$ and $3125$ are coprime, dividing by $2^{z}$ modulo $3125$ means multiplying by the inverse of $2^{z}$. Euler's theorem says that $a^{\varphi(m)}\equiv1\pmod m$ for every $a$ coprime to $m$, where $\varphi(m)$ counts the numbers from $1$ to $m$ that are coprime to $m$. ^[6] Here $\varphi(3125)=3125-625=2500$, because exactly the $625$ multiples of $5$ are not coprime to $3125$. ^[7] Writing $z=2500c+t$ with $0\le t\lt2500$, we get $2^{z}\equiv2^{t}$ and $2^{2500-t}\cdot2^{t}=2^{2500}\equiv1$, so

$$2^{-z}\equiv2^{\,2500-t}\pmod{3125},\qquad t=z\bmod2500.$$

For $z=4$ this is $2^{2496}\equiv586$, and indeed $16\cdot586=9376=3\cdot3125+1$. For $N=10^{12}$ we have $z\equiv2497\pmod{2500}$, so the multiplier is simply $2^{3}=8$. Such powers are computed by repeated squaring. A result starts at $1$ and a base at $2$; while the exponent is positive, the result is multiplied by the base if the exponent is odd, then the base is squared and the exponent halved, rounding down. For the exponent $13=8+4+1$ the base runs through $2,4,16,256$, and the result collects $2\cdot16\cdot256=2^{13}$. An exponent of at most $2500$ has at most $12$ binary digits, so the loop runs at most $12$ times. ^[8]

9. Joining the two remainders

Let $R$, with $0\le R\lt3125$, be the remainder of $N!/10^{z}$ modulo $3125$. The numbers below $100000$ with this remainder are $x=R+3125k$ for $k=0,1,\dots,31$, and since $3125=97\cdot32+21$, they satisfy $x\equiv R+21k\pmod{32}$. As $21\cdot29=609=19\cdot32+1$, the number $29$ is the inverse of $21$ modulo $32$, so $x\equiv0\pmod{32}$ holds for exactly one of the $32$ candidates, the one with

$$k\equiv-29R\pmod{32}.$$

By Sections 3 and 4 this $x$ is $f(N)$ whenever $N\ge8$. The inverse replaces the sieving search along $R,R+3125,R+6250,\dots$ by one step. ^[4]

10. The whole computation for $N=20$

Legendre's formula gives $z=4$, and Section 5 gives $U(20)=A(20)\,A(4)$. Both arguments are below $3125$, so $q=0$ and the table entries $P(20)=2401$ and $P(4)=24$ are used directly: $U(20)\equiv2401\cdot24=57624$${}\equiv1374\pmod{3125}$. Section 8 gives $R\equiv1374\cdot586\equiv2039$ modulo $3125$. Then $2039=63\cdot32+23$, so $k\equiv-29\cdot23\equiv5\pmod{32}$, and $x=2039+5\cdot3125=17664$, which is $f(20)$ from the statement.

How the Code Works

The C++ program ^[10] starts in main with init_P_table, which fills the global array P_table with the table $P(0),\dots,P(3125)$ of Section 7: entry $i$ is entry $i-1$ times $i$ modulo $3125$ when $5\nmid i$, and a copy of entry $i-1$ otherwise. get_P_mod_3125(n) splits $n=3125q+r$ and returns P_table[r] times $1$ or $3124$ modulo $3125$, depending on whether $q$ is even or odd; this is $A(n)$ by Section 7. factorial_no_5_mod_3125(n) is the recursion of Section 5: it returns $1$ for $n=0$ and otherwise multiplies get_P_mod_3125(n) by its own value at n / 5; for $10^{12}$ it makes $19$ calls, one per base-$5$ digit and a last one with $n=0$. count_factors_5 evaluates Legendre's formula by adding n / 5 and dividing $n$ by $5$ until $n$ is $0$, and power is the square-and-multiply loop of Section 8.

f(n) combines these: two_exp is $t=z\bmod2500$, inv_2 is $2^{2500-t}\bmod3125$, and res3125 is $R$. For $n\ge10$, solve_crt computes target as $(32-R\bmod32)\bmod32$, which is $-R\bmod32$, then $k$ as $29$ times target modulo $32$, and returns $R+3125k$, as in Section 9. For $n\lt10$, f instead multiplies $1,\dots,n$, removes the trailing zeros by dividing by $10$ and returns the remainder modulo $100000$; this covers the cases $n\le7$ of Section 4, in which the remainder modulo $32$ is not $0$. The C++ version computes the modular part even for $n\lt10$ and then discards it. main computes $f(9)$, $f(10)$ and $f(20)$, prints them next to the values of the statement and stops with exit code $1$ if one of them differs; otherwise it computes $f(10^{12})$ and prints it with setw(5) and setfill('0'), so a leading zero would be kept. The file also defines a helper modInverse and the constants MOD and M1, which the program never uses.

The Python program builds the same table inside solve and defines the other functions there, with the same names; its f tests $n\lt10$ first, and the result is printed with the format 05d. The Java program has the static methods power, initP, getPmod, factNo5, count5, crt and f, does its arithmetic in long, with int only for loop counters and array indices, and prints with %05d. Neither of them checks the statement values. Every product the programs form is below $3125^2\lt10^7$, so the 64-bit integers of C++ (int64, an alias of long long) and Java cannot overflow, and Python's integers have no fixed size; only $N$, its first few quotients by powers of $5$ and $z$ need more than 32 bits. The recursion is $19$ calls deep, far below Python's default limit of $1000$.

Complexity and Verification

Building the table takes $3125$ loop steps, $2500$ of them multiplications, once. For one value of $N$, the recursion and the loop of count_factors_5 each take one step per base-$5$ digit of $N$, that is $\lfloor\log_5N\rfloor+1=18$ steps for $10^{12}$; power needs at most $12$ rounds, and joining the remainders takes a few operations. Counting arithmetic operations and stored integers, time and memory are therefore $O(5^5+\log N)$; the memory is the table of $3126$ entries and a recursion stack of $19$ calls. All three programs finish at once.

The C++ program checks $f(9)$, $f(10)$ and $f(20)$ from the statement before it computes the answer, and all three programs print the same five digits. For this article we also ran checks that are not part of the published solutions:

For every $8\le N\le10^5$, a separate implementation of Sections 5 to 9 agrees with $N!/10^{z}$ kept as an exact integer.
For every $N\le10^9$, the function f of the C++ file agrees with the safe product of Section 2, and so do the Java f for $N\le10^7$ and the Python f for $N\le10^6$. Without the branch for small $n$, the joining step of Section 9 gives the right value for every $8\le N\le10^9$ and a wrong one for every $N\le7$, as Section 4 predicts.
A second method works modulo $100000$ directly, without the Chinese remainder theorem, the inverse of $2$ or the sign of Section 6. Every $m\le N$ is written uniquely as $2^i5^jt$ with $t$ coprime to $10$, so $N!/10^{z}$ is $2^{d(N)}$ times the product of $C(\lfloor N/(2^i5^j)\rfloor)$ over the pairs $(i,j)$ with $2^i5^j\le N$, where $C(n)$ is the product of the numbers up to $n$ that are coprime to $10$; for $N=10^{12}$ there are $371$ such pairs. This method agrees with the C++ f on $720$ other values of $N$ up to about $10^{12}$, among them $10^k$ for $k\le11$ and the $100$ integers within $50$ of $10^{12}$, and with all three programs at $N=10^{12}$.
The last digits of $f(0),\dots,f(29)$ agree with the final non-zero digits of $n!$ listed in OEIS A008904. ^[9]

Footnotes and References

The references below supply the background facts used above. The recursion over the base-$5$ digits, the table with the sign of complete blocks, the bound on $d(N)$ and the operation counts are derived explicitly in this article.

Project Euler 160 — Factorial Trailing Digits. The official statement defines f(N) as the last five digits before the trailing zeros in N!, gives 9! = 362880 and 10! = 3628800 with f(9) = f(10) = 36288 and 20! = 2432902008176640000 with f(20) = 17664, and asks for f(1 000 000 000 000).
Wikipedia — Trailing zero. The number of trailing zeros of n! is the multiplicity of the prime factor 5 in n!, ⌊n/5⌋ + ⌊n/5²⌋ + ⋯, because there are at least as many factors 2; for example, 32! ends in seven zeros.
Wikipedia — Legendre's formula. The exponent of a prime p in n! is the sum of ⌊n/pⁱ⌋ over i ≥ 1, a sum with only finitely many non-zero terms.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N and prescribed remainders there is exactly one integer x with 0 ≤ x < N that leaves these remainders. The search by sieving runs through a₁, a₁ + n₁, a₁ + 2n₁, … until the remainder modulo n₂ is right; Section 9 does this in one step.
Wikipedia — Wilson's theorem. By Gauss's generalisation, the product of the positive integers less than m and coprime to m is one less than a multiple of m when m is 4, a power of an odd prime or twice a power of an odd prime, and one more than a multiple of m otherwise. The elementary proof for a prime p arranges the factors other than ±1 in pairs whose products are congruent to 1, the argument of Section 6.
Wikipedia — Euler's theorem. If a and n are coprime positive integers, then a^φ(n) ≡ 1 (mod n), so exponents of a can be reduced modulo φ(n).
Wikipedia — Euler's totient function. For a prime p, φ(pᵏ) = pᵏ − pᵏ⁻¹, because the numbers from 1 to pᵏ that are not coprime to pᵏ are the pᵏ⁻¹ multiples of p; for 3125 = 5⁵ this gives 2500.
Wikipedia — Modular exponentiation. The right-to-left binary method multiplies the result by the base when the lowest bit of the exponent is 1, then squares the base and shifts the exponent; its running time is O(log e) for the exponent e.
OEIS A008904 — final nonzero digit of n!. Its terms begin 1, 1, 2, 6, 4, 2, 2, 4, 2, 8, 8, 8, 6, 8, 2, and the entry gives formulas for them in terms of the base-5 digits of n.
C++ — Euler160.cpp. The linked immutable C++ revision contains init_P_table, get_P_mod_3125, factorial_no_5_mod_3125, count_factors_5, solve_crt, f and the checks of f(9), f(10) and f(20) in main. It is the implementation record for this article.

Problem 160 source code

C++

#include <algorithm>
#include <cmath>
#include <future>
#include <iomanip>
#include <iostream>
#include <vector>

using namespace std;

using int64 = long long;

const int64 MOD = 100000;
const int64 M1 = 32;
const int64 M2 = 3125;
const int64 PHI_M2 = 2500;

int64 P_table[M2 + 1];

int64 power(int64 base, int64 exp, int64 mod) {
  int64 res = 1;
  base %= mod;
  while (exp > 0) {
    if (exp % 2 == 1)
      res = (res * base) % mod;
    base = (base * base) % mod;
    exp /= 2;
  }
  return res;
}

int64 modInverse(int64 n, int64 mod) {
  return power(n, mod > 2000 ? PHI_M2 - 1 : mod / 2, mod);
}

void init_P_table() {
  P_table[0] = 1;
  for (int64 i = 1; i <= M2; ++i) {
    if (i % 5 != 0) {
      P_table[i] = (P_table[i - 1] * i) % M2;
    } else {
      P_table[i] = P_table[i - 1];
    }
  }
}

int64 get_P_mod_3125(int64 n) {
  int64 q = n / M2;
  int64 r = n % M2;
  int64 res = (q % 2 == 0) ? 1 : (M2 - 1);
  res = (res * P_table[r]) % M2;
  return res;
}

int64 factorial_no_5_mod_3125(int64 n) {
  if (n == 0)
    return 1;
  int64 part1 = get_P_mod_3125(n);
  int64 part2 = factorial_no_5_mod_3125(n / 5);
  return (part1 * part2) % M2;
}

int64 count_factors_5(int64 n) {
  int64 count = 0;
  while (n > 0) {
    count += n / 5;
    n /= 5;
  }
  return count;
}

int64 solve_crt(int64 res3125) {
  int64 inv21 = 29;
  int64 target = (32 - (res3125 % 32)) % 32;
  int64 k = (target * inv21) % 32;
  return res3125 + k * M2;
}

int64 f(int64 n) {
  int64 fact_no_5 = factorial_no_5_mod_3125(n);
  int64 num_5s = count_factors_5(n);

  int64 two_exp = num_5s % PHI_M2;
  int64 inv_2 = power(2, PHI_M2 - two_exp, M2);

  int64 res3125 = (fact_no_5 * inv_2) % M2;

  if (n >= 10) {
    return solve_crt(res3125);
  } else {
    int64 res = 1;
    for (int i = 1; i <= n; ++i)
      res *= i;
    while (res % 10 == 0)
      res /= 10;
    return res % 100000;
  }
}

int main() {
  init_P_table();

  auto handle9 = std::async(std::launch::async, f, 9);
  auto handle10 = std::async(std::launch::async, f, 10);
  auto handle20 = std::async(std::launch::async, f, 20);

  int64 r9 = handle9.get();
  int64 r10 = handle10.get();
  int64 r20 = handle20.get();

  std::cout << "Validation:" << std::endl;
  std::cout << "f(9)  = " << r9 << " (Expected: 36288)" << std::endl;
  std::cout << "f(10) = " << r10 << " (Expected: 36288)" << std::endl;
  std::cout << "f(20) = " << r20 << " (Expected: 17664)" << std::endl;

  bool ok = (r9 == 36288 && r10 == 36288 && r20 == 17664);
  if (!ok) {
    std::cerr << "Validation Failed!" << std::endl;
    return 1;
  }

  std::cout << "Validation Passed. Computing f(10^12)..." << std::endl;

  int64 N = 1000000000000LL;
  int64 result = f(N);

  std::cout << "f(10^12) = " << std::setfill('0') << std::setw(5) << result
            << std::endl;

  return 0;
}

Python

# Problem 160: Factorial trailing digits
# Find the last 5 non-zero digits of 10^12!

def solve():
    MOD = 100000
    M2 = 3125
    PHI_M2 = 2500
    
    def power(base, exp, mod):
        result = 1
        base %= mod
        while exp > 0:
            if exp % 2 == 1:
                result = (result * base) % mod
            base = (base * base) % mod
            exp //= 2
        return result
    
    # Build P_table: product of integers 1..M2 coprime to 5, mod 3125
    P_table = [0] * (M2 + 1)
    P_table[0] = 1
    for i in range(1, M2 + 1):
        if i % 5 != 0:
            P_table[i] = (P_table[i-1] * i) % M2
        else:
            P_table[i] = P_table[i-1]
    
    def get_P_mod_3125(n):
        q = n // M2
        r = n % M2
        res = 1 if q % 2 == 0 else (M2 - 1)
        return (res * P_table[r]) % M2
    
    def factorial_no_5_mod_3125(n):
        if n == 0: return 1
        return (get_P_mod_3125(n) * factorial_no_5_mod_3125(n // 5)) % M2
    
    def count_factors_5(n):
        count = 0
        while n > 0:
            count += n // 5
            n //= 5
        return count
    
    def solve_crt(res3125):
        inv21 = 29  # 21^(-1) mod 32
        target = (32 - (res3125 % 32)) % 32
        k = (target * inv21) % 32
        return res3125 + k * M2
    
    def f(n):
        if n < 10:
            res = 1
            for i in range(1, n + 1):
                res *= i
            while res % 10 == 0:
                res //= 10
            return res % MOD
        
        fact_no_5 = factorial_no_5_mod_3125(n)
        num_5s = count_factors_5(n)
        two_exp = num_5s % PHI_M2
        inv_2 = power(2, PHI_M2 - two_exp, M2)
        res3125 = (fact_no_5 * inv_2) % M2
        return solve_crt(res3125)
    
    N = 10**12
    result = f(N)
    print(f"{result:05d}")

solve()

Java

public class Euler160 {
    static final long M2 = 3125, PHI = 2500;
    static long[] P = new long[(int) M2 + 1];

    static long power(long b, long e, long m) {
        long r = 1;
        b %= m;
        while (e > 0) {
            if (e % 2 == 1)
                r = r * b % m;
            b = b * b % m;
            e /= 2;
        }
        return r;
    }

    static void initP() {
        P[0] = 1;
        for (int i = 1; i <= (int) M2; i++)
            P[i] = i % 5 != 0 ? (P[i - 1] * i) % M2 : P[i - 1];
    }

    static long getPmod(long n) {
        long q = n / M2, r = n % M2;
        return ((q % 2 == 0 ? 1 : M2 - 1) * P[(int) r]) % M2;
    }

    static long factNo5(long n) {
        return n == 0 ? 1 : (getPmod(n) * factNo5(n / 5)) % M2;
    }

    static long count5(long n) {
        long c = 0;
        while (n > 0) {
            c += n / 5;
            n /= 5;
        }
        return c;
    }

    static long crt(long r) {
        long t = (32 - (r % 32)) % 32;
        long k = (t * 29) % 32;
        return r + k * M2;
    }

    static long f(long n) {
        if (n < 10) {
            long r = 1;
            for (int i = 1; i <= n; i++)
                r *= i;
            while (r % 10 == 0)
                r /= 10;
            return r % 100000;
        }
        long fn5 = factNo5(n), n5 = count5(n), inv2 = power(2, PHI - n5 % PHI, M2);
        return crt((fn5 * inv2) % M2);
    }

    public static void main(String[] args) {
        initP();
        System.out.printf("%05d%n", f(1000000000000L));
    }
}