格子点の列挙 in Haskell DouKaku?

格子点の列挙

二次元平面上の格子点(X,Y座標がともに整数の点)を、原点から近い順に列挙してください。

同じ距離の点はどういう順番でも構いませんが、可能であればX軸に一番近い第一象限の点から原点を中心として反時計回りの順に列挙してください。列挙の方法は、1行に一つの点の、X,Y座標を出力することとします。

サンプル出力

0, 0
1, 0
0, 1
-1, 0
0, -1
1, 1
-1, 1
1, -1
-1, -1
2, 0

最低でも1000件まで列挙できることを確認してください。また「反時計回り」の条件も満たしている場合は、1000番目の頂点が何かも併せて答えてください。

このお題はかもさんの投稿を元にしています。ご協力ありがとうございました。

Perl がなかったので。力技。

use strict;

my $PI = atan2(1, 1) * 4;
my $MAXR = int(sqrt(1000 / $PI) + sqrt(2) + 0.5);
my @res = ();

for (my $i = 0; $i <= $MAXR; $i++) {
    for (my $j = 0; $j <= $MAXR; $j++) {
        my $r = sqrt($i * $i + $j * $j);
        push(@res, [$i, $j, $r, atan2($j, $i)]);
        ($i != 0) && push(@res, [-$i, $j, $r, atan2($j, -$i)]);
        ($j != 0) && push(@res, [$i, -$j, $r, atan2(-$j, $i) + 2 * $PI]);
        ($i * $j != 0) && push(@res, [-$i, -$j, $r, atan2(-$j, -$i) + 2 * $PI]);
    }
}

foreach my $p (sort {($a->[2] <=> $b->[2]) || ($a->[3] <=> $b->[3])} @res) {
    printf("%3d, %3d\n", splice(@{$p}, 0, 2));
}

post comment

Posted feedbacks - Haskell

効率はよくないけど、ナイーブな実装
1000番目(1-origin) は(-8,16)

grids :: [(Integer,Integer)]
grids = concat [ f (x,y) | z <- [0..], y <- [0..z], x <- [y..z], x^2+y^2==z ]
  where
    f (0,0)             = [(0,0)]
    f (x,0)             = [(x,0),(0,x),(-x,0),(0,-x)]
    f (x,y) | x == y    = [(x,y),(-x,y),(-x,-y),(x,-y)]
            | otherwise = [(x,y),(y,x),(-y,x),(-x,y),(-x,-y),(-y,-x),(y,-x),(x,-y)]

{-
*Main> mapM_ print $ take 1000 $ grids
(0,0)
(1,0)
(0,1)
(-1,0)
(0,-1)
(1,1)
(-1,1)
(-1,-1)
(1,-1)
(2,0)
... 中略 ...
(16,8)
(8,16)
(-8,16)               -- 1000番目(1-origin)
*Main>
-}

前出のコードでは

[ (x,y) | z <- [0..], y <- [0..z], x <- [y..z], x^2 + y^2 == z ]

で仰角０°から４５°の間で原点から近いほうから，かつ，仰角の小さい方か
ら格子点がならびます．上のコードでは y の上限 z になっていますが，これ
は√(z/2)を超えない自然数であればよいはずで，手抜きした分，余分な値を
チェックすることになってしまっています．

この部分を改良して，さらに monadic スタイルで書いたコードにしてみました．

import Control.Monad
grids :: [(Integer,Integer)]
grids = do { z <- [0..]
           ; let z' = floor (sqrt (fromInteger z / 2))
           ; y <- [0..z']
           ; x <- [y..z]
           ; guard (x^2 + y^2 == z)
           ; f (x,y)
           }
  where
    f (0,0)             = [(0,0)]
    f (x,0)             = [(x,0),(0,x),(-x,0),(0,-x)]
    f (x,y) | x == y    = [(x,y),(-x,y),(-x,-y),(x,-y)]
            | otherwise = [(x,y),(y,x),(-y,x),(-x,y),(-x,-y),(-y,-x),(y,-x),(x,-y)]

x の上限も√z でいいはずだ。
さらに効率アップ(って、最初からこう書くべきだった)

import Control.Monad
grids :: [(Integer,Integer)]
grids = do { z <- [0..]
           ; let z'  = floor (sqrt (fromInteger z / 2))
           ; let z'' = floor (sqrt (fromInteger z)) 
           ; y <- [0..z']
           ; x <- [y..z'']
           ; guard (x^2 + y^2 == z)
           ; f (x,y)
           }
  where
    f (0,0)             = [(0,0)]
    f (x,0)             = [(x,0),(0,x),(-x,0),(0,-x)]
    f (x,y) | x == y    = [(x,y),(-x,y),(-x,-y),(x,-y)]
            | otherwise = [(x,y),(y,x),(-y,x),(-x,y),(-x,-y),(-y,-x),(y,-x),(x,-y)]

0<=x<=yの格子点をn^2<x^2+y^2<=(n+1)^2の有限領域に切って、
そこをソートしていきます。
span2は有限領域を切り出すために定義しました。
(いい名前が思いつかなかった)

いささか汚いのが不満ですが、投稿してしまいます。

import Data.List (nub, groupBy, sortBy)

span2 :: (a -> Bool) -> ([a] -> Bool) -> [[a]] -> ([[a]],[[a]])
span2 p1 p2 zss = (xss, yss'++zss')
  where
    (yss,zss') = span p2 zss
    (xss,yss') = unzip $ map (span p1) yss

ans = iter 0 triangle
  where
    triangle = [[(x,y,x*x+y*y) | x <- [y..]] | y <- [0..]]

    iter n zss = xs ++ iter (n+1) zss'
      where
        (yss,zss') = span2 p1 p2 zss
        p1 (_,_,d)       = d <= n*n
        p2 ((_,y,_):xss) = y <= ceiling (fromInteger n / sqrt 2)
        xs = concat
           $ map mirror
           $ groupBy (\(_,_,d1) (_,_,d2) -> d1 == d2)
           $ sortBy (\(_,_,d1) (_,_,d2) -> compare d1 d2)
           $ concat yss
        -- sortByが安定であることを仮定

    mirror = nub
           . starling (++) (reverse.r3)
           . starling (++) (reverse.r2)
           . starling (++) (reverse.r1)
           . map (\(x,y,_)->(x,y))
      where
        r1 = map (\(x,y) -> (y,x))
        r2 = map (\(x,y) -> (-x,y))
        r3 = map (\(x,y) -> (x,-y))
        starling f g x = f x (g x)

ああ。ご指摘のとおりです。
結局仰角でsortすることにしてしまいました。

import Data.List
import Data.Complex

binapp o f x y = f x `o` f y
cmpapp :: Ord b => (a -> b) -> a -> a -> Ordering
cmpapp = binapp compare

grid :: Integer -> [(Integer,Integer)]
grid z = let z'  = floor $ sqrt (fromInteger z / 2)
             z'' = floor $ sqrt (fromInteger z)
         in  sortBy g $ concat [ f (x,y) | y <- [0..z'], x <- [y..z''], x^2+ y^2 == z ]
  where 
    f (0,0)             = [(0,0)] 
    f (x,0)             = [(x,0),(0,x),(-x,0),(0,-x)]
    f (x,y) | x == y    = [(x,y),(-x,y),(-x,-y),(x,-y)]
            | otherwise = [(x,y),(y,x),(-y,x),(-x,y),(-x,-y),(-y,-x),(y,-x),(x,-y)]
    g = cmpapp (h . phase . uncurry (:+) . e)
    h x = x + if x < 0 then 2*pi else 0
    e = uncurry (binapp (,) fromInteger)

grids :: [(Integer,Integer)]
grids = concatMap grid [0..]

{-
*Main> mapM_ print $ take 1000 $ grids
(0,0)
(1,0)
(0,1)
(-1,0)
(0,-1)
(1,1)
(-1,1)
(-1,-1)
(1,-1)
(2,0)
... 中略 ...
(16,8)
(8,16)
(-8,16)                      -- 1000番目(1-origin)
*Main>
-}

#3312 はログインせずに投稿していました。投稿し直し。
---

Haskell でラマヌジャン数を小さい順に列挙するときなどに
使われる手法を用いています。

import List (groupBy)

main = mapM_ print
  $ ((0,0):)
  $ concat
  $ map (rotate . map fst)
  $ groupBy (\(_,d1) (_,d2) -> d1 == d2)
  $ msort [[((x,y),x^2+y^2)| x <- [1..]] | y <- [0..]]

msort ((z:zs):zss) = z: fsort zs (msort zss)
fsort zs@(z@(_,d1):zs') ws@(w@(_,d2):ws')
  | d1 <= d2 = z: fsort zs' ws
  | d1 >  d2 = w: fsort zs  ws'

rotate zs = zs' ++ map (\(x,y) -> (-x,-y)) zs' where
  zs' = zs ++ map (\(x,y) -> (-y,x)) zs

三角波を使って生成。
1,000番目は 8, -14 でいいのかな。

triangle_wave :: Int -> Int -> Int
triangle_wave i n = abs(i `mod` (n * 4) - n * 2) - n

n_grids :: Int -> [(Int, Int)]
n_grids 0 = [(0, 0)]
n_grids n = map (\x -> (triangle_wave x n, triangle_wave (x - n) n)) [0..(n * 4 - 1)]

grids :: [(Int, Int)]
grids = concat $ map n_grids [0..]

main :: IO ()
main = mapM_ print $ take 1000 grids

遅延評価を使わないバージョン。

(1,0)-->(1,1)-->(1,2)-->(1,3)-->...
     +->(2,0)-->(2,1)->(2,2)->(2,3)->...
             +->(3,0)->(3,1)->(3,2)->(3,3)->...
                ...
という tree を走査する。

oceanさんの #3256 と多分本質部分は同じ。（Python 勉強してないので自信なし）
「議席数をドント方式で」の#1439（基はにしおさんの#1214）とも通じるところがある。

import Data.Graph.Inductive.Internal.Heap (empty, insert, splitMin)
import List (unfoldr, groupBy)

main = mapM_ print
  $ ((0,0):)
  $ concat
  $ map (rotate . map fst)
  $ groupBy (\(_,d1) (_,d2) -> d1 == d2)
  $ unfoldr f
  $ insertPs [(1,0)] empty

f heap = Just ((p,d), heap'') where
  ((d,_), p@(x,y), heap') = splitMin heap
  children = [(x,y+1)] ++ if (y > 0) then [] else [(x+1,0)]
  heap'' = insertPs children heap'

insertPs ps heap = foldr insert heap
  $ map (\(x,y) -> ((x^2+y^2,y), (x,y))) ps

rotate zs = let
  zs' = zs ++ map (\(x,y) -> (-y,x)) zs
  zs'' = zs' ++ map (\(x,y) -> (-x,-y)) zs'
  in zs''

インスタンス宣言のとこが汚いですが。第一象限で近い順に格子点を求め、それらを90度ずつ回転させてコピーします。

import Data.List
newtype Pt = Pt (Int, Int) deriving (Eq, Show)
instance Ord Pt where
    Pt (x1,y1) <  Pt (x2,y2)
        | x1^2 + y1^2 ==  x2^2 + y2^2 && y1==y2 = x1 > x2
        | x1^2 + y1^2 ==  x2^2 + y2^2 = y2 > y1
        | otherwise = x1^2 + y1^2 <  x2^2 + y2^2
    Pt (x1,y1) <= Pt (x2,y2) = x1^2 + y1^2 <= x2^2 + y2^2
mirror (Pt (x,y)) = [Pt (x,y), Pt (-y,x), Pt (-x,-y), Pt (y,-x)]
minimums l = partition (==m) l
    where m = minimum l
getMin plain = (concatMap mirror $ sort ms)
               ++ getMin (concatMap f ms ++ cand)
    where f (Pt (1,y)) = [Pt (2,y), Pt (1,y+1)]
          f (Pt (x,y)) = [Pt (x+1,y)]
          (ms, cand) = minimums plain
koushiten = Pt (0,0) : getMin [(Pt (1,0))]

DouKaku?

格子点の列挙

Posted feedbacks - Haskell

Index

Feed

Other

Link

Pathtraq