Hatena::Grouphaskell

Haskellで遊ぶよ

|

2009-06-24

「すべての奇数の合成数は、1つの整数の自乗の2倍と1つの素数の和で表される」?

10:51

Goldbach が考えた命題の反例を上げる問題。(一般に Goldbach 予想と呼ばれているものではない)


合成数というのは、9,15,21,25のように、2つ以上の素数の積となるもの。

例えば9は、

9 = 7 + 2 * 1^2

となり、「1つの整数の自乗の2倍と1つの素数の和」である。

この反例を見つける Haskell プログラム。

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

oddCompos :: [Integer]
oddCompos = filter isCompos [3,5..]
    where
        isCompos x = not.null $ [n| n <- [3,5..floor(sqrt(fromIntegral(x))) ], x `mod` n == 0]

main = print $ head $ filter isNotGoldbach oddCompos
    where
        isNotGoldbach x = all (/=x) [n*n*2+p | n <- takeWhile (\n->n*n*2<x) [1..], p <- takeWhile (\p->n*n*2+p<=x) primes]

結果。

5777

最後のリスト内包表記で、n*n*2 を3回もやっているのがなんとかならないかなー。リストの内包表記は、その中に出てくる値の状態を保持できないのがたまにキズ。なんと、let が使えたらしい。

また、一番最後の primes というところで、毎回最初から素数を求めてるんじゃないかと思ったけど、これは trace してみたところ、ちゃんと必要なところだけ計算しているみたいだった。(各 x に対して最初からやり直しているような書き方だけど、hugs でも最適化されている)

rst76rst762009/06/24 22:05関数 isNotGoldbach の中で n に対して p は一意に定まるので let で書けます。
これを直せばインタプリタでも十分速いはず。
n*n*2 の計算も一応まとめると、↓のようになるかと。
isNotGoldbach x = all (/=x) [n2+p | n2 <- takeWhile (<x) $ map (\n->n*n*2) [1..], let p = head $ dropWhile (<x-n2) primes]

edvakfedvakf2009/06/24 23:58そんなところで let が使えたのですね!! 便利です。
おっしゃる通り、かなり速くなりました。ありがとうございます。

トラックバック - http://haskell.g.hatena.ne.jp/edvakf/20090624

2009-06-11

takeWhile2とdropWhile2

00:46

takeWhile               :: (a -> Bool) -> [a] -> [a]
takeWhile p []          =  []
takeWhile p (x:xs) 
            | p x       =  x : takeWhile p xs
            | otherwise =  []

dropWhile               :: (a -> Bool) -> [a] -> [a]
dropWhile p []          =  []
dropWhile p xs@(x:xs')
            | p x       =  dropWhile p xs'
            | otherwise =  xs

これを参考に、リストのある要素と、その次の要素の2つを引数に取る関数を作る。

takeWhile2              :: (a -> a -> Bool) -> [a] -> [a]
takeWhile2 p []         =  []
takeWhile2 p (_:[])     =  []
takeWhile2 p (x:x':xs')
            | p x x'    =  x : (takeWhile2 p (x':xs'))
            | otherwise =  []

dropWhile2              :: (a -> a -> Bool) -> [a] -> [a]
dropWhile2 p []         =  []
dropWhile2 p (_:[])     =  []
dropWhile2 p xs@(x:x':xs')
            | p x x'    =  dropWhile2 p (x':xs')
            | otherwise =  xs

こういうことがやりたかったので。

Main> takeWhile2 (\x y -> y-x<50) $ dropWhile2 (\x y -> y-x<10) $ map (^2) [1..]
[25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,441,484,529,576,625]

学んだこと

xs@(x:xs')

全体としてパターン xs、部分としてパターン x:xs' にマッチする引数。


前の要素と比べる

ただ、takeWhile2でちょっと悩ましいのはtakeされるのが 評価関数の第一パラメタにわたる値なのか大にパラメタにわたる値なのか、どっちがいいのかは結構その場その場で変わってきそうなところのような気もしますね…

ずらしてみよう - 取り急ぎブログです

ごもっともです。

takeWhile2'             :: (a -> a -> Bool) -> [a] -> [a]
takeWhile2' p []        =  []
takeWhile2' p (_:[])    =  []
takeWhile2' p (x:x':xs')
            | p x x'    =  x : (takeWhile2' p (x':xs'))
            | otherwise =  [x]

dropWhile2'             :: (a -> a -> Bool) -> [a] -> [a]
dropWhile2' p []        =  []
dropWhile2' p (_:[])    =  []
dropWhile2' p xs@(x:x':xs')
            | p x x'    =  dropWhile2' p (x':xs')
            | otherwise =  x':xs'
Main> takeWhile2 (\x y -> y-x<20) $ dropWhile2 (\x y -> y-x<10) $ scanl1 (+) [1..]
[45,55,66,78,91,105,120,136,153,171]
Main> takeWhile2' (\x y -> y-x<20) $ dropWhile2' (\x y -> y-x<10) $ scanl1 (+) [1..]
[55,66,78,91,105,120,136,153,171,190]

(最初に書いてたやつはちょっと末端がおかしかったので密かに修正)

トラックバック - http://haskell.g.hatena.ne.jp/edvakf/20090611

2009-06-03

逆ポーランド記法電卓

07:18

苦労して作ったら、案の定ずっと良いのがあった。涙目。


自作版。訳あって Rational で結果を返したかったので。

import Char (isDigit,digitToInt,isSpace)
import Ratio
import System
import Debug.Trace

stringToRatio :: String -> Rational
stringToRatio xs = (foldl1 (\s x -> s*10 + x) $ map (toInteger.digitToInt) xs) % 1

isOp :: Char -> Bool
isOp x = x `elem` "+-*/"

toOp :: Char -> (Rational -> Rational -> Rational)
toOp '+' = (+)
toOp '-' = (-)
toOp '*' = (*)
toOp '/' = (/)

dropWhileSpace = dropWhile isSpace

parse :: ([Rational],String) -> ([Rational],String)
parse (mem, []) = (mem, [])
parse (mem, (x:xs))
    | (trace.show) (mem,x) False = (mem,x:xs)   -- for debugging
    | isOp x    = let ([arg2,arg1],mem') = splitAt 2 mem
                      op = toOp x
                   in parse ( (op arg1 arg2):mem', dropWhileSpace xs)
    | isDigit x = let (ns,xs') = span isDigit xs
                   in parse ( (stringToRatio (x:ns)):mem, dropWhileSpace xs')

poland :: String -> Rational
poland str = head $ fst $ parse ([], dropWhileSpace str)

main = do
        str  <-  getLine
        print $ poland str 

実行。

% runghc poland.hs
1 2 + 3 4 + *
([],'1')
([1%1],'2')
([2%1,1%1],'+')
([3%1],'3')
([3%1,3%1],'4')
([4%1,3%1,3%1],'+')
([7%1,3%1],'*')
21%1

Wikipedia 版

calc :: String -> [Float]
calc = foldl f [] . words
  where 
    f (x:y:zs) "+" = (y + x):zs
    f (x:y:zs) "-" = (y - x):zs
    f (x:y:zs) "*" = (y * x):zs
    f (x:y:zs) "/" = (y / x):zs
    f xs y = read y : xs

すごい簡潔。。words という素晴しい関数があったのか。read y : xs の部分も思いつかなかったなあ。。

って全部やん。


次の目標。禁則処理をちゃんとしたい。

自分のコードで parse の返り値を Maybe にすると、再帰部分まで書き換えるのが大変… ここは素直に Wikipedia 版を使おうかな。


とりあえずやってみた

import Char (isDigit)
import Ratio
import System

calc :: String -> Maybe [Rational]
calc = foldl poland (Just []) . words
  where
    poland :: Maybe [Rational] -> String -> Maybe [Rational]
    poland stack string = stack >>= flip f string -- (\stk -> f stk string) 

    f :: [Rational] -> String -> Maybe [Rational]
    f (x:y:zs) "+" = Just (y+x:zs)
    f (x:y:zs) "-" = Just (y-x:zs)
    f (x:y:zs) "*" = Just (y*x:zs)
    f (x:y:zs) "/" = Just (y/x:zs)
    f xs y         = if all isDigit y
                then Just (toRational (read y) : xs)
                else Nothing

main = do
        str  <-  getLine
        print $ calc str

Wikipedia 版で、禁則処理を入れるとこうなるのかなあ。気になる点が、

  1. f 関数の中で "Just" を使いすぎてて汚い。
  2. foldl に渡す初期値を Just [] としなければいけないのはダサい。

1番目は、f を純関数 (というのか?) にして禁則処理を別で定義すればいけるんだろうか?

2番目は難しそう。foldl の性質上、返る値が Maybe なら入れる値も Maybe でないといけない。

トラックバック - http://haskell.g.hatena.ne.jp/edvakf/20090603

2009-06-02

07:19

(++) :: [a] -> [a] -> [a]
[]     ++ ys = ys
(x:xs) ++ ys = x : (xs ++ ys)
  [1,2,3] ++ [4]
= 1 : ([2,3] ++ [4])
= 1 : (2 : ([3] ++ [4]))
= 1 : (2 : (3 ++ ([] ++ [4])))
= 1 : (2 : (3 : [4]))
= 1 : (2 : [3,4])
= 1 : [2,3,4]
= [1,2,3,4]

ステップ数は、(++) の第一引数の要素数 x 2 + 1

大きなリストの最後に何かくっつける場合は影響あるかも。


(!!)                :: [a] -> Int -> a
xs     !! n | n < 0 =  error "Prelude.!!: negative index"
[]     !! _         =  error "Prelude.!!: index too large"
(x:_)  !! 0         =  x
(_:xs) !! n         =  xs !! (n-1)
  [1,2,3,4] !! 2
= [2,3,4] !! 1
= [3,4] !! 0
= 3

ステップ数は、(!!) の第二引数 + 1


(++) も (!!) も O(N) の関数。ついでに length なんかも O(N)。リスト操作系はだいたいそう。ちょっと注意が必要。


探索

10:46

いいのみつけた。

探索

深さ優先探索
dfs :: (a -> [a]) -> a -> [a]
dfs f x = x:(f x >>= dfs f)
幅優先探索
bfs :: (a -> [a]) -> a -> [a]
bfs f = bfs' . (:[])
  where bfs' [] = []
        bfs' xs = xs ++ bfs' (xs >>= f)

-- 一行で書くと
bfs f = concat . takeWhile (not.null) . iterate (>> f) . (:[])
Programming:玉手箱:その他

こっちもすごく参考になる。

「ああいうことになってしまった原因は、与えられた問題に特化したコードを書こうとする姿勢にあると思われます。」

ハイすいません。あとで、そうなってしまった diff を書き直してみよう。。

トラックバック - http://haskell.g.hatena.ne.jp/edvakf/20090602

2009-06-01

diff

09:49

Haskell で diff を書こうと頑張ったんだけど、出入力とかのあたりに来たら急にヤル気が落ちていったので、とりあえずここらで書いておこうと思う。

アルゴリズムについては以下を読んで学んだ。リンク先の説明でもまだ分かりにくいのだけれど、自前では説明したくても面倒くさすぎる。

一つ目のリンクにはコード例が載っていて、(そのままでは動かないのだけど) 参考になった。

まず、以下が単純に O(NP) アルゴリズムによる探索をする部分。

onp :: (Eq a) => [a] -> [a] -> Int 
onp listA listB
    | delta < 0 = startONP listB listA
    | otherwise = head [p | p <- [0..], n == fp delta p ]
    where
        delta = n - (length listA)
        n = length listB
        fp = fpGen listA listB

-- generates fp 
-- fp returns y coordinate of the furthest point reachable on diagonal k, with detour p
fpGen :: (Eq a) => [a] -> [a] -> (Int -> Int -> Int)
fpGen listA listB = fp
    where 
        fp k p 
            | p == -1       = 0
            | k < - p       = 0
            | k > delta + p = 0
            | k == delta    = snake k $ max ((fp (k-1) p)     + 1) (fp (k+1) p)
            | k < delta     = snake k $ max ((fp (k-1) p)     + 1) (fp (k+1) (p-1))
            | k > delta     = snake k $ max ((fp (k-1) (p-1)) + 1) (fp (k+1) p)
        delta = (length listB) - (length listA)
        snake = snakeGen listA listB

-- generates snake
-- snake returns, for given k and y0, y coord of the point along the diagonal 
--      as far as the elems in A and B at the point are identical
snakeGen :: (Eq a) => [a] -> [a] -> (Int -> Int -> Int)
snakeGen listA listB = snake
    where
        snake k y0 = (+y0).length.takeWhile (==True) $ zipWith (==) (drop (y0-k) listA) (drop y0 listB)

main = print $ onp [1,2,3,4,5,6,7,8] [1,3,2,4,5,7,8,8,9,6]

与えられたリストによって、fp 関数 (例では配列だけど、Haskell では遅延評価があるので関数のほうが都合がいい)、と snake 関数を生成している。

上から4行目。

head [p | p <- [0..], n == fp delta p ]

普通なら fp(k,p) において、p を 0 から増やしつつ、k を [-p, delta + p] の範囲で変化させるロジックを書かなくてはならないのだけど、今回は「後ろから」評価することにしたため、

fp k p

とするだけでいい。探索の終了条件は、fp delta p が n と等しくなるとき。なので、head とリスト内包表記で書いた。

ここでちょっと問題があって、p が 0 から増やされる度に同じ引数で fp が計算されることがある。これはとりあえず課題として置いておく。

State モナドとかメモ化とか読んでみたけど、まだ全然理解できない。


次に、探索したルートを巡れるようにする。

そこで、k と y の値のタプルを繋げるリストを、Chain 型として定義した。

-- Chain type --
type Chain = [(Int,Int)]

stepRight,stepDown,slide :: Chain -> Chain
stepRight ((k,y):ns) = (k + 1 , y + 1):(k,y):ns
stepDown  ((k,y):ns) = (k - 1, y):(k,y):ns
slide     ((k,y):ns) = (k, y + 1):(k,y):ns

progress :: Chain -> Int
progress [] = 0
progress ((_,y):_) = y

diagonal :: Chain -> Int
diagonal [] = 0
diagonal ((k,_):_) = k

toXY :: Chain -> [(Int,Int)]
toXY [] = []
toXY ((k,y):ns) = (y-k,y):(toXY ns)
-- end --

onp :: (Eq a) => [a] -> [a] -> [(Int,Int)]
onp listA listB
    | delta < 0 = zipSwap $ onp listB listA
    | otherwise = reverse.toXY.head.dropWhile ((/=n).progress) $ map (fp delta) [0..]
    where
        delta = n - (length listA)
        n = length listB
        fp = fpGen listA listB

-- generates fp 
-- fp returns y coordinate of the furthest point reachable on diagonal k, with detour p
fpGen :: (Eq a) => [a] -> [a] -> (Int -> Int -> Chain)
fpGen listA listB = fp
    where 
        fp k p 
            | p < 0         = []
            | k < - p       = []
            | k > delta + p = []
            | k == delta    = snake k $ further (fp (k-1) p) (fp (k+1) p)
            | k < delta     = snake k $ further (fp (k-1) p) (fp (k+1) (p-1))
            | k > delta     = snake k $ further (fp (k-1) (p-1)) (fp (k+1) p)
        delta = (length listB) - (length listA)
        snake = snakeGen listA listB

-- analogy of "max (a+1) b" for Node
further :: Chain -> Chain -> Chain
further left upper
    | left == [] && upper == []          = [(0,0)]
    | left == []                         = stepDown upper
    | upper == []                        = stepRight left
    | progress left + 1 < progress upper = stepDown upper
    | otherwise                          = stepRight left

-- generates snake
-- snake returns, for given k and y0, y coord of the point along the diagonal 
--      as far as the elems in A and B at the point are identical
snakeGen :: (Eq a) => [a] -> [a] -> (Int -> Chain -> Chain)
snakeGen listA listB = snake
    where
        snake k chain = foldl (\chn _ -> slide chn) chain [1..len]
            where
                len = let y0 = progress chain
                      in  length.takeWhile (==True) $ zipWith (==) (drop (y0-k) listA) (drop y0 listB)

fp、snake、max という関数を、Chain 型を使って定義する。例えば fp は、(与えられた k と p に対して)「到達できる最遠点の y 座標」を返していたのだが、「到達できる最遠点までのルート」を返すようにした。

これで、onp 関数は、終点に到達するまでのルートの座標を返すようになった。


その座標から unified format に直す部分を考えていて、面倒になってきたので一旦ここで止めとく。

またヤル気が出たら始めるということで。

トラックバック - http://haskell.g.hatena.ne.jp/edvakf/20090601
|