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

Flatten Hidden 
ええと。
[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個のペアがあればいいのだから、
　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)
となって、同じ結論になりました。

ソート用比較関数を見直し。

1000番目の座標は他の人と同じく(-8, 16)になりました。

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

def f((x0, y0, z0), (x1, y1, z1)):
  if z0 < z1: return -1
  if z0 > z1: return 1
  if y0 >= 0:
    if y1 < 0 or x0 > x1: return -1
  else:
    if y1 < 0 and x0 < x1: return -1
  return 1

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

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

for x, y, z in sorted(l, cmp=f)[:n]:
  print x, y

効率はよくないけど、ナイーブな実装
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>
-}

前出のコードでは

[ (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)]

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

fを作用させる以前に同距離にある点に関して、出力順がお題を満たしていない気がします。
((5,0)と(4,3) など)
ああ。ご指摘のとおりです。
結局仰角でsortすることにしてしまいました。

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

距離の二乗を返してるから distance って関数名はちょっと嘘ですね。
ユークリッド距離でなくても、距離の公理を満たせば distance と呼んでもいいんじゃないかしらん。
それはそうなんですが、二乗すると三角不等式が成り立っていないという。
そうでした、(1/2)^2+(1/2)^2 < 1^2 ですもんね。
探索範囲を狭めました。
ついでに反時計回りになるようにしました。

評価の方法ですが、
円内の格子点に正方形を対応させると円の被覆ができるので、
これを使って格子点の数を円の面積で見積もるというのが大まかな方針です。

細かい計算は以下。
軸の周りで一対一対応にならなくてややこしいことになってます。

C: 原点を中心とする半径 r の閉円板
f: 単位正方形に対し、その頂点のうち最も原点に近いものを対応させる写像

C に属する格子点の数を K とすると、
f による像が C に入っているような単位正方形の数は
 原点に写るものが 4 個
 原点以外で軸上に写るものが 8r 個
 それ以外の点に写るものが K-4r-1 個
合わせて K+4r+3 個となります。
これらの単位正方形の和集合を S とします。

C は S の部分集合なので、S の面積は C の面積以上になります。
S は K+4r+3 個の単位正方形からなるので、その面積は K+4r+3 です。
また C の面積は pi*r^2 なので
pi*r^2 <= K + 4r + 3
が得られます。K が格子点の数だったので、格子点は少なくとも
pi*r^2 - 4r -3 個以上あることがわかります。

というわけで、格子点の数が n 以上となる半径 r を見つけるには
pi*r^2 - 4r - 3 >= n
を解けばよく、その結果
r >= (sqrt(4 + (3 + n) * pi) + 2) / pi
となります。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
(defun point< (p1 p2)
  (destructuring-bind ((x1 y1) (x2 y2)) `(,p1 ,p2)
     (let ((a1 (+ (* x1 x1) (* y1 y1)))
           (a2 (+ (* x2 x2) (* y2 y2))))
       (or (< a1 a2)
           (and (= a1 a2)
                (if (>= y1 0)
                    (or (> 0 y2) (> x1 x2))
                  (and (> 0 y2) (< x1 x2))))))))

(defun list-lattice-points (n)
  (loop with a = (ceiling (/ (+ (sqrt (+ 4 (* (+ 3 n) pi))) 2) pi))
    for i from (- a) to a
    nconc (loop with b = (floor (sqrt (- (* a a) (* i i))))
            for j from (- b) to b collect (list i j)) into points
    finally (format t "~V{~{~D, ~D~}~%~}"
                    n (sort points #'point<))))

円周をたどって。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 << ", "