Haskell(六) Project Euler 练习1-26

1. Haskell(一)入门
2. Haskell(二)函数式编程
3. Haskell(三) Monad
4. Haskell(四)总结和工具链
5. Haskell(五) 总结和展望
6. Haskell(六) Project Euler 练习1-26

本文会选择一些有意思的 Project Euler 的题目，学习怎么用 Haskell 写算法，并且逐渐学习相关语言和数学知识。前 100 题可以直接分享答案，后续的题目就只给暗示和加密的答案，可以发邮件获取解密方式。

代码仓库：https://github.com/learnerLj/projecteuler

P1

显然 3 和 5 的倍数，其实是 3 和 5 的倍数，去除重复算的 15 的倍数。而且 3 的倍数的和，还可以化简成求和公式。

sumOfMultiples limit = sumDivisible 3 + sumDivisible 5 - sumDivisible 15
    where
        cumulativeSum n = n * (n + 1) `div` 2
        sumDivisible n = n * cumulativeSum (limit `div` n)

p2 偶数项斐波那契数列

https://projecteuler.net/problem=2

偶数项的斐波那契数列，实际上也有递推公式。

sumEvenLessThanImproved :: Integer -> Integer
sumEvenLessThanImproved limit = sum $ takeWhile (< limit) evenFibList
  where
    evenFibList = 2 : 8 : zipWith (\x y -> x + 4 * y) evenFibn1 evenFibn2
    evenFibn1 = evenFibList
    evenFibn2 = tail evenFibList

或者先求出来斐波那契数列，然后筛选偶数项也行。

sumEvenLessThan :: Integer -> Integer
sumEvenLessThan limit = sum $ filter even $ takeWhile (< limit) fibList
  where
    fibList =
      let fibn1 = fibList
          fibn2 = tail fibList
       in 1 : 2 : zipWith (+) fibn1 fibn2

P3 质因数分解

https://projecteuler.net/problem=3

质因数分解就先求出来质数，然后再一个一个匹配，看是否是某个数的因数。

primes :: [Integer]
primes = sieve [2..]
    where
    sieve [] = []
    sieve (p:xs) = p: sieve [x| x<-xs, x `mod` p /=0]

primeFactors :: Integer -> [Integer]
primeFactors n =primeFactors' n primes
    where
        primeFactors' 1 _ = []
        primeFactors' _ [] = []
        primeFactors' m (p:xs)
            | m `mod` p ==0 = p: primeFactors' (divide m p) xs
            | otherwise = primeFactors' m xs
        divide m p
            | m `mod` p ==0 = divide (m `div` p) p
            | otherwise = m

largestPrimeFactor :: Integer -> Integer
largestPrimeFactor n= last $ primeFactors n

P4 回文数的分解

这里学习如何用 Haskell 写循环，对于不熟悉的人来说，还是有点挑战的。首先遍历 x 和 y，y 从大到小，如果 y 到达下界，就 x-1，重置 y，直到 x 到达下界返回。

这里回文数比如整除 11，然后 11 是质数，所以回文数之一必然整除 11，所以每次 y -11。

isPalindrome :: Int -> Bool
isPalindrome n = n == reverseNum n
    where
    reverseNum :: Int -> Int
    reverseNum = go 0
        where
        go acc 0 = acc
        go acc x = let (q, r) = x `quotRem` 10 in go (acc * 10 + r) q


largestPalindromeProduct :: Int -> Int
largestPalindromeProduct n = go 0 upperBound (divisible11 upperBound)
  where
    lowerBound = 10 ^ (n - 1)
    upperBound = 10 ^ n - 1
    go maxPalin x y
      | x < lowerBound = maxPalin
      | y < lowerBound = go maxPalin (x - 1) (divisible11 (x-1))
      | isPalindrome currentPalin = go (max maxPalin currentPalin) (x-1) (divisible11 (x-1))
      | otherwise = go maxPalin x (y - 11)
      where currentPalin = x * y
    divisible11 m = m - m `mod` 11

P5 最小公倍数

最小公倍数 x，事实上是覆盖 a 和 b 的质数因子的最小数。比如 6=2*3, 8= 2*2*2，那么最小公倍数要有 3 个 2，1 个 3，也即 2*2*2*3=24。多个数也是一样的。

lcmRangeImproved :: Integer -> Integer
lcmRangeImproved n = product $ map maxPower primesInRange
  where
    primesInRange = takeWhile (<= n) primes
    maxPower p = p ^ maxExponent p
    maxExponent p = last $ takeWhile (\x -> p ^ x <= n) [1 ..]
primes :: [Integer]
primes = sieve [2 ..]
  where
    sieve [] = []
    sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p /= 0]

P9 毕达哥拉斯三元组

毕达哥拉斯三元组，满足特殊代数关系，类似平方和。另外，三元组的倍数仍然是三元组。所以有如下数学推导，三边满足

$$a=m^2-n^2, b=2mn, c=m^2+n^2\\ \mathrm{where}~m>n$$

那么 $a+b+c=2m(m+n)$，所以 $m$ 和 $n$ 实际上只有很少的选择，再然后考虑去除倍数的影响，确保 $m$ 和 $m+n$ 互素。

findPythagoreanTripletImproved :: Integer -> Maybe [(Integer, Integer, Integer)]
findPythagoreanTripletImproved totalSum
    | odd totalSum = Nothing
    | null candicates = Nothing
    | otherwise = Just $ map getSides candicates
  where
    candicates =
        let relativePrimeFac (x, y) = (x `div` d, y `div` d, d ^ 2)
              where
                d = gcd x y
         in [relativePrimeFac (a, b) | (a, b) <- factors $ totalSum `div` 2, b < 2 * a]
    getSides (m, k, d) = (d * a, d * b, d * c)
      where
        n = k - m
        a = m ^ 2 - n ^ 2
        b = 2 * m * n
        c = m ^ 2 + n ^ 2
factors :: Integer -> [(Integer, Integer)]
factors x = [(a, x `div` a) | a <- [2 .. (floor . sqrt $ fromInteger x)], x `mod` a == 0]

P12 三角形数的因子个数

三角形数有特征的，满足$n(n+1)/2$，然后因子个数实际是质因子的组合，对于每类质因子$a$，如果有 $m_a$个，那就可以选择 0 到 $m_a$ 个 $a$ 组成一个因子，所以 $m+1$ 种选择。为了复用质因数列表，于是用了 State Monad。

type PrimeState = State [Int]
firstDivsorsN2 :: Int -> Int
firstDivsorsN2 num = evalState (findFirstDivisor num) (sieve [2 ..])
  where
    findFirstDivisor :: Int -> PrimeState Int
    findFirstDivisor n = head <$> filterM (fmap (> n) <$> divisorNum) triNum

    divisorNum :: Int -> PrimeState Int
    divisorNum m = product <$> (map (succ . length) . group <$> primeFactors m)

    primeFactors :: Int -> PrimeState [Int]
    primeFactors num = gets $ \ps -> factorize num ps
      where
        factorize 1 _ = []
        factorize _ [] = []
        factorize m (p : ps)
          | m `mod` p == 0 = p : factorize (m `div` p) (p : ps)
          | otherwise = factorize m ps

    sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p /= 0]
    triNum = [m * (m + 1) `div` 2 | m <- [1 ..]]

P14 最长 Collatz 序列

序列长度的递推关系很清晰，但是最好要记录下状态，可以在已有的状态上递推，所以用 State Monad。之后比较长度即可。

{-# LANGUAGE TupleSections #-}

module P14 (longestUnder2) where

import Control.Monad.State
import Data.List (foldl', maximumBy)
import qualified Data.Map as Map
import Data.Ord (comparing)

-- | Finds the starting number under a given limit which produces the longest Collatz sequence.
--   'n' is the upper limit for the starting number.
longestUnder :: Integer -> Integer
longestUnder n = fst $ maximumBy (comparing snd) $ zip [1 .. n] (map collatzLen [1 .. (n - 1)])
  where
    -- \| Computes the length of the Collatz sequence for a given starting number.
    collatzLen 1 = 1 -- Base case: sequence length is 1 for starting number 1.
    collatzLen start
      | even start = 1 + collatzLen (start `div` 2) -- If 'start' is even.
      | otherwise = 1 + collatzLen (3 * start + 1) -- If 'start' is odd.

-- Type aliases for readability.
type Cache = Map.Map Int Int

type CollatzState = State Cache

-- | An optimized version of 'longestUnder' using the State Monad for caching.
longestUnder2 :: Int -> Int
longestUnder2 m = fst $ maximumBy (comparing snd) $ evalState (mapM collatzPair [1 .. m - 1]) Map.empty
  where
    -- \| Creates a pair (starting number, sequence length) for each number.
    collatzPair :: Int -> CollatzState (Int, Int)
    collatzPair num = (num,) <$> collatzLen num

    -- \| Computes the Collatz sequence length with caching.
    collatzLen :: Int -> CollatzState Int
    collatzLen 1 = return 1 -- Base case: sequence length is 1 for starting number 1.
    collatzLen n = do
      cache <- get
      maybe (updateCache n) return (Map.lookup n cache)
      where
        -- \| Updates the cache with the new sequence length if not already cached.
        updateCache num = do
          len <- collatzLen $ if even num then num `div` 2 else 3 * num + 1
          modify (Map.insert num (1 + len))
          return (1 + len)

-- | Main function to print the starting number under 1,000,000 with the longest Collatz sequence.
main :: IO ()
main = print $ longestUnder2 1000000

P15 走格子路径的方式

latticePaths2 side = dp !! side !! side
  where
    dp = [[lattice x y | y <- [0 .. side]] | x <- [0 .. side]]
    lattice 0 _ = 1
    lattice _ 0 = 1
    lattice x y = dp !! (x - 1) !! y + dp !! x !! (y - 1)

latticePaths3 :: (Integral a) => a -> a
latticePaths3 side = numerator $ product [(side + i) % i | i <- [1 .. side]]

可以直接看代码，每个点的来源是左边或者上面的点，所以动态规划就可以了。不过呢，其实也是有组合公式，因为每个点的路径，实际上是 d(down)r(right)，这样 ddrddrr 的序列组合，这样就到达了(4,3)。7 个位置选出 3 个给 r 就是一种走法。

P18 路径之和

step 函数是从最后一行和倒数第二行开始，计算得到从最后一行到倒数第二行的每个数字最大的和，然后这个列表继续和倒数第三行操作，最后就剩下一个数。记 $s_{i,j}$ 为从上到下到第 $i$ 行第 $j$ 列的数 $n_{i,j}$ 的路径之和。那么

$s_{i,j}=\max{ s_{i-1,j-1}, s_{i-1,j}}+n_{i,j}$，对了如果下标超出范围，值都是 0。

-- | Computes the maximum path sum in a triangle.
--   The triangle is represented as a list of lists, where each list is a row in the triangle.
maxPathSum :: [[Int]] -> Int
maxPathSum [] = 0 -- Base case: if the triangle is empty, the max sum is 0.
maxPathSum tri = head $ foldr1 step tri
  where
    -- \| 'step' combines two rows of the triangle into one by choosing the maximum adjacent numbers.
    --   It is used as the function in a right fold to accumulate the maximum path sum.
    step [] _ = [] -- Base case for the step function when the last row is empty.
    step _ [] = [] -- Base case for the step function when the next row is empty.
    step lastRow@(x : xs) nextRow@(y : ys) =
      (x + max y (head ys)) : step xs ys -- Calculate the max sum for each position in the row.

-- 'foldr1' is used to apply 'step' function starting from the bottom of the triangle,
-- eventually reducing the entire triangle to a single row containing only the max path sum.

Problem 21 - 40

P21 Amicable Numbers

这里实际想要证明，因子之和，其实有一个公式。显然对于质数 $p$ 有：

$$\sigma \left( p \right) =p+1$$

对于质因子只有 $p$ 的数 $P^a$有

$$\sigma \left( p^a \right) =\sum_{i=0}^a{p^i}=\frac{p^{a+1}-1}{p-1}$$

这是因为因子只能是几个 1 或者 p 的倍数。

考虑如果一个数有 2 个不同质因数 p 和 q 组成，分别由 a 和 b 个，那么质因数组成的方式是 p 选出若干个，q 选出若干个，实际上就是下面多项式相乘的方式

$$\left( 1+p+\cdots p^a \right) \left( 1+q+\cdots q^b \right) =\sigma \left( p^a \right) \sigma \left( q^b \right)$$

那么拓展到若干个不同质因数的情况，知道是乘起来啦。

下面代码中，关键在于求真因子的和，然后记忆化了这个和。

module P21 (sumAcicable, sumAcicable2) where

import Control.Monad.State
import Data.List (group, nub)
import qualified Data.Map as Map

-- | Calculates the sum of amicable numbers under a given limit.
sumAcicable :: Int -> Int
sumAcicable num = sum $ acicable (num - 1)
  where
    -- \| Generates a list of amicable numbers up to a given limit.
    acicable n = [a | a <- [2 .. n], let b = sumDivisors a, a == sumDivisors b, a /= b]

    -- \| Calculates the sum of proper divisors of a number.
    sumDivisors n = foldr (\(x, y) acc -> acc + x + y) 1 (factors n)

    -- \| Generates a list of factors of a number along with their complementary factors.
    factors x = [(a, if a * a == x then 0 else x `div` a) | a <- [2 .. (floor . sqrt $ fromIntegral x)], x `mod` a == 0]

-- | An alternative implementation using the State Monad for caching.
sumAcicable2 :: Int -> Int
sumAcicable2 num = evalState (acicable (num - 1)) Map.empty
  where
    -- \| Generates a list of prime numbers up to half the given limit.
    primes = sieve [2 .. num `div` 2]

    -- \| Calculates the sum of amicable numbers using a State Monad for caching.
    acicable :: Int -> State (Map.Map Int Int) Int
    acicable n = sum <$> mapM evaluateDivisors [2 .. n]

    -- \| Evaluates and returns an amicable number if conditions are met, using cached results.
    evaluateDivisors :: Int -> State (Map.Map Int Int) Int
    evaluateDivisors a = do
      b <- sumDivisors a
      bDivSum <- sumDivisors b
      if bDivSum == a && a /= b then return a else return 0

    -- \| Computes the sum of divisors of a number, with caching.
    sumDivisors :: Int -> State (Map.Map Int Int) Int
    sumDivisors n = do
      cache <- get
      case Map.lookup n cache of
        Just s -> return s
        Nothing -> do
          let facs = primeFactors n primes
              exponents = map length $ group facs
              deltas = zipWith (\p a -> (p ^ (a + 1) - 1) `div` (p - 1)) (nub facs) exponents
              s = product deltas - n
          modify (Map.insert n s)
          return s

    -- \| Generates the prime factors of a number.
    primeFactors :: Int -> [Int] -> [Int]
    primeFactors 1 _ = []
    primeFactors _ [] = []
    primeFactors m (p : ps)
      | m `mod` p == 0 = p : primeFactors (m `div` p) (p : ps)
      | otherwise = primeFactors m ps

    -- \| Generates a list of prime numbers using the sieve method.
    sieve :: [Int] -> [Int]
    sieve [] = []
    sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p /= 0]

-- | Main function to print the sum of amicable numbers under 10000.
main = print $ sumAcicable2 10000

P26 小数循环节和循环群

循环群的角度

分母和 10 互质时

做除法的过程，实际上是计算 $10^k \mod n$ 的值的过程，这个值就是作为第 k 次除法的余数。当余数为 1 的时候，也即$10^k \mod n=1$， 说明就构成了循环。但是这不一定是第一次循环，因为余数可能在其他数字，比如 2 的时候，出现过 2 次，那么这个循环就要更短。

先考虑当 10 和 n 互质的时候，构成循环群$\{10^k\mathrm{mod}n|k\in \mathbb{N} ^*\}$，而且那么循环群的阶是最小使它余 1 的数，而且必然是 $\varphi \left( n \right)$的因子。显然如果$\varphi \left( n \right)$是质数，那么阶就是$\varphi \left( n \right)$。

这是因为群是闭合的，任意两个群元素的乘积仍然是群的元素。如果我们考虑 $a^k \mod m$ 的所有可能的幂，这个集合在模 $m$ 的乘法下也是闭合的。那么如果 $a$ 生成了一个循环群，那么这个群的阶（即群中元素的个数）就是最小的正整数 $k$，使得 $a^k \equiv 1 \pmod{m}$。这个阶 $k$ 必须整除任何 $a^x \equiv 1 \pmod{m}$ 成立的 $x$，包括 $\phi(m)$

由于 $a^{\phi(m)} \equiv 1 \pmod{m}$ 总是成立的（根据欧拉定理），这意味着群的阶必须是 $\phi(m)$ 的因子。因为，如果存在某个 $k < \phi(m)$，使得 $a^k \equiv 1 \pmod{m}$，那么 $\phi(m)$ 必须是 $k$ 的倍数，以保证群的循环性和阶的定义。

简而言之，幂的循环至多在达到 $\phi(m)$ 时完成一次循环，而实际的循环可能在更早的点完成，因此实际的循环长度（群的阶）是 $\phi(m)$ 的因子。

计算阶的办法 除了直接从 1 开始遍历，找到最小正整数使得 $10^k \mod m=1$，还可以进行简化，尤其是 k 非常大的情况。
当$m = pq$，其中$p$和$q$是不同的质数，并且$a$与$p$和$q$都互质。

假设$a^k \mod p$的乘法群阶为$n_p$，$a^k \mod q$的乘法群阶为$n_q$。根据欧拉定理，我们知道$a^{\phi(p)} \equiv 1 \mod p$和$a^{\phi(q)} \equiv 1 \mod q$，因此$n_p | \phi(p)$且$n_q | \phi(q)$。

$a^k \mod m$的乘法群阶，设为$n_m$，是最小的正整数$n$，使得$a^n \equiv 1 \mod m$。由于$m = pq$，根据中国剩余定理，如果我们可以同时满足$a^n \equiv 1 \mod p$和$a^n \equiv 1 \mod q$，那么$a^n \equiv 1 \mod m$也成立。

为了找到$n_m$，我们需要找到最小的$n$，使得$a^n \equiv 1 \mod p$且$a^n \equiv 1 \mod q$。这意味着$n$必须同时是$n_p$和$n_q$的倍数。因此，$n_m = \text{lcm}(n_p, n_q)$。

那就到了另外一个问题了，如何快速求得 $10^k \mod p^b$ 的阶？我们只知道阶一定是 $\phi(p^b)=p^{b-1}(p-1)$ 的因子。

分母和 10 不互质的时候

当 10 和 n 不互质的时候，n 的因数包含 2 或者 5，设 $n=2^a*5^b*m$，m 与 10 互质，a,b 是自然数且不全为 0。此时

$$10^k \mod (2^a5^b)$$

当 k 足够大时，必然为 0。我们只需要证明，有限小数乘以无限循环小数，不会改变循环节长度。这个需要证明，但是我不会。

也就是不必管 2 或者 5 作为因子的情况，提取出与 10 互质的因子即可。（证明我还不会）。

算法

首先除净 2 和 5，值一定是$p^aq^b$这样结构组合（部分指数可以为 0），此时$\phi(p^b)$ 小于 $p^b$，那么计算最小公倍数，肯定小于$p^a q^b-1$。也就是说，对于质数 p，它的阶必然大于任何比他小的数的阶。所以只要找到 1000 以内最大的质数 $p_{max}$，并且和直到 1000，除净之后且比$p_{max}$大的数比较即可。

这又涉及到问题**相邻两个质数 p, q 之间的数，除以 2 之后，可能比 p 大吗？**这实际上等价于是否存在$p<x<q$使得 $2p<x<q$。这我不知道，但是在 1000 以内是不可能的。

那么对于**两个相邻质数 p,q 之间，存在$x=p_{1}^{a_1}p_{2}^{a_2}\cdots p_{m}^{a_m}$不是 2 或 5 的倍数，也不是质数，但是 10^k mod x 的阶比 10^k mod p 的阶大吗？**也就是

$$\mathrm{gcd}\left( \varphi \left( p_{1}^{a_1} \right) ,\cdots ,\varphi \left( p_{m}^{a_m} \right) \right)$$

可能比 $p-1$ 大吗？

这个数必然小于这些欧拉函数的乘积

$$\varphi \left( x \right) =\varphi \left( p_{1}^{a_1} \right) \cdots \varphi \left( p_{m}^{a_m} \right)$$

且$\varphi \left( x \right)<x$。而且除了质数 2 之外，p-1 必然是偶数，那么必然存在公因数 2，回到上一个问题，至少在 1000 以内是不可能。在除以 2 之后，一定比 $p$小。

所以，循环节最长的数，一定是质数。假如能证明：

1. 相邻两个质数 p, q 之间的数，除以 2 之后，一定比 p 小。
2. 有限小数乘以无限循环小数，不会改变循环节长度。

这似乎已经有相关定理，可以证明第二点，而且可以得到非循环部分的长度。

常规直观做法

循环小数一定是无限的，而且循环节中的数字可能重复。按照编码方面的经验，不读取到最后一个数字，是无法直到这一串数字，是否是循环节。所以必须设定一个最大可能的循环节长度。

其次，即使是循环节，那么不一定是最短的循环节，可能是最短循环节的倍数。比如说 001 循环，也会满足 001001 的循环。所以当找到一个循环节时，它的因子的长度，不能是循环节。

循环的部分，从任意地方选区循环节的长度，它的周期不变。就像 001 循环，实际上只看后半部分，也可以看作时 010 的循环。

所以，只要跳过足够的小数位，就一定是循环的部分，并且周期不变。那只要按着周期的长度，去找两段匹配上即可。

-- Generates the decimal representation of 1/deno after the decimal point.
afterPoint :: Int -> [Int]
afterPoint = afterPoint' 1
  where
    afterPoint' num deno =
      let (d, r) = (10 * num) `divMod` deno -- Multiplies numerator by 10, divides by denominator, and keeps remainder.
       in d : afterPoint' r deno -- Recursively continues with remainder as new numerator.

-- Calculates the length of the recurring cycle in the decimal representation of 1/deno.
recurCycle :: Int -> Int
recurCycle deno = recurCycle' onlyCycle maxPossibleLen
  where
    maxSkip = myLog2 deno -- Determines the number of digits to skip based on log base 2 of denominator.
    onlyCycle = drop maxSkip $ afterPoint deno -- Skips initial digits that don't participate in the cycle.
    maxPossibleLen = deno - 1 -- The maximum possible cycle length is deno - 1.
    recurCycle' xs len
      | all (== 0) (take maxPossibleLen xs) = 0 -- If all digits are zero, cycle length is 0.
      | whenLen len && not (any whenLen (factors len)) = len -- Checks if current length satisfies cycle conditions.
      | otherwise = recurCycle' xs (len - 1) -- Otherwise, decreases length and tries again.
      where
        whenLen l = (l /= 0) && listEq (take l xs) (drop l xs) -- Checks if two subsequences of length l are equal.
    myLog2 :: Int -> Int
    myLog2 num = floor $ logBase 2 (fromIntegral num) -- Calculates the floor of log base 2 of a number.
    listEq [] _ = True -- Base case for equality check: empty list is equal to any list.
    listEq (x : xs) (y : ys) = (x == y) && listEq xs ys -- Recursive case: checks if heads are equal and proceeds to tails.

-- Generates the list of factors of a number by finding divisors up to the square root and their complements.
factors :: Int -> [Int]
factors num = let halfFactors = filter ((== 0) . mod num) [2 .. floor . sqrt $ fromIntegral num] in halfFactors ++ map (div num) (reverse halfFactors)

考虑循环节长度特征的做法

循环小数的循环部分，实际上是一个数乘以几何级数，可以记做

$$0.\left( a_1\cdots a_n \right) =a_1\cdots a_n\sum_{i=1}^{\infty}{10^{-in}}$$

根据求和公式，可以得到

$$0.\left( a_1\cdots a_n \right) =\frac{a_1\cdots a_n\left( 10^{-n}-\lim_{i\rightarrow +\infty} 10^{-in} \right)}{1-10^{-n}} \\ =\frac{a_1\cdots a_n\left( 1-\lim_{i\rightarrow +\infty} 10^{-in} \right)}{10^n-1}$$

观察到，循环节的长度 n 和分母有关，或者说由分母决定。只要能把分母转化成$10^n-1$的形式，那么一定是循环小数。也就是说，对于分数 $a/b$，找到正整数$n$，使得$b|10^n-1$，那么这个 $n$ 就是循环节。

代码如下，需要避免整数溢出：

-- Improved version of recurCycle that utilizes an infinite list of numbers composed entirely of 9s to find the cycle length.
recurCycleImprove :: Integer -> Int
recurCycleImprove deno
  | any ((== 0) . mod deno) [2, 5] = 0 -- Returns 0 immediately if deno is a multiple of 2 or 5.
  | otherwise = maybe 0 (+ 1) (findIndex ((== 0) . (`mod` deno)) all9s) -- Finds the cycle length using an index in the all9s list.

-- Generates an infinite list of integers where each integer is composed entirely of 9s.
all9s :: [Integer]
all9s = map (\n -> 10 ^ n - 1) [1 :: Integer ..]

背景知识

欧拉定理描述

首先学习欧拉定理（数论），欧拉定理表明，若 $n,a$ 为正整数，且互质，那么有

$$a^{\varphi (n)}\equiv 1\left( \mathrm{mod}\;n \right)$$

其中 $\varphi \left( n \right)$ 表示小于 n 的正整数中和 n 互质的数，而且满足

$$\varphi \left( n \right) =n\prod{\left( 1-\frac{1}{p} \right)}$$

其中 p 取 n 的所有质因子。比如说$n=12=2^2\times 3$ ，互质的数为$\left\{ 3,5,7,11 \right\}$，恰好

$$12\times \left( 1-\frac{1}{2} \right) \times \left( 1-\frac{1}{3} \right) =4$$

欧拉定理证明

相关背景知识，可见：初等数论笔记 Part 1： 欧拉定理

欧拉函数的证明，分成各种情况：

素数幂

$n=p^a$ 显然只有 $kp$ 与 n 有除 1 以外的公因数，k 取 1 到 $p^{a-1}-1$ 共 $p^{a-1}$个数，所以

$$\varphi \left( p^a \right) =p^a-p^{a-1}=p^a\left( 1-\frac{1}{p} \right)$$

互质的整数的乘积

如果 $m$ 和 $n$ 是互质的整数，则 $\varphi(mn) = \varphi(m)\varphi(n)$。这是因为由于 $m$ 和 $n$ 互质，小于 $mn$ 的整数中与 $mn$ 互质的整数，可以通过取小于 $m$ 的与 $m$ 互质的整数，与小于 $n$ 的与 $n$ 互质的整数，的所有可能组合来构造。

我们需要考虑由 $m$ 和 $n$ 形成的整数对格点，并分析这些点中有多少是与 $mn$ 互质的。具体步骤如下：

1. 整数对的构造：
   考虑由 $m$ 和 $n$ 形成的 $mn$ 个整数对 $(i, j)$，其中 $1 \leq i \leq m$ 和 $1 \leq j \leq n$。如果$ij$ 与$mn$互质，称作互质点，这些整数对可以视为一个 $m \times n$ 的矩形格点图。
2. 互质点的计数：
   对于每个整数对 $(i, j)$，我们需要判断 $mn$ 是否与 $ij$ 互质。因为 $m$ 和 $n$ 互质，显然 $i$ 必须与$m$互质，$j$与 $n$互质，才会有 $ij$与$mn$互质，否则就存在公因数了。所以互质点最多$\varphi(m)\varphi(n)$个。
3. 独立互质点的乘积：
   在 $m$ 中与 $m$ 互质的整数数量为 $\varphi(m)$，在 $n$ 中与 $n$ 互质的整数数量为 $\varphi(n)$。这些互质的整数可以独立地在 $m$ 和 $n$ 的范围内选择，而且不会影响对方的选择，所以互质点最少有$\varphi(m)\varphi(n)$个。

因此，与 $mn$ 互质的整数对 $(i, j)$ 的总数就是 $\varphi(m)$ 和 $\varphi(n)$ 的乘积，即 $\varphi(m) \times \varphi(n)$。综上所述，我们得出 $\varphi(mn) = \varphi(m)\varphi(n)$。

一般情况

对于任意整数 n，存在唯一质数分解 $n = p_1^{k_1} p_2^{k_2} \cdots p_r^{k_r}$。由于这些素数幂相互之间是互质的，根据第 2 步，我们可以写出

$$\varphi(n) = \varphi(p_1^{k_1})\varphi(p_2^{k_2}) \cdots \varphi(p_r^{k_r})$$

根据第 1 步，我们知道每个 $\varphi(p_i^{k_i})$。将它们代入上面的公式，我们得到

$$\varphi(n) = p_1^{k_1}\left(1 - \frac{1}{p_1}\right) p_2^{k_2}\left(1 - \frac{1}{p_2}\right) \cdots p_r^{k_r}\left(1 - \frac{1}{p_r}\right)$$

这可以进一步简化为

$$\varphi(n) = n\left(1 - \frac{1}{p_1}\right)\left(1 - \frac{1}{p_2}\right) \cdots \left(1 - \frac{1}{p_r}\right)$$

循环群

循环群是指能由单个元素所生成的群。设$(G, \cdot )$为一个群，若存在一个元素 ，使得 $G=\left\langle \,g\,\right\rangle =\left\{g^{k}|\;k\in \mathbb {Z} \right\}$，那么$(G, \cdot )$ 形成循环群。群 $G$ 内任意一个元素所生成的群都是循环群，而且是 $G$ 的子群。

假设 $(G,*)$ 是一个 群（group），若 $H$ 是 $G$ 的一个非空子集（subset）且同时 $H$ 与相同的二元运算 $*$ 亦构成一个群，则 $(H,*)$ 称为 $(G,*)$ 的一个 子群（subgroup）。