H-99: Ninety-Nine Haskell Problemsより [31..41]、以下が今までの軌跡です。
99 Haskell Problems より、[1..10]
99 Haskell Problems より、[11..20]
99 Haskell Problems より、[21..28]
99 Haskell Problems より、[31] でもその前に
31. 素数判定。
こちらで議論したのはサイドストーリーになってしまいました、当初私のロジックとしては（何故か）素数列を作ってしまえばその数がリストに入ってるかどうか自明じゃん、と思ったところからずれていました。
明らかですが、素数判定に素数列は大雑把に言えば要りません、その数の平方根までで割れるかどうかを見れば良いので。
解答からアイディアをもらって書き直しましたがそれでも解答のパフォーマンスの方が良いですね。
ここで味噌はall を使うことによって、素数が”必ず入っている”candidates をsqrt(n) まで取った数で割ってみた結果が一個でも割りきれてしまえば((/= 0) . mod n) == False その瞬間にFalse を返すところですね。
all の評価方法は必ずしもHaskell だけのものというわけではありませんが、上手な使い方だと思います。
以下が実装になります、少々見にくいですが、当初の実装から型とかいわゆるコーナーケースを変えたりしたのでその痕跡が残っているだけです。
Prob31.lhs
> module Prob31 where
Determine whether a given integer number is prime.
see for the implementation:
http://qiita.com/bra_cat_ket/items/6a99a9b01682886607d0
> merge :: (Ord a) => [a] -> [a] -> [a]
> merge xx@(x:xs) yy@(y:ys)
> | x < y = x : (merge xs yy)
> | x == y = x : (merge xs ys)
> | x > y = y : (merge xx ys)
> diff :: (Ord a) => [a] -> [a] -> [a]
> diff xx@(x:xs) yy@(y:ys)
> | x < y = x : (diff xs yy)
> | x == y = diff xs ys
> | x > y = diff xx ys
> -- primes, nonPrimes :: [Integer]
> primes, nonPrimes :: Integral a => [a]
> primes = 2:3:5: (diff [7,9..] nonPrimes)
> nonPrimes = foldr1 merge' $ map helper $ tail primes
> where
> merge' (x:xs) ys = x : (merge xs ys)
> helper p = [n*p | n <- [p, p+2..]]
*Prob31> 7917 `elem` takeWhile (<7919) primes
False
> -- isPrime :: Integer -> Bool
> isPrime :: Integral a => a -> Bool
> -- isPrime n
> -- | odd n = n `elem` takeWhile (< n+2) primes
> -- | otherwise = False
> isPrime n
> | odd n = all ((/= 0) . mod n) $ takeWhile (<= upperLimit) primes
> | otherwise = False
> where upperLimit = floor . sqrt $ fromIntegral n
Of course, we don't have to make primes to check whether each input is prime or not!
So this is clearly a super inefficient implementation.
The below implementation is from
https://wiki.haskell.org/99_questions/Solutions/31
All primes larger than 2,3 can be written as
6*m + 1, 6*m + 5
forms.
In addition, all we have to do is to check whether the input n has its factor up to sqrt(n)
> isPrime' :: Integral a => a -> Bool
> isPrime' n
> | n < 2 = False
> | n == 2 && n == 3 = True
> -- isPrime' n
> -- | n < 4 = n /= 1
> isPrime' n = all ((/=0) . mod n) $ takeWhile (<= upperLimit) candidates
> where
> candidates = (2:3: [x+i | x <- [6,12..], i <- [-1,1]])
> upperLimit = floor . sqrt $ fromIntegral n
The following is, on the other hand, super slow.
> -- primes' :: [Integer]
> primes' :: Integral a => [a]
> primes' = 2 : [n| n<-[3,5..], isPrime' n]
*Prob31> isPrime 982451653
True
(0.19 secs, 38,213,232 bytes)
*Prob31> isPrime' 982451653
True
(0.03 secs, 6,751,536 bytes)
32. 最大公約数をユークリッドの互除法を使って求めよ。
例題にあった負の数の扱いに少し頭を捻りましたが、アルゴリズムの解答はwiki に載っていたので楽ちんでした。
ここではオリジナルとされている引き算タイプを選択。
Prob32.lhs
> module Prob32 where
Determine the greatest common divisor of two positive integer numbers. Use Euclid's algorithm.
https://en.wikipedia.org/wiki/Euclidean_algorithm
The following implementation is the subtraction-based version which was Euclid's original version.
> myGCD :: Integral a => a -> a -> a
> myGCD a b
> | b < 0 = myGCD a (-b)
> myGCD a b
> | a == b = a
> | b > a = myGCD b a
> | b < a = myGCD (a-b) b
*Prob32> [myGCD 36 63, myGCD (-3) (-6), myGCD (-3) 6]
[9,3,3]
*Prob32> myGCD 1071 462
21
33. 互いに素。
イヤに簡単な問題が続きますね、というか31. が一見さんには難しかったので比較の問題かもしれません。
32. で作った最大公約数を見てあげたら良いだけです。
Prob33.lhs
> module Prob33 where
Determine whether two positive integer numbers are coprime. Two numbers are coprime if their greatest common divisor equals 1.
> import Prob32 (myGCD)
> coprime :: Integral a => a -> a -> Bool
> coprime n m = (myGCD n m == 1)
*Prob33> coprime 35 64
True
34. 与えられた正の整数m に対して、[1..m-1] の中でm と互いに素な数字の数を返すオイラーの関数を計算せよ。
これまた前問の結果を使って互いに素かどうか判定したリストの長さを数えたら良いだけになります。
Prob34.lhs
> module Prob34 where
Calculate Euler's totient function phi(m).
Euler's so-called totient function phi(m) is defined as the number of positive integers r (1 <= r < m) that are coprime to m.
> import Prob33 (coprime)
*Prob33> length $ filter (coprime 10) [1..10]
4
> totient :: Integral a => a -> Int
> totient m
> | m > 0 = length $ filter (coprime m) [1..m]
*Prob34> map totient [1..20]
[1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,8]
https://oeis.org/A000010
35. 素因数分解。
お決まりですね、せっかく31. で効率の良い素数列の生成を勉強したので使わない手はありません。
ここでの実装は単純に素数列の先頭からタメシワリをしてみてそれを(:) で積んでいくというアルゴリズムを選択したので出力は例のそれとは逆順になります、これも必要ならばreverse でひっくり返したら良いです。
(追記ここから)
地下鉄で読み物してたらcons(:) を使う場所次第じゃんと思ったので試してみたらreverse いらない実装出来ました。
直観的にはこの２つの対応はfoldl とfoldr の対応と一緒なので、生の再起じゃなくてこれら高階関数使って書けるはずだと思うのですが、ちょっとうまく行かなかったですね、もう少し考えてみます。
(追記ここまで)
Prob35.lhs
> module Prob35 where
Determine the prime factors of a given positive integer. Construct a flat list containing the prime factors in ascending order.
> import Prob31 (primes)
*Prob35> takeWhile ((>) . floor . sqrt . fromIntegral $ 315) primes
[2,3,5,7,11,13]
> divides :: Integral a => a -> a -> Bool
> divides n m = (n `mod` m == 0)
>
> primeFactors :: Integral a => a -> [a]
> primeFactors n = helper n primes []
>
> helper :: Integral a => a -> [a] -> [a] -> [a]
> helper n pList@(p:ps) factors
> | n == 1 = factors
> | n `divides` p = helper n' pList (p:factors)
> | otherwise = helper n ps factors
> where n' = n `div` p
*Prob35> primeFactors 315
[7,5,3,3]
> primeFactors' :: Integral a => a -> [a]
> primeFactors' n = let
> helper' n pList@(p:ps) factors
> | n == 1 = factors
> | n `divides` p = p : (helper' n' pList factors)
> | otherwise = helper' n ps factors
> where n' = n `div` p
> in
> helper' n primes []
*Prob35> primeFactors' 315
[3,3,5,7]
36. 素因数分解、"指数"表記で。
今度はフォーマットを変えて、指数表記で出力させるやつですね。
とっさにData.List.group とzip のアイディアが浮かんだのでちょっと実験してすぐ実装出来ました。
(追記ここから)
前問35. の改変に伴って指数表記も昇順に出来ます。
(追記ここまで)
Prob36.lhs
> module Prob36 where
Determine the prime factors of a given positive integer.
Construct a list containing the prime factors and their multiplicity.
> import Prob35 (primeFactors, primeFactors')
> import Data.List (group)
*Prob36> Data.List.group $ primeFactors 315
[[7],[5],[3,3]]
*Prob36> map length it
[1,1,2]
*Prob36> map head $ Data.List.group $ primeFactors 315
[7,5,3]
*Prob36> map length $ Data.List.group $ primeFactors 315
[1,1,2]
*Prob36> zip [7,5,3] it
[(7,1),(5,1),(3,2)]
> primeFactorsMult :: Integral a => a -> [(a, Int)]
> primeFactorsMult n = zip factors' powers
> where
> factors = primeFactors n
> factors' = map head $ group $ factors
> powers = map length $ group $ factors
> primeFactorsMult' :: Integral a => a -> [(a, Int)]
> primeFactorsMult' n = zip factors' powers
> where
> factors = primeFactors' n
> factors' = map head $ group $ factors
> powers = map length $ group $ factors
*Prob36> primeFactorsMult' 315
[(3,2),(5,1),(7,1)]
37. 効率のよいオイラーの関数を実装せよ。
ここでいう効率のよい、とは問題文中に与えられている”公式”で与えられるものです。
本当は証明してから実装するべきなんでしょうけれど、、、
具体的な式表示があるとHaskell の強みの”そのまま”書ける感じでごり押せます。
解答例のリスト内包表記は思いつきたかった、、、
Prob37.lhs
> module Prob37 where
> import Prob34 (totient)
> import Prob36 (primeFactorsMult)
Calculate Euler's totient function phi(m) (improved).
See problem 34 for the definition of Euler's totient function.
If the list of the prime factors of a number m is known in the form of problem 36 then the function phi(m) can be efficiently calculated as follows:
Let
[(p1, m1), (p2, m2), (p3, m3)..]
be the list of prime factors (and their multiplicities) of a given number m.
Then phi(m) can be calculated with the following formula:
(p1-1)*p1^(m1-1) * (p2-1)*p2^(m2-1) * ..
*Prob37> primeFactorsMult 315
[(7,1),(5,1),(3,2)]
*Prob37> map (\(p,m) -> (p-1)*p^(m-1)) it
[6,4,6]
*Prob37> product it
144
> phi :: Integral a => a -> a
> phi n = product $ map (\(p,m) -> (p-1)*p^(m-1)) $ primeFactorsMult n
*Prob31> isPrime 12073
True
*Prob37> phi 12073
12072
(0.06 secs, 19,111,352 bytes)
*Prob37> totient 12073
12072
(2.18 secs, 604,187,792 bytes)
Using list comprehension, here is a beautiful solution:
> phi' :: Integral a => a -> a
> phi' n = product [(p-1)*p^(m-1) | (p,m) <- primeFactorsMult n]
*Prob37> totient 10090
4032
(1.59 secs, 438,197,056 bytes)
*Prob37> phi 10090
4032
(0.01 secs, 4,160,272 bytes)
*Prob37> phi' 10090
4032
(0.01 secs, 2,613,072 bytes)
38. 以上２つのオイラー関数を比較せよ。
前問の中でやりました。
リスト内包表記のが一番良さそう。
Compare the two methods of calculating Euler's totient function.
Use the solutions of problems 34 and 37 to compare the algorithms. Take the number of reductions as a measure for efficiency. Try to calculate phi(10090) as an example.
39. 指定する範囲内の素数のリスト。
31.で素数列の生成を頑張ったかいがここでありました！
というわけでさくっと具体的に実験してみました、Prelude の関数だけで切り取れます。
Prob39.lhs
> module Prob39 where
A list of prime numbers.
Given a range of integers by its lower and upper limit, construct a list of all prime numbers in that range.
> import Prob31 (primes)
*Prob39> dropWhile (<10) $ takeWhile (<20) primes
[11,13,17,19]
(0.01 secs, 3,617,392 bytes)
> primesR :: Integral a => a -> a -> [a]
> primesR lowerL upperL = dropWhile (< lowerL) $ takeWhile (< upperL) primes
40. ゴールドバッハ予想。
任意の偶数は２つの素数和で表せられるという有名な予想です、計算機的な”証明”とも言えるでしょうか？1
アイディアとしては与えられた偶数より小さいところまで素数列を取ってきて、そのlast を取ってそれに再帰的に前の数を足してみてチェックしていくという方法を取りました。
解答からはこれまた大変美しいリスト内包表記。
アイディアは一緒ですが書き方を知ってるとこんなにすっきりコンパクトに書けます。
しかもパフォーマンスは遜色ない、素晴らしい。
Prob40.lhs
> module Prob40 where
Goldbach's conjecture.
Goldbach's conjecture says that every positive even number greater than 2 is the sum of two prime numbers.
Example: 28 = 5 + 23.
It is one of the most famous conjecture in number theory that has not been proved to be correct in the general case.
It has been numerically confirmed up to very large numbers (much larger than we can go with our Prolog system).
Write a predicate to find the two prime numbers that sum up to a given even integer.
*Prob40> takeWhile (<28) primes
[2,3,5,7,11,13,17,19,23]
*Prob40> init $ takeWhile (<28) primes
[2,3,5,7,11,13,17,19]
*Prob40> last $ takeWhile (<28) primes
23
> import Prob31
> import Prob39 (primesR)
> goldbach :: Integral a => a -> (a,a)
> goldbach n
> | odd n = error "input should be an even number"
> | even n = helper pList n
> where
> pList = takeWhile (<n) primes
> helper pp@(p:ps) n
> | biggest + p == n = (p, biggest)
> | biggest + p > n = helper newPrimeList n
> | biggest + p < n = helper ps n
> where biggest = last pp
> newPrimeList = init pp
> goldbachL :: Integral t => t -> [(t, t)]
> goldbachL n = [(x,y) | x<-pr, y<-pr, x<y, x+y == n]
> where pr = primesR 2 (n-2)
> goldbach' :: Integral a => a -> (a,a)
> goldbach' = head . goldbachL
*Prob40> goldbach' 1000000
(17,999983)
(15.05 secs, 5,005,535,400 bytes)
*Prob40> goldbach 1000000
(17,999983)
(14.35 secs, 4,895,268,680 bytes)
41. 与えられた範囲の偶数についてゴールドバッハ分解を行え。また分解の小さい方が例えば50より大きいような場合を[2..3000] について返すような関数を書け。
二つ目はまぁまぁ面白そうですが、ひとつ目は簡単ですね。
今までの流れでリスト内包表記だろうと思ったのでさっくり書けます。
二つ目はfilter を使うことを思いつけば、あとは実際に使うフィルターの実装だけでこの場合は単に射影してそれが例えば50より大きいかどうかを判定したら良いのでラムダ（無名）関数で書けます。
Prob41.lhs
> module Prob41 where
Given a range of integers by its lower and upper limit, print a list of all even numbers and their Goldbach composition.
In most cases, if an even number is written as the sum of two prime numbers, one of them is very small.
Very rarely, the primes are both bigger than say 50.
Try to find out how many such cases there are in the range 2..3000.
> import Prob40 (goldbach)
> goldbachList :: Integral a => a -> a -> [(a,(a,a))]
> goldbachList lowerL upperL
> | lowerL > upperL = error "wrong input"
> | otherwise = [(x, goldbach x)| x <- [lowerL .. upperL], even x]
> goldbachList' :: Integral t => t -> t -> t -> [(t, (t, t))]
> goldbachList' lowerL upperL threshold
> = filter (\(_,(b,_)) -> b > threshold) $ goldbachList lowerL upperL
*Prob41> goldbachList 9 20
[(10,(3,7)),(12,(5,7)),(14,(3,11)),(16,(3,13)),(18,(5,13)),(20,(3,17))]
(0.01 secs, 2,575,184 bytes)
*Prob41> goldbachList' 3 2000 50
[(992,(73,919)),(1382,(61,1321)),(1856,(67,1789)),(1928,(61,1867))]
(2.27 secs, 793,645,824 bytes)
-
もちろん形式証明のそれとは関係ありません、ここでは単に”充分”大きな数まで実験したということを計算機的と言ってみたまでです。また反例がこのしらみつぶし的な方法で見つかっても、これはゴールドバッハ予想に対する否定的な証明になるのであながち間違いでは無いと思います。 ↩