DouKaku?beta

challenge 格子点の列挙

二次元平面上の格子点(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 がなかったので。力技。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
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 - Flatten

Nested Hidden 
1000番目は[1250,-1250]かな？

[[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1], [1, 1], [-1, 1], [-1, -1], ...]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
def lattice(limit)
  ary = [[0,0]]
  i = 1
  n = 1
  while i < limit
    ary << [n,0] << [0,n] << [-n,0] << [0,-n]
    ary << [n,n] << [-n,n] << [-n,-n] << [n,-n]
    i+=8
    n+=1
  end
  ary
end


ええと。
[1, 2]が[2, 2]より近い格子点だと言うことを忘れていませんか？

1000番目は(-8, 16)？Tkinterで視覚的に表示
するコードもつけました。debug()という関数を
実行してください。
 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
30
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
import itertools

def distance2(a):
    x, y = a
    return x * x + y * y

def rotate_90(x, y):
    return -y, x

def grid_points():
    yield (0, 0)
    def grid_points_in_first_quadrant():
        table = [] # index: y, value: x ( y >= 0 and x >= 1 )
        while 1:
            def candidates():
                for y, x in enumerate(table):
                    yield x + 1, y
                yield 1, len(table)
            x, y = min(candidates(), key=distance2)
            yield x, y
            if y < len(table):
                table[y] = x
            else:
                assert y == len(table)
                table.append(x)
    for _, it in itertools.groupby(grid_points_in_first_quadrant(), key=distance2):
        points = list(it)
        for _ in xrange(4):
            for x, y in points:
                yield x, y
            for i in xrange(len(points)):
                points[i] = rotate_90(*points[i])

def debug(): # show points visually
    import Tkinter as Tk
    import math

    root = Tk.Tk()
    canvas = Tk.Canvas(root, width=400, height=400)
    canvas.pack(fill=Tk.BOTH, expand=True)

    def draw_point(x, y):
        y = -y
        n = 10
        m = 5
        offset = 200
        canvas.create_oval(x * n - m + offset, y * n - m + offset, x * n + m + offset, y * n + m + offset, fill="red")

    it = grid_points()
    def command():
        x, y = it.next()
        draw_point(x, y)
        print (x, y), math.hypot(x, y)

    button = Tk.Button(text="Next", command=command)
    button.pack(fill=Tk.X)

    root.mainloop()

def main():
    for x, y in itertools.islice(grid_points(), 1000):
        print "%d, %d" % (x, y)
#   debug()

if __name__ == '__main__':
    main()

1000番目は (-8, 16) だと思います。 配列 x, y のダミーはかなり変ですが、もう直す気がしないのでそのままにしておきます。 
 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
30
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
#!/bin/bash

set -ue

FMT="%d, %d\n"

print_lattice () {
  # distance^2
  local -i d
  # cache of square of integers
  local -ia scache=(0)
  # #scache
  local -i n_scache="${#scache[@]}"

  local -i i j
  local -ia x
  local -ia y

  for ((d = 0;; d++)); do
    # set dummy data(due to "set -u")
    x=(-1)
    y=(-1)

    if [ $d -gt ${scache[n_scache-1]} ]; then
      scache[n_scache++]=$(( n_scache * n_scache ))
    fi

    for ((i = n_scache-1; i >= 0; i--)); do
      for ((j = 0; j < n_scache; j++)); do
        if let "scache[i]+scache[j]-d"; then
          continue
        fi

        x=("${x[@]}" $i)
        y=("${y[@]}" $j)
      done
    done

    if [ ${#x[@]} -gt 1 ]; then
      # origin
      if ! let "x[1] + y[1]"; then
        printf "$FMT" 0 0
      fi
      # quadrant 1
      for ((i = 1; i < ${#x[@]}; i++)); do
        if [ ${x[i]} -ne 0 ]; then
          printf "$FMT" ${x[i]} ${y[i]}
        fi
      done
      # quadrant 2
      for ((i = ${#x[@]} -1 ; i > 0; i--)); do
        if [ ${y[i]} -ne 0 ]; then
          printf "$FMT" -${x[i]} ${y[i]}
        fi
      done
      # quadrant 3
      for ((i = 1; i < ${#x[@]}; i++)); do
        if [ ${x[i]} -ne 0 ]; then
          printf "$FMT" -${x[i]} -${y[i]}
        fi
      done
      # quadrant 4
      for ((i = ${#x[@]} -1 ; i > 0; i--)); do
        if [ ${y[i]} -ne 0 ]; then
          printf "$FMT" ${x[i]} -${y[i]}
        fi
      done
    fi
  done
}

print_lattice


無駄が多いですが、普通に解いてみました。
ところで、反時計回りだとサンプルの出力は以下のようになるべきでは?  

0, 0
1, 0
0, 1
-1, 0
0, -1
1, 1
-1, 1
-1, -1
1, -1
2, 0
※(1, -1)と(-1, -1)が逆のはず。

ちなみに1,000番目は -8, 16 です。
 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
30
(use srfi-1)
(use srfi-42)
(use math.const)

(define (make-point-list n)
  (let ((a (ceiling->exact (/ (sqrt (* 2 n)) 2))))
    (list-ec (: x (- a) (+ a 1)) (: y (- a) (+ a 1))
             (make-rectangular x y))))

(define (sort-point-list lst)
  (define (%angle z)
    (let1 t (angle z)
      (if (< t 0)
          (+ (* 2 pi) t)
          t)))
  (sort lst (lambda (z1 z2)
              (let ((r1 (magnitude z1))
                    (r2 (magnitude z2)))
                (if (= r1 r2)
                    (< (%angle z1) (%angle z2))
                    (< r1 r2))))))

(define (lattice-point-list n)
  (for-each (lambda (z)
              (format #t "~d, ~d~%"
                      (inexact->exact (real-part z))
                      (inexact->exact (imag-part z))))
            (take (sort-point-list (make-point-list n)) n)))

;; (lattice-point-list 1000)を実行すると答えが出ます


xとyのループ範囲の計算方法を#3258から拝借しました(n2の値)。
ただし、なぜこうなるのかは理解できていません。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from math import sqrt

def f0(x, y):
  return (0 if x >= 0 else 1) if y >= 0 else (2 if x < 0 else 3)

def f(a, b):
  return -1 if a[0] < b[0] else 1 if a[0] > b[0] else cmp(a[1], b[1]) if a[1] != b[1] else (-1 if a[2] > b[2] else 1) if a[1] <= 1 else (-1 if a[2] < b[2] else 1)

n = 1000
n2 = int(sqrt(n * 2) / 2)

l = [(sqrt(x**2 + y**2), f0(x, y), x, y) for x in range(-n2, n2+1) for y in range(-n2, n2+1)]

for t in sorted(l, cmp=f)[:n]:
  print t[2], t[3]

簡単に説明すると、xとyが[-sqrt(n)/2, sqrt(n)/2]の範囲の正方形を考えた場合、その中の格子点の数はn個以上あります。

なので、この正方形の外接円内に含まれる格子点の数もn個以上あるはずです。よってこの円内に含まれる格子点を列挙して、ソートして最初のn個を表示すれば問題は解けます。

ただ、円内に含まれる格子点の列挙は面倒なので、この円に外接する正方形を考え、そこに含まれる格子点を列挙することで問題を解いています(ただし無駄が多くなる)。その正方形は x, yが[-sqrt(2n)/2, sqrt(2n)/2]となるので、#3258のような計算式となります。 

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

 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
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>
-}

Squeak Smalltalk で。

1. 円形のブラシを描く機能（#dotOfSize:）があるので、これで描いた円形のドット数を数え（primCountBits）て、同心円内に列挙すべき格子点をすべて含む最少の矩形格子を探し、その x （もしくは y）の最大値（max）を求めます。

2. -max から max まで整数で、重複を許す二項の順列組み合わせを生成（asDigitsToPower: 2 do: [...]）し、それぞれの項目を x、y に持つ、ポイントオブジェクト（x@y）としてコレクション（points）に収めます。

3. pointsを、各要素の原点からの距離（r）、距離が同じなら x 軸となす角度（theta）の小さい順に並べ変えます。

4. 前から n 個をとり（first: n）、各ポイントの x、y を出力。終了後、ファイルを開いて表示（edit）し、最後の出力を結果として返して（^last）います。 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
| n count max out points last |
n := 1000.
count := [:m | (Form dotOfSize: (m + 1 * 2)) primCountBits - (m + 1 * 4) + 1].
max := (1 to: Float infinity) detect: [:r | (count value: r) >= n].
out := FileStream forceNewFileNamed: 'doukaku65.out.txt'.
points := OrderedCollection new: max * max.
(max negated to: max) asDigitsToPower: 2 do: [:xy | points add: xy first @ xy second].
points := points asArray sort: [:a :b | a r < b r or: [a r = b r and: [a theta < b theta]]].
points := points first: n. 
points do: [:pt | out nextPutAll: (last := '{1}, {2}' format: {pt x. pt y}); cr].
out edit.
^last   "=> '-8, 16' "


格子点のリストを適当に生成しておいて、
原点からの距離でソートすることで近い順に並べてます。
反時計回りにはなってません。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
(defun distance (point)
  (let ((x (first point)) (y (second point)))
    (+ (* x x) (* y y))))

(defun list-lattice-points (n)
    (loop with a = (ceiling (sqrt (/ n 2)))
      for i from (- a) to a
      nconc (loop for j from (- a) to a collect (list i j)) into points
      finally (format t "~V{~{~D, ~D~}~%~}"
                      n (sort points #'<= :key #'distance))))


前出のコードでは

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

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

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
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)]


距離の二乗を返してるから distance って関数名はちょっと嘘ですね。

詳細な説明ありがとうございます。

私も先の投稿後、自分なりに考えてみました。

最初は1000個のペアがあればいいのだから、
　n = 1000
　n2 = sqrt(n) / 2
でいいかと考えましたが、この場合の最大値であるsqrt(2 * n2 ** 2)より小さくなる
(x, 0)あるいは(0, y)の組み合わせが存在するので
　n2 = sqrt(2 * (sqrt(n) / 2) ** 2)
とすべきで、変形すれば
　n2 = sqrt(n / 2)
となって、同じ結論になりました。


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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
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)]


ユークリッド距離でなくても、距離の公理を満たせば distance と呼んでもいいんじゃないかしらん。

円周をたどって。1000個目は -8, 16。
 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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <vector>
#include <iostream>
using namespace std;

void lattice(int num) {
	vector<pair<int, int>> v;
	vector<pair<int, int>>::iterator ite;
	int n = 1, r = 0, x_start = 0;
	
	if (num <= 0) return;
	cout << "0, 0" << endl;

	while (n < num) {
		if (r == x_start * x_start) ++x_start;
		int x = x_start;
		int y = 0;
		int next_r = x * x;
		while (x > 0) {
			int d = x * x + y * y;
			if (d <= r) {
				++x; ++y;
			} else {
				if (d < next_r) {
					v.clear();
					v.push_back(pair<int, int>(x, y));
					next_r = d;
				} else if (d == next_r) {
					v.push_back(pair<int, int>(x, y));
				}
				--x;
			}
		}
		
		for (ite = v.begin(); ite != v.end() && n < num; ++ite, ++n)
			cout << ite->first << ", " << ite->second << endl;
		for (ite = v.begin(); ite != v.end() && n < num; ++ite, ++n)
			cout << -ite->second << ", " << ite->first << endl;
		for (ite = v.begin(); ite != v.end() && n < num; ++ite, ++n)
			cout << -ite->first << ", " << -ite->second << endl;
		for (ite = v.begin(); ite != v.end() && n < num; ++ite, ++n)
			cout << ite->second << ", " << -ite->first << endl;

		v.clear();
		r = next_r;
	}
}

int main(void) {
	lattice(1000);
	return 0;
}


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

いささか汚いのが不満ですが、投稿してしまいます。
 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
30
31
32
33
34
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)


fを作用させる以前に同距離にある点に関して、出力順がお題を満たしていない気がします。
((5,0)と(4,3) など)

(n,0) ～ (n,n)
を配列に追加したら距離でソートして
距離が n 以下の部分を取り出して全周に展開して表示
…を繰り返すことで n<0x8000 あたりまで表示出来るはず。
コードはとりあえず1000番目を越えたら止まる設定ですが…。
1000番目は(-8,16)でした。

最初配列の先頭から要素を取り出すのに array_shift を使ったら
配列の残りのキーが0からの連番に上書きされてしまう事に気がつかなくて
困ったことに。
逆順にソートして array_pop なら良かったのか…。

 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
<?php
$c=0;
printf("%4d: %4d,%4d\n",++$c,0,0);
$a=array();
for($i=1; /*$i<0x8000*/ $c<=1000;++$i)
{	for($j=0;$j<=$i;++$j)
		$a[$i*$i+$j*$j][$i]=$j;
	ksort($a);
	for($d=$i*$i;key($a)<=$d;)
	{	$b=current($a);
		unset($a[key($a)]);
		foreach($b as $x => $y)
		{	if($y && $x!=$y)
				$b[$y]=$x;
		}
		krsort($b);
		for($j=0;$j<4;++$j)
		{	foreach($b as $x => $y)
			{	if($j&1) { $t=$x; $x=$y; $y=$t; }
				if(($j+1)&2) $x=-$x;
				if($j&2) $y=-$y;
				printf("%4d: %4d,%4d\n",++$c,$x,$y);
			}
		}
	}
}
?>

求める個数を引数として与えると、格子点を行列を返す。最後の三行のコメントを外すと求めた格子点をグラフ上に表示する。 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
function r = latticepointsaroundorigin(num)
% r = latticepointssorted(n) retruns a list of lattice points,
% in the order of  distance from the origin and polar angle.
% Uncomment the last three lines to show the locations of the points.
% (ja.doukaku.org Q65)
rg = ceil(sqrt(num)/2);
[x y] = meshgrid(-rg:rg);
pts = [x(:) y(:) x(:).^2+y(:).^2 atan2(-y(:),-x(:))];
r = sortrows(pts, [3,4]);
r = r(1:num,1:2);
% figure
% axis([-rg rg -rg rg])
% text(r(:,1),r(:,2),num2cell(1:num))

範囲の計算が変なまま出してしまったので再投稿。

同じく最後三行のコメントを外すと求めた格子点図示する。等距離の場合の順序をx軸を起点とした偏角の順にソートするという条件も満たしている。1000個目は [-8, 16]。 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
function r = latticepointsaroundorigin(num)
% r = latticepointssorted(n) retruns a list of lattice points,
% in the order of  distance from the origin and polar angle.
% Uncomment the last three lines to show the locations of the points.
% (ja.doukaku.org Q65)
rg = ceil(sqrt(num/pi));
[x y] = meshgrid(-rg:rg);
pts = [x(:) y(:) sqrt(x(:).^2+y(:).^2) atan2(-y(:),-x(:))];
r = sortrows(pts, [3,4]);
r = r(1:num,1:2);
% figure
% axis([-rg rg -rg rg])
% text(r(:,1),r(:,2),num2cell(1:num))


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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14