From e943355e2b22e02fb49bcd202dab5a266f791479 Mon Sep 17 00:00:00 2001
From: Nayuki Minase
Date: Tue, 28 Jul 2015 05:08:08 +0000
Subject: [PATCH] P120, P123, P124, P128, P129, P132, P133, P250, P301, P429:
Added Python solutions.
---
python/eulertest.py | 10 +++
python/p120.py | 36 +++++++++
python/p123.py | 32 ++++++++
python/p124.py | 24 ++++++
python/p128.py | 185 ++++++++++++++++++++++++++++++++++++++++++++
python/p129.py | 58 ++++++++++++++
python/p132.py | 29 +++++++
python/p133.py | 64 +++++++++++++++
python/p250.py | 30 +++++++
python/p301.py | 53 +++++++++++++
python/p429.py | 54 +++++++++++++
11 files changed, 575 insertions(+)
create mode 100644 python/p120.py
create mode 100644 python/p123.py
create mode 100644 python/p124.py
create mode 100644 python/p128.py
create mode 100644 python/p129.py
create mode 100644 python/p132.py
create mode 100644 python/p133.py
create mode 100644 python/p250.py
create mode 100644 python/p301.py
create mode 100644 python/p429.py
diff --git a/python/eulertest.py b/python/eulertest.py
index 1f150f9..e4b0eda 100644
--- a/python/eulertest.py
+++ b/python/eulertest.py
@@ -100,7 +100,17 @@ ANSWERS = {
115: "168",
116: "20492570929",
117: "100808458960497",
+ 120: "333082500",
+ 123: "21035",
+ 124: "21417",
+ 128: "14516824220",
+ 129: "1000023",
+ 132: "843296",
+ 133: "453647705",
218: "0",
+ 250: "1425480602091519",
+ 301: "2178309",
+ 429: "98792821",
}
diff --git a/python/p120.py b/python/p120.py
new file mode 100644
index 0000000..be3c227
--- /dev/null
+++ b/python/p120.py
@@ -0,0 +1,36 @@
+#
+# Solution to Project Euler problem 120
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+# For a given a, what is the n that maximizes the remainder, and what is the value of this remainder?
+#
+# Let's simplify one term, mod a^2:
+# (a+1)^n = 1^n + (n choose 1) a 1^(n-1) + (n choose 2) a^2 1^(n-2) + ... (by the binomial theorem)
+# = 1 + an + 0. (remaining addends are 0 because they have a to the power of 2 or more, mod a^2)
+# Similarly for the other term, mod a^2:
+# (a-1)^n = (-1)^n + (n choose 1) a (-1)^(n-1) + ...
+# = (-1)^(n-1) (-1 + an + 0)
+# = if n is even then (1 - an) else (an - 1).
+# Therefore, adding the two terms:
+# (a+1)^n + (a-1)^n
+# = if n is even then 2 else 2an.
+#
+# We can always make 2an >= 2 by taking n=1, for example. So we can disregard the "n is even" case.
+# Maximizing 2an mod a^2 for n is the same as maximizing 2n mod a for n.
+# If a is even, then the maximum achievable value is a - 2 by setting n = a/2 - 1.
+# Else a is odd, then the maximum achievable value is a - 1 by setting n = (a - 1) / 2.
+#
+# In conclusion, if a is even, the maximum remainder is a(a-2);
+# otherwise a is odd, the maximum remainder is a(a-1).
+def compute():
+ ans = sum(i * (i - (2 if i % 2 == 0 else 1)) for i in range(3, 1001))
+ return str(ans)
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p123.py b/python/p123.py
new file mode 100644
index 0000000..b145af6
--- /dev/null
+++ b/python/p123.py
@@ -0,0 +1,32 @@
+#
+# Solution to Project Euler problem 123
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+# Using the result from Project Euler #120, we know that
+# (a-1)^n + (a+1)^n mod a^2 = if n is even then 2 else 2an.
+# Since 2 is tiny, we can skip the n is even case.
+# a is the n'th (1-based) prime number, so a > n. In fact for n >= 5,
+# we have a > 2n, so we can take 2an directly without moduloing it by a^2.
+def compute():
+ # Sieve of Eratosthenes
+ isprime = [True] * 1000001
+ isprime[0] = isprime[1] = False
+ for i in range(len(isprime)):
+ if isprime[i]:
+ for j in range(i * i, len(isprime), i):
+ isprime[j] = False
+ primes = [n for (n, x) in enumerate(isprime) if x]
+
+ for n in range(5, len(primes), 2):
+ rem = n * primes[n - 1] * 2
+ if rem > 10000000000:
+ return str(n)
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p124.py b/python/p124.py
new file mode 100644
index 0000000..0e9a7f6
--- /dev/null
+++ b/python/p124.py
@@ -0,0 +1,24 @@
+#
+# Solution to Project Euler problem 124
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+def compute():
+ # Modification of the Sieve of Eratosthenes
+ rads = [1] * 100001
+ for i in range(2, len(rads)):
+ if rads[i] == 1:
+ for j in range(i, len(rads), i):
+ rads[j] *= i
+
+ data = [(rads[i], i) for i in range(len(rads))]
+ data.sort()
+ return str(data[10000][1])
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p128.py b/python/p128.py
new file mode 100644
index 0000000..b83c1ae
--- /dev/null
+++ b/python/p128.py
@@ -0,0 +1,185 @@
+#
+# Solution to Project Euler problem 128
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+import eulerlib, itertools
+
+
+# Let's do mathematical analysis to drastically reduce the amount of
+# logic we need to implement and calculation the computer needs to do.
+# We begin with a couple of definitions.
+#
+# Ring number: Each cell belongs in a hexagonal ring,
+# numbered starting from 0 at the center like this:
+# 3
+# 3 3
+# 3 2 3
+# 3 2 2 3
+# 2 1 2
+# 3 1 1 3
+# 2 0 2
+# 3 1 1 3
+# 2 1 2
+# 3 2 2 3
+# 3 2 3
+# 3 3
+# 3
+#
+# Corner/edge cell: Within a ring, each cell is
+# either a corner cell or an edge cell, as shown:
+# C
+# E E
+# E C E
+# C E E C
+# C C C
+# E C C E
+# E C E
+# E C C E
+# C C C
+# C E E C
+# E C E
+# E E
+# C
+#
+# Basic observations:
+# - Except for the degenerate ring 0, each ring k has 6k cells.
+# The kth ring has exactly 6 corner cells and 6(k - 1) edge cells.
+# - In the code we will skip the PD (prime difference) calculation for
+# rings 0 and 1 because the existence of ring 0 breaks many patterns.
+# - Doing the PD calculation for rings 0 and 1 by hand (n = 1 to 7
+# inclusive), we find that PD(n) = 3 for and only for n = 1, 2.
+#
+# Now let's analyze the characteristics of all cells in rings 2 or above.
+# It's hard to justify these assertions rigorously, but they are true from
+# looking at the spiral diagram.
+#
+# - Corner cells along the upward vertical direction and the edge cells
+# immediately to the right of this vertical column are the most interesting,
+# so we will save these cases for last.
+#
+# - Claim: Except for cells immediately right of the upward corner column,
+# no edge cell satisfies PD(n) = 3. Proof: Take an arbitrary edge cell n
+# not immediately to the right of the upward corner column...
+# - The two neighbors in the same ring have a difference of 1 compared to n,
+# which is not a prime number.
+# - The two neighbors in the previous (inward) ring are consecutive numbers,
+# so exactly one of them has an even absolute difference with n. Because
+# n is in ring 2 or above, the difference with any neighboring number in the
+# previous ring is at least 6. Thus an even number greater than 2 is not prime.
+# - Similarly, the two neighbors in the next (outward) ring are consecutive numbers.
+# One of them has an even difference with n, and this number is also at least 6,
+# so one neighbor is definitely not prime.
+# - Therefore with at least 4 neighbors that do not have a prime difference, PD(n) <= 2.
+# Example of an edge cell n = 11 in ring 2, which is straight left of the origin:
+# 10
+# 24 03
+# 11
+# 25 04
+# 12
+#
+# - Claim: No corner cell in the other 5 directions satisfies PD(n) = 3.
+# Proof: Take an arbitrary corner cell n in the non-upward direction...
+# - Two of its neighbors (in the same ring) have a difference of 1,
+# which is not prime.
+# - One neighbor is in the previous ring (inward) while three neighbors
+# are in the next ring (outward).
+# - Let the inner ring neighbor be k and the outer ring's middle neighbor
+# be m. The three outer ring neighbors are {m - 1, m, m + 1}.
+# - Then n - k + 6 = m - n. Also, {m - 1, m + 1} have the same parity,
+# and {k, m} have the same other parity.
+# - Either both {|k - n|, |m - n|} are even or both {|m - 1 - n|, |m + 1 - n|} are even.
+# In any case, all these differences are at least 6, so the even numbers are not prime.
+# - Therefore with at least 4 neighbors that do not have a prime difference, PD(n) <= 2.
+# Example of a corner cell n = 14 in ring 2, which is straight below the origin:
+# 05
+# 13 15
+# 14
+# 28 30
+# 29
+#
+# - Now let's consider an arbitrary upward corner cell n in ring k, with k >= 2.
+# We shall give variables to all its neighbors like this:
+# d
+# e f
+# n
+# b c
+# a
+# - a is in the previous ring, {b, c} are in the same ring as n,
+# and {d, e, f} are in the next ring.
+# - Equations derived from the structure of the hexagonal spiral:
+# n = 3k(k - 1) + 2.
+# a = n - 6(k - 1).
+# b = n + 1.
+# c = n + 6k - 1 = d - 1.
+# d = n + 6k.
+# e = n + 6k + 1 = d + 1.
+# f = n + 6k + 6(k + 1) - 1 = n + 12k + 5.
+# - Hence we get these absolute differences with n:
+# |a - n| = 6(k - 1). (Not prime because it's a multiple of 6)
+# |b - n| = 1. (Not prime)
+# |c - n| = 6k - 1. (Possibly prime)
+# |d - n| = 6k. (Not prime because it's a multiple of 6)
+# |e - n| = 6k + 1. (Possibly prime)
+# |f - n| = 12k + 5. (Possibly prime)
+# - Therefore for each k >= 2, we need to count how many numbers
+# in the set {6k - 1, 6k + 1, 12k + 5} are prime.
+# Example of a corner cell n = 8 in ring 2, which is straight above the origin:
+# 20
+# 21 37
+# 08
+# 09 19
+# 02
+#
+# - Finally let's consider an arbitrary edge cell immediately to the right of the
+# upward vertical column. Suppose the cell's value is n and it is in ring k,
+# with k >= 2. Give variables to all its neighbors like this:
+# f
+# c e
+# n
+# a d
+# b
+# - {a, b} are in the previous ring, {c, d} are in the current ring, and {e, f} are in
+# the next ring. The ascending ordering of all these numbers is (a, b, c, d, n, e, f).
+# - Equations derived from the structure of the hexagonal spiral:
+# n = 3k(k + 1) + 1.
+# a = n - 6k - 6(k - 1) + 1 = n - 12k + 7.
+# b = n - 6k.
+# c = n - 6k + 1.
+# d = n - 1.
+# e = n + 6(k + 1) - 1 = n + 6k + 5.
+# f = n + 6(k + 1).
+# - Hence we get these absolute differences with n:
+# |a - n| = 12k - 7. (Possibly prime)
+# |b - n| = 6k. (Not prime because it's a multiple of 6)
+# |c - n| = 6k - 1. (Possibly prime)
+# |d - n| = 1. (Not prime)
+# |e - n| = 6k + 5. (Possibly prime)
+# |f - n| = 6(k + 1). (Not prime because it's a multiple of 6)
+# - Therefore for each k >= 2, we need to count how many numbers
+# in the set {6k - 1, 6k + 5, 12k - 7} are prime.
+# Example of an edge cell n = 19 in ring 2:
+# 37
+# 08 36
+# 19
+# 02 18
+# 07
+def compute():
+ TARGET = 2000 # Must be at least 3
+ count = 2 # Because n = 1 and 2 satisfy PD(n) = 3
+ for ring in itertools.count(2):
+ if eulerlib.is_prime(ring * 6 - 1) and eulerlib.is_prime(ring * 6 + 1) and eulerlib.is_prime(ring * 12 + 5):
+ count += 1
+ if count == TARGET:
+ return str(ring * (ring - 1) * 3 + 2)
+ if eulerlib.is_prime(ring * 6 - 1) and eulerlib.is_prime(ring * 6 + 5) and eulerlib.is_prime(ring * 12 - 7):
+ count += 1
+ if count == TARGET:
+ return str(ring * (ring + 1) * 3 + 1)
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p129.py b/python/p129.py
new file mode 100644
index 0000000..d83cc50
--- /dev/null
+++ b/python/p129.py
@@ -0,0 +1,58 @@
+#
+# Solution to Project Euler problem 129
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+import itertools
+
+
+# Let n >= 1 be arbitrary but assume that it is coprime with 10.
+# We want to find the smallest k such that R(k) = 0 mod n, and we'll show that 1 <= k <= n.
+#
+# Let "the sequence" of n values be (R(1) mod n, R(2) mod n, R(3) mod n, ..., R(n) mod n).
+# For the sake of contradiction, assume that none of the values in the sequence are 0.
+#
+# Each number in the sequence is an integer in the range [1, n).
+# The range has n - 1 elements, but there are n elements in the sequence.
+# Hence by the pigeonhole principle, there exist two distinct indexes
+# in the sequence where the elements have the same value.
+#
+# Suppose the two distinct indexes (1-based) are i and j.
+# So the two values in question are R(i) mod n and R(j) mod n.
+# Suppose WLOG that j > i. Then clearly R(j) - R(i) = 0 mod n,
+# and so R(j) - R(i) = 1...10...0 = R(j - i) * 10^i = 0 mod n.
+#
+# Since 10 is coprime with n, 10 (and its powers) are invertible modulo n.
+# Multiply everything in the equation by 10^-i, and we get R(j - i) = 1...1 = 0 mod n.
+#
+# We know 1 <= j - i <= n - 1. Then R(i - j) mod n, which is 0, is in the sequence.
+# This contradicts our assumption that none of (R(1), R(2), ... R(n)) is 0 mod n.
+#
+# Therefore if we want to find an n whose solution k is such that
+# k > 1000000, then we need to have n > 1000000.
+def compute():
+ LIMIT = 10**6
+ for n in itertools.count(LIMIT):
+ if least_divisible_repunit(n) > LIMIT:
+ return str(n)
+
+
+# Returns the smallest k such that R(k) is divisible by n.
+def least_divisible_repunit(n):
+ if n % 2 == 0 or n % 5 == 0:
+ return 0
+ k = 1
+ s = 1 # Loop invariant: Equal to R(k) mod n
+ p = 1 # Loop invariant: Equal to 10^k mod n
+ while s % n != 0:
+ k += 1
+ p = p * 10 % n
+ s = (s + p) % n
+ return k
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p132.py b/python/p132.py
new file mode 100644
index 0000000..16dd5a3
--- /dev/null
+++ b/python/p132.py
@@ -0,0 +1,29 @@
+#
+# Solution to Project Euler problem 132
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+import eulerlib
+
+
+def compute():
+ sum = 0
+ count = 0
+ i = 2
+ while count < 40:
+ if eulerlib.is_prime(i) and repunit_mod(10**9, i) == 0:
+ sum += i
+ count += 1
+ i += 1
+ return str(sum)
+
+
+def repunit_mod(k, m):
+ return (pow(10, k, m * 9) - 1) // 9
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p133.py b/python/p133.py
new file mode 100644
index 0000000..e6d8e90
--- /dev/null
+++ b/python/p133.py
@@ -0,0 +1,64 @@
+#
+# Solution to Project Euler problem 133
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+# Repunit formula: R(k) = (10^k - 1) / 9. (Using geometric series)
+#
+# For the rest of the argument, let n be an arbitrary integer that is coprime with 10.
+#
+# Let k = A(n) be the smallest positive integer such that R(k) = 0 mod n.
+# From problem #129, we know k exists and satisfies 1 <= k <= n.
+#
+# Lemma: For each natural number m, R(m) = 0 mod n if and only if m is a multiple of k.
+# Proof:
+# Backward direction:
+# Assume m is a multiple of k. Then factorize m = jk, where j is an integer.
+# Look at R(m) = R(jk) = 1...1 ... 1...1 (j groups of k 1's) = 10...010...010...01 * R(k) (informally)
+# = (sum of 10^(ik) for i = 0 to s-1) * R(k).
+# We already have R(k) = 0 mod n, thus (sum of 10^(ik) for i = 0 to s-1) * R(k) = R(m) = 0 mod n.
+# Forward direction (by converse):
+# Assume m is not a multiple of k. Suppose for contradiction that R(m) = 0 mod n.
+# Similar the previous argument, we can zeroize blocks of k 1's while preserving the value of R(m) mod n.
+# Namely, we delete the top k 1's by subtracting R(k) * 10^(m-k), which is 0 mod n because R(k) = 0 mod n.
+# After repeated deletion of the most significant 1's, we can get m' = m mod k, so that 0 < m' < k.
+# (m' != 0 because we assumed m is not a multiple of k.) But with R(m') = R(m) = 0 mod n, and m' < k,
+# this contradicts the definition of k = A(n), the smallest value such that R(k) = 0 mod n.
+# Hence the supposition that R(m) = 0 mod n is false.
+#
+# Does there exist an x such that R(10^x) is a multiple of n? By the lemma, this is true if and only if
+# there exists an x such that 10^x is a multiple of k. This means k must be a product of 2's and 5's.
+#
+# Actually, we don't need to compute k = A(n) to perform this test. If k = 2^a * 5^b, then all sufficiently large
+# powers of 10 are a multiple of k. (If k has other prime factors, then no power of 10 is a multiple of k.)
+# We know 1 <= k < n, so in this problem 1 <= k < 10^5. For k in this range, the largest exponent among a and b is 16
+# (for the number 2^16 = 65536). (In general, the largest exponent is floor(log2(limit)); in this case limit = 10^5.)
+# So we only need to test if 10^16 is a multiple of k, equivalent to testing if R(10^16) is a multiple of n.
+def compute():
+ # Sieve of Eratosthenes
+ isprime = [True] * 100001
+ isprime[0] = isprime[1] = False
+ for i in range(len(isprime)):
+ if isprime[i]:
+ for j in range(i * i, len(isprime), i):
+ isprime[j] = False
+ primes = [n for (n, x) in enumerate(isprime) if x]
+
+ ans = 0
+ for p in primes:
+ if p == 2 or p == 5 or not has_divisible_repunit(p):
+ ans += p
+ return str(ans)
+
+
+# Tests whether there exists a k such that R(10^k) is a multiple of p
+def has_divisible_repunit(p):
+ return (pow(10, 10**16, p * 9) - 1) // 9 % p == 0
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p250.py b/python/p250.py
new file mode 100644
index 0000000..af4f7e3
--- /dev/null
+++ b/python/p250.py
@@ -0,0 +1,30 @@
+#
+# Solution to Project Euler problem 250
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+def compute():
+ # Use dynamic programming
+ MOD = 10**16
+ subsets = [0] * 250 # subsets[i] is {the number of subsets with sum equal to i mod 250} mod 10^16
+ subsets[0] = 1
+
+ for i in range(1, 250250 + 1):
+ temp = pow(i, i, 250)
+ newsubsets = list(subsets)
+ for j in range(len(subsets)):
+ k = (j + temp) % 250
+ newsubsets[k] = (subsets[j] + subsets[k]) % MOD
+ subsets = newsubsets
+
+
+ ans = (subsets[0] - 1) % MOD
+ return str(ans)
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p301.py b/python/p301.py
new file mode 100644
index 0000000..305fb4b
--- /dev/null
+++ b/python/p301.py
@@ -0,0 +1,53 @@
+#
+# Solution to Project Euler problem 301
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+# In a game of Nim where both players play with perfect strategy, if the current state is a collection (multiset) of piles
+# with sizes {n1, n2, ..., n_m}, then the current player will lose if and only if n1 XOR n2 XOR ... XOR n_m = 0.
+# In this problem, we specialize the condition to just n XOR 2n XOR 3n = 0.
+#
+# Facts:
+# 3n = 2n + n.
+# a ^ b = 0 iff a = b. (from digital logic)
+# a + b = (a ^ b) + ((a & b) << 1). (from digital logic)
+# Hence:
+# n ^ 2n ^ 3n = 0 (what we want)
+# iff n ^ 2n = 3n (from the second fact)
+# iff n ^ 2n = (n ^ 2n) + ((n & 2n) << 1) (from the third fact)
+# iff (n & 2n) << 1 = 0 (by cancelling on both sides)
+# iff n & 2n = 0 (left-shifting doesn't change zeroness)
+# iff the binary representation of n does not have consecutive '1' bits.
+#
+# How many binary strings of length i have no consecutive 1's?
+# First partition the set into strings that begin with a 0 and strings that begin with a 1.
+# For those that begin with a 0, the rest of the string can be any string of length i-1 that doesn't have consecutive 1's.
+# For those that begin with a 1, the rest of the string can be any string of length i-1 that begins with a 0 and doesn't have consecutive 1's.
+# Let numStrings(i, j) be the number of bit strings of length i that begin with the bit j and have no consecutive 1's. Then:
+# numStrings(1, 0) = 1. (base case)
+# numStrings(1, 1) = 1. (base case)
+# numStrings(i, 0) = numStrings(i-1, 0) + numStrings(i-1, 1). (for i >= 2)
+# numStrings(i, 1) = numStrings(i-1, 0). (for i >= 2)
+# This corresponds to a shifted Fibonacci sequence, because:
+# numStrings(i, 0) = numStrings(i-1, 0) + numStrings(i-2, 0). (substitute)
+# numStrings(1, 0) = 1. (base case)
+# numStrings(2, 0) = 2. (derived)
+# So numStrings(i, 0) = fibonacci(i + 1).
+# What we want is numStrings(30, 0) + numStrings(30, 1) = numStrings(31, 0) = fibonacci(32).
+#
+# Actually, that answer considers numbers in the range [0, 2^30), which is not exactly what we want.
+# According to the problem statement, we need to exclude 0 and include 2^30. But both are losing positions, so the adjustments cancel out.
+def compute():
+ a = 0
+ b = 1
+ for i in range(32):
+ a, b = b, a + b
+ return str(a)
+
+
+if __name__ == "__main__":
+ print(compute())
diff --git a/python/p429.py b/python/p429.py
new file mode 100644
index 0000000..28eb150
--- /dev/null
+++ b/python/p429.py
@@ -0,0 +1,54 @@
+#
+# Solution to Project Euler problem 429
+# by Project Nayuki
+#
+# http://www.nayuki.io/page/project-euler-solutions
+# https://github.com/nayuki/Project-Euler-solutions
+#
+
+
+# Let n be an arbitrary positive integer. Suppose n is factorized as p1^k1 * p2^k2 * ... * {p_m}^{k_m},
+# where the p's are prime and distinct (thus the k's are as large as possible).
+# Let {p1^k1, p2^k2, ..., {p_m}^{k_m}} be the set of "maximal prime powers".
+# Then all the unitary divisors of n are exactly all the subsets of maximal prime powers,
+# where each subset is viewed as a product of its elements.
+#
+# For n!, its prime factorization uses and only uses all prime numbers from 1 to n (inclusive).
+# For each prime p, the number n! has exactly floor(n/p) + floor(n/p^2) + floor(n/p^3) + ... factors of p.
+# Thus we can calculate the p's and k's quite easily.
+#
+# To solve the remaining parts of the problem, we use dynamic programming.
+# Suppose we have found all the unitary divisors that are products of maximal prime powers less than {p_i}^{k_i},
+# and suppose this set is {a, b, c}. Then when we include {p_i}^{k_i} into consideration, we double the size of the set
+# because now {a * {p_i}^{k_i}, b * {p_i}^{k_i}, c * {p_i}^{k_i}} are also unitary divisors.
+def compute():
+ LIMIT = 10**8
+ MOD = 1000000009
+
+ # Sieve of Eratosthenes
+ isprime = [True] * (LIMIT + 1)
+ isprime[0] = isprime[1] = False
+ for i in range(len(isprime)):
+ if isprime[i]:
+ for j in range(i * i, len(isprime), i):
+ isprime[j] = False
+ primes = [n for (n, x) in enumerate(isprime) if x]
+
+ ans = 1
+ for p in primes:
+ power = count_factors(LIMIT, p)
+ ans *= 1 + pow(p, power * 2, MOD)
+ ans %= MOD
+ return str(ans)
+
+
+# Returns the number of factors of p (prime) in factorial(n).
+def count_factors(n, p):
+ if n == 0:
+ return 0
+ else:
+ return n // p + count_factors(n // p, p)
+
+
+if __name__ == "__main__":
+ print(compute())
--
2.21.0