Newton法で平方根をもとめる

cf. http://hpcgi2.nifty.com/1to100pen/wiki/wiki.cgi?p=%CB%E8%C6%FCHaskell 2006-05-21

最近流行りの SICP の 1.1.7 Newton法からお題を拝借。

おわっと,流行だったのか。もしかして乗り遅れた?

上の [1..100]>>=pen さんのもそうだし,hasko さんとこ(id:hasko:20060520:1148133213)のコメントにもあるように,先に予測値のリストを作ってしまうことにしよう。

guess x = iterate (\g -> (g + x/g) / 2.0) 1.0

goodEnough eps old new = abs (1.0 - old/new) < eps

find2 f i (x:xs) | f i x = x
                 | otherwise = find2 f x xs

squareRoot x = find2 (goodEnough 0.0001) x (guess x) 

guess で予測値のリストを作って,goodEnough で判定。判定は新旧予測値の比を利用し,引数で指定できるようにした。 で,あとは適当な値をリストから探すだけなんだけど,ちょうどいい関数が見あたらなかったので find2 を定義してみた。 結果。

*Main> squareRoot 4
2.000000000000002
*Main> squareRoot 1.0e-6
1.0000000000000117e-3
*Main> squareRoot 1.0e29
3.162277660171076e14
*Main> (squareRoot 1.0e-6) ^ 2
1.0000000000000235e-6
*Main> (squareRoot 1.0e29) ^ 2
1.0000000000017057e29

ふむ,よさげ。

Haskellプログラマ進化論

えーと,どこで見つけたかわかんなくなっちゃったけど。

 The Evolution of a Haskell Programmer

階乗を求める fac 関数にこんなに書き方があるなんて。

俺はさしずめ Another senior Haskell programmer あたりかな。まだまだ先は遠いな。Interpretive Haskell programmer とか Origamist Haskell programmer とかなんて,何をやってるのかさっぱりだ。
でも最後の Tenured professor は可笑しい。

練習問題

入門Haskell―はじめて学ぶ関数型言語」 p.86 より。昨日より1ページ戻ったけど。

(1) repeat,cycle,iterate をそれぞれ自分で定義しなさい。

これは簡単。再帰を使ってこうだ。

myRepeat a = a:map id (myRepeat a)

myCycle a = a ++ myCycle a

myIterate f a = a:map f (myIterate f a)

結果は省略。

(2) repeat,cycle,iterate のうち,どれか1つを自分で定義し,残りの2つはその定義した関数をもとに定義しなさい。ほかの Prelude 関数を使ってもかまいません。

こっちはちょっと難しい。とくに iterate でかなりつまった。
まずは repeat を使って他の2つを定義。

myCycle2 = (foldr (++) []) . myRepeat

myIterate2 f a = map (\x -> x a) (rep f)
  where rep = (scanl (.) id) . myRepeat

結果。

*Main> take 50 $ myCycle "abc"
"abcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcab"
*Main> take 50 $ myIterate2 (+1) 1
[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,3
0,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50]

よし,うまくいってるみたいだ。myCycle2 はこないだ覚えた foldr で無限リスト。
myIterate2 のほうは,myRepeat で関数 f のリストを作っておいて,scanl1 で合成している。これで [id, f, (f . f), (f . f . f), …] というふうなリストができる。かなりトリッキーな気がするけど,いいんだろうか。ま,結果はうまくいってるみたいだけど。

myCycle を使って他の2つを定義。

myRepeat3 a = myCycle [a]

myIterate3 f a = map (\x -> x a) (cyc f)
  where cyc g = scanl (.) id (myCycle [g])

myRepeat3 は大した工夫ではなし,それができればあとは myIterate3 は上の myIterate2 とおなじ。

*Main> take 50 $ myRepeat3 'a'
"aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
*Main> take 50 $ myIterate3 (+1) 1
[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,3
0,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50]

最後。myIterate を使って。

myRepeat4 = myIterate id

myCycle4 = (foldr (++) []) . (myIterate id)
*Main> take 50 $ myRepeat4 'a'
"aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
*Main> take 50 $ myCycle4 "abc"
"abcabcabcabcabcabcabcabcabcabcabcabcabcabcabcabcab"

追記:
myCycle2 と myCycle4 は foldr じゃなくて concat を使った方が簡単だな。

myCycle2 = concat . myRepeat

myCycle4 = concat . (myIterate id)

原始的なピタゴラス数(こんどこそ)

cf. id:takatoh:20060519:pythagorean

昨日のダメダメなのを修正するぞ。
いや,実はアップしてからすぐに気が付いたんだけど時間がないから後にしよう,なんて考えたら,案の定ツッコまれた。
というかアップする前に気づけ。

というわけで,昨日のでは条件が足りなかったようだ。必要な条件は次の3つ。1つ目だけは昨日のにも入っていた。

  • m > n
  • m と n は互いに素
  • m – n が奇数

で,こうなった。

pPythagorean a = [(m*m-n*n, 2*m*n, z) | m <- [1..a], n <- [1..a], m>n, gcd m n == 1, odd (m-n), let z = m*m+n*n, z <= a]

結果。

*Main> pPythagorean 100
[(3,4,5),(5,12,13),(15,8,17),(7,24,25),(21,20,29),(9,40,41),(35,12,37),(11,60,61
),(45,28,53),(33,56,65),(13,84,85),(63,16,65),(55,48,73),(39,80,89),(77,36,85),(
65,72,97)]

これでOKかな。

練習問題

入門Haskell―はじめて学ぶ関数型言語」 p.87 より。

lattice 関数を改良し,現在の2次元格子点限定から任意のn次元格子点に対応するように Integer -> Integer -> [ [Integer] ] の型 を持つようにしなさい。また,それを利用して all_lattice も n次元格子点をすべて順番に列挙するように改良しなさい。

※注意:ヒントにtypo というか文字の取り違えがある。紛らわしい。

これは簡単。

lattice 1 n = [[n]]
lattice d n = [a:b | a <- [0..n], b <- lattice (d-1) (n-a)]

all_lattice d = concat [lattice d i | i <- [0..]]

結果。

*Main> take 50 $ all_lattice 2
[[0,0],[0,1],[1,0],[0,2],[1,1],[2,0],[0,3],[1,2],[2,1],[3,0],[0,4],[1,3],[2,2],[
3,1],[4,0],[0,5],[1,4],[2,3],[3,2],[4,1],[5,0],[0,6],[1,5],[2,4],[3,3],[4,2],[5,
1],[6,0],[0,7],[1,6],[2,5],[3,4],[4,3],[5,2],[6,1],[7,0],[0,8],[1,7],[2,6],[3,5]
,[4,4],[5,3],[6,2],[7,1],[8,0],[0,9],[1,8],[2,7],[3,6],[4,5]]
*Main> take 50 $ all_lattice 3
[[0,0,0],[0,0,1],[0,1,0],[1,0,0],[0,0,2],[0,1,1],[0,2,0],[1,0,1],[1,1,0],[2,0,0]
,[0,0,3],[0,1,2],[0,2,1],[0,3,0],[1,0,2],[1,1,1],[1,2,0],[2,0,1],[2,1,0],[3,0,0]
,[0,0,4],[0,1,3],[0,2,2],[0,3,1],[0,4,0],[1,0,3],[1,1,2],[1,2,1],[1,3,0],[2,0,2]
,[2,1,1],[2,2,0],[3,0,1],[3,1,0],[4,0,0],[0,0,5],[0,1,4],[0,2,3],[0,3,2],[0,4,1]
,[0,5,0],[1,0,4],[1,1,3],[1,2,2],[1,3,1],[1,4,0],[2,0,3],[2,1,2],[2,2,1],[2,3,0]
]

OK。…と思ったら型が合わない。

*Main> :t lattice
lattice :: (Num a, Num a1, Enum a) => a1 -> a -> [[a]]

ああ,そうか。Int でもかまわないものな。それに,引数に整数以外を与えてしまうとおかしなことになる。
これは宣言してやればいいだけだ。

lattice :: Integer -> Integer -> [[Integer]]

リスト内包表記と無限リスト

(i,j) の無限リストを得ようとして,素朴にこんな風にやっても期待通りにはいかない。

Prelude> [(i,j) | i <- [0..], j <- [0..]]
[(0,0),(0,1),(0,2),(0,3),(0,4),(0,5),(0,6),(0,7),(0,8),(0,9),(0,10),(0,11),(0,12
),(0,13),(0,14),(0,15),(0,16),(0,17),(0,18),(0,19),(0,20),(0,21),(0,22),(0,23),(
0,24),(0,25),(0,26),(0,27),(0,28),(0,29),(0,30),(0,31),(0,32),(0,33),(0,34),(0,3
5),(0,36),(0,37),(0,38),(0,39),(0,40),(0,41),(0,42),(0,43),(0,44),(0,45),(0,46),
(以下略)

いつまでたっても i が 2 にならない。i が 2 になる前に, i が 1 で j が 1.. の (i,j) をすべて列挙しようとするからだ。

こういう場合は,まず i+j == 1 であるような (i,j) を列挙し,次に i+j == 2 の (i,j) を……というふうにする。こんな感じ。

all_lattice = concat [lattice k | k <- [0..]] lattice k = [(i, k-i) | i <- [0..k]]
*Main> take 100 $ all_lattice
[(0,0),(0,1),(1,0),(0,2),(1,1),(2,0),(0,3),(1,2),(2,1),(3,0),(0,4),(1,3),(2,2),(
3,1),(4,0),(0,5),(1,4),(2,3),(3,2),(4,1),(5,0),(0,6),(1,5),(2,4),(3,3),(4,2),(5,
1),(6,0),(0,7),(1,6),(2,5),(3,4),(4,3),(5,2),(6,1),(7,0),(0,8),(1,7),(2,6),(3,5)
,(4,4),(5,3),(6,2),(7,1),(8,0),(0,9),(1,8),(2,7),(3,6),(4,5),(5,4),(6,3),(7,2),(
8,1),(9,0),(0,10),(1,9),(2,8),(3,7),(4,6),(5,5),(6,4),(7,3),(8,2),(9,1),(10,0),(
0,11),(1,10),(2,9),(3,8),(4,7),(5,6),(6,5),(7,4),(8,3),(9,2),(10,1),(11,0),(0,12
),(1,11),(2,10),(3,9),(4,8),(5,7),(6,6),(7,5),(8,4),(9,3),(10,2),(11,1),(12,0),(
0,13),(1,12),(2,11),(3,10),(4,9),(5,8),(6,7),(7,6),(8,5)]

原始的なピタゴラス数

互いに素である3数からなるピタゴラス数を原始的な(または素な)ピタゴラス数というらしい。

 cf. http://ja.wikipedia.org/wiki/%E3%83%94%E3%82%BF%E3%82%B4%E3%83%A9%E3%82%B9%E6%95%B0

これには公式があって,内包表記で次のように書ける。公式について詳しくは上のリンク先を見よ。
頭の p は primitive の p。それから let ~ も使ってみた。

pPythagorean a = [(m*m-n*n, 2*m*n, z) | m <- [1..a], n <- [1..a], m>n, let z = m*m+n*n, z <= a]
*Main> pPythagorean 100
[(3,4,5),(8,6,10),(5,12,13),(15,8,17),(12,16,20),(7,24,25),(24,10,26),(21,20,29)
,(16,30,34),(9,40,41),(35,12,37),(32,24,40),(27,36,45),(20,48,52),(11,60,61),(48
,14,50),(45,28,53),(40,42,58),(33,56,65),(24,70,74),(13,84,85),(63,16,65),(60,32
,68),(55,48,73),(48,64,80),(39,80,89),(28,96,100),(80,18,82),(77,36,85),(72,54,9
0),(65,72,97)]

z <= a を let z = ~ よりも前に書いてはいけない。Not in scope: `z’ だと怒られる。

リスト内包表記

リストを作るとき,要素を列挙したり範囲を指定するほかに,計算をしながら作ることもできる。
たとえば100以下の偶数のリストは

*Main> [x | x <- [1..100], x `mod` 2 == 0]
[2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56
,58,60,62,64,66,68,70,72,74,76,78,80,82,84,86,88,90,92,94,96,98,100]

| の左側がリストの要素,右側がその条件になる。条件には次のことが書ける。

  • パターン <- リスト
  • 真偽式
  • let パターン = 式

上の例では,パターン <- リスト と真偽式と使っている。3番目の let ~は「入門Haskell」には詳しい説明がないけど,「あまり使わない」とか書いてあるからそういうことにしとこう。

以下,練習。
三角数

triangular a = [n*(n+1)/2 | n <- [1..a]]
*Main> triangular 50
[1.0,3.0,6.0,10.0,15.0,21.0,28.0,36.0,45.0,55.0,66.0,78.0,91.0,105.0,120.0,136.0
,153.0,171.0,190.0,210.0,231.0,253.0,276.0,300.0,325.0,351.0,378.0,406.0,435.0,4
65.0,496.0,528.0,561.0,595.0,630.0,666.0,703.0,741.0,780.0,820.0,861.0,903.0,946
.0,990.0,1035.0,1081.0,1128.0,1176.0,1225.0,1275.0]

平方数

square a = [n*n | n <- [1..a]]
*Main> square 50
[1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529
,576,625,676,729,784,841,900,961,1024,1089,1156,1225,1296,1369,1444,1521,1600,16
81,1764,1849,1936,2025,2116,2209,2304,2401,2500]

ピタゴラス数

pythagorean a = [(x, y, z) | x <- [1..a], y <- [1..a], z <- [1..a], x < y, x*x + y*y == z*z]
*Main> pythagorean 50
[(3,4,5),(5,12,13),(6,8,10),(7,24,25),(8,15,17),(9,12,15),(9,40,41),(10,24,26),(
12,16,20),(12,35,37),(14,48,50),(15,20,25),(15,36,39),(16,30,34),(18,24,30),(20,
21,29),(21,28,35),(24,32,40),(27,36,45),(30,40,50)]

ピタゴラス数って結構あるんだな。

foldr と無限リスト

昨日の IO () さんのコメントから。

 cf. id:takatoh:20060514:myInits

えーと,つまり foldl は無限リストを扱えないけど, foldr は扱えるってことだな。foldl を使った版の myInits がダメだったのは last じゃなくて foldl のせいだってことだ。

それより,foldr は無限リストを扱えるって?

Prelude> foldr (:) [] [1..]
[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,3
0,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,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83
,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107

うわ,ホントだ。なんでだ。リストの右(後ろ)から評価するんじゃないのか。
うーん,”Yet Another Haskell Tutorial”を読んでみようか。英語は苦手なんだがなぁ……と思ったら id:moi2 さんのところに翻訳がある。ありがたい。該当個所はこのあたりか。

この中に foldr の評価過程の例がある。

foldr (-) 1 [4,8,5]
==> 4 – (foldr (-) 1 [8,5])
==> 4 – (8 – foldr (-) 1 [5])
==> 4 – (8 – (5 – foldr (-) 1 []))
==> 4 – (8 – (5 – 1))
==> 4 – (8 – 4)
==> 4 – 4
==> 0

上の例をこれと同じようにしてみると

foldr (:) [] [1..]
==> 1:(foldr (:) [] [2..])
==> 1:2:(foldr (:) [] [3..])
==> 1:2:3:(foldr (:) [] [4..])
==> 1:2:3:4:(foldr (:) [] [5..])
(以下無限に続く)

こんな感じになるはずで,無限リストを扱えるようには思えない。
いや,そうか,リストだから遅延評価されるんだ。だからうまくいくってことか。
ということは結果がリストにならないような場合にはうまくいかないはず。

Prelude> foldr (+) 0 [1..]
*** Exception: stack overflow

やっぱり。