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

image-20240114210638244

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

1
2
3
4
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

image-20240114210657311

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

1
2
3
4
5
6
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

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

1
2
3
4
5
6
7
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

image-20240114210721042

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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 回文数的分解

image-20240114210737606

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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 最小公倍数

image-20240114211157599

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

1
2
3
4
5
6
7
8
9
10
11
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 毕达哥拉斯三元组

image-20240114211714499

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

那么 ,所以 实际上只有很少的选择,再然后考虑去除倍数的影响,确保 互素。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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 三角形数的因子个数

image-20240114213717924

三角形数有特征的,满足,然后因子个数实际是质因子的组合,对于每类质因子,如果有 个,那就可以选择 0 到 组成一个因子,所以 种选择。为了复用质因数列表,于是用了 State Monad。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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 序列

image-20240114214356024

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
{-# 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 走格子路径的方式

image-20240114214650028

1
2
3
4
5
6
7
8
9
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 路径之和

image-20240114215043832

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

,对了如果下标超出范围,值都是 0。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
-- | 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

image-20240114215746348

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

对于质因子只有 的数

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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 小数循环节和循环群

image-20240122011802524

循环群的角度

分母和 10 互质时

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

先考虑当 10 和 n 互质的时候,构成循环群,而且那么循环群的阶是最小使它余 1 的数,而且必然是 的因子。显然如果是质数,那么阶就是

这是因为群是闭合的,任意两个群元素的乘积仍然是群的元素。如果我们考虑 的所有可能的幂,这个集合在模 的乘法下也是闭合的。那么如果 生成了一个循环群,那么这个群的阶(即群中元素的个数)就是最小的正整数 ,使得 。这个阶 必须整除任何 成立的 ,包括

由于 总是成立的(根据欧拉定理),这意味着群的阶必须是 的因子。因为,如果存在某个 ,使得 ,那么 必须是 的倍数,以保证群的循环性和阶的定义。

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

计算阶的办法 除了直接从 1 开始遍历,找到最小正整数使得 ,还可以进行简化,尤其是 k 非常大的情况。
,其中是不同的质数,并且都互质。

假设的乘法群阶为的乘法群阶为。根据欧拉定理,我们知道,因此

的乘法群阶,设为,是最小的正整数,使得。由于,根据中国剩余定理,如果我们可以同时满足,那么也成立。

为了找到,我们需要找到最小的,使得。这意味着必须同时是的倍数。因此,

那就到了另外一个问题了,如何快速求得 的阶?我们只知道阶一定是 的因子。

分母和 10 不互质的时候

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

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

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

算法

首先除净 2 和 5,值一定是这样结构组合(部分指数可以为 0),此时 小于 ,那么计算最小公倍数,肯定小于。也就是说,对于质数 p,它的阶必然大于任何比他小的数的阶。所以只要找到 1000 以内最大的质数 ,并且和直到 1000,除净之后且比大的数比较即可。

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

那么对于**两个相邻质数 p,q 之间,存在不是 2 或 5 的倍数,也不是质数,但是 10^k mod x 的阶比 10^k mod p 的阶大吗?**也就是

可能比 大吗?

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

。而且除了质数 2 之外,p-1 必然是偶数,那么必然存在公因数 2,回到上一个问题,至少在 1000 以内是不可能。在除以 2 之后,一定比 小。

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

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

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

常规直观做法

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
-- 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)

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

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

根据求和公式,可以得到

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

代码如下,需要避免整数溢出:

1
2
3
4
5
6
7
8
9
-- 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 的正整数中和 n 互质的数,而且满足

其中 p 取 n 的所有质因子。比如说 ,互质的数为,恰好

欧拉定理证明

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

欧拉函数的证明,分成各种情况:

素数幂

显然只有 与 n 有除 1 以外的公因数,k 取 1 到 个数,所以

互质的整数的乘积

如果 是互质的整数,则 。这是因为由于 互质,小于 的整数中与 互质的整数,可以通过取小于 的与 互质的整数,与小于 的与 互质的整数,的所有可能组合来构造。

我们需要考虑由 形成的整数对格点,并分析这些点中有多少是与 互质的。具体步骤如下:

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

因此,与 互质的整数对 的总数就是 的乘积,即 。综上所述,我们得出

一般情况

对于任意整数 n,存在唯一质数分解 。由于这些素数幂相互之间是互质的,根据第 2 步,我们可以写出

根据第 1 步,我们知道每个 。将它们代入上面的公式,我们得到

这可以进一步简化为

循环群

循环群是指能由单个元素所生成的群。设为一个群,若存在一个元素 ,使得 ,那么 形成循环群。群 内任意一个元素所生成的群都是循环群,而且是 的子群。

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