Qiita Teams that are logged in
You are not logged in to any team

Community
Service
Qiita JobsQiita ZineQiita Blog
1
Help us understand the problem. What is going on with this article?
@bra_cat_ket

# 99 Haskell Problems より、[31..41]

More than 5 years have passed since last update.

H-99: Ninety-Nine Haskell Problemsより [31..41]、以下が今までの軌跡です。

# 31. 素数判定。

こちらで議論したのはサイドストーリーになってしまいました、当初私のロジックとしては（何故か）素数列を作ってしまえばその数がリストに入ってるかどうか自明じゃん、と思ったところからずれていました。

ここで味噌はall を使うことによって、素数が”必ず入っている”candidates をsqrt(n) まで取った数で割ってみた結果が一個でも割りきれてしまえば((/= 0) . mod n) == False その瞬間にFalse を返すところですね。

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
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. 最大公約数をユークリッドの互除法を使って求めよ。

ここではオリジナルとされている引き算タイプを選択。

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 と互いに素な数字の数を返すオイラーの関数を計算せよ。

Euler's totient function

これまた前問の結果を使って互いに素かどうか判定したリストの長さを数えたら良いだけになります。

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 でひっくり返したら良いです。

(追記ここから)

(追記ここまで)

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 のアイディアが浮かんだのでちょっと実験してすぐ実装出来ました。

(追記ここから)

(追記ここまで)

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. 効率のよいオイラーの関数を実装せよ。

ここでいう効率のよい、とは問題文中に与えられている”公式”で与えられるものです。

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. ゴールドバッハ予想。

アイディアとしては与えられた偶数より小さいところまで素数列を取ってきて、その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] について返すような関数を書け。

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)
``````

1. もちろん形式証明のそれとは関係ありません、ここでは単に”充分”大きな数まで実験したということを計算機的と言ってみたまでです。また反例がこのしらみつぶし的な方法で見つかっても、これはゴールドバッハ予想に対する否定的な証明になるのであながち間違いでは無いと思います。

1
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
1. We will deliver articles that match you
By following users and tags, you can catch up information on technical fields that you are interested in as a whole
2. you can read useful information later efficiently
By "stocking" the articles you like, you can search right away
A Haskell, Maxima lover, PhD in physics. About to leave qiita, perhaps.