Tuesday, January 13, 2009

2平方数の和

ある素数pを
p = a^2 + b^2
の形に表そうという問題があります。
p = 1 (mod 4)の時に解があるということは知っていましたが、解き方はしりませんでした。総当たりで解いていたのです。どうやらこれは非常に効率の悪い方法だったようです。

本ってすばらしいですね。その方法が書いてありました。a^2+b^2=kpとなるa,bからスタートしてkを1まで減らすという作戦です。

1. a^2 = -1 (mod p) となるaを求める。
 フェルマーの定理と乱数を利用。

2. a^2 + 1^2 = 0 (mod p) なので、b = 1とおけば、a^2+b^2 = kp。k=1なら終了。

3. u = a (mod k), v = b (mod k)とする。u^2+v^2 = a^2+b^2 (mod k)となる。
  q = (u^2+v^2)/kとする。
 
4. 定義より(u^2+v^2)(a^2+b^2) =qkkp
  また式の変形で、(u^2+v^2)(a^2+b^2) = (ua+vb)^2 + (va-ub)^2

5. それぞれの項がkで割れる筈なので
a = (ua+vb)/k
b = (va-ub)/k
k = q
として2に戻る。

こんな計算方法、よく見つけたなあ。すごい。
これはかなり速いんですね。以下、私のコード。


from random import randint

def square_mod(p):
while True:
a = randint(2, p - 1)
b = pow(a, (p - 1) / 4) % p
if (b*b) % p == p - 1:
return b

def two_square_sum(p):
if p % 4 != 1: return [0, 0]

a = square_mod(p)
b = 1
k = (a * a + b * b) / p

while k > 1:
u, v = a % k, b % k

# abs(u) <= k/2
# abs(v) <= k/2
if u > k / 2: u -= k
if v > k / 2: v -= k

# u * u + v * v = 0 (mod k)
q = (u * u + v * v) / k

# (u * u + v * v) * (a * a + b * b) = (u * a + v * b)^2 + (v * a - u * b)^2 = (q * k) * (k * p)
c, d = (u * a + v * b) / k, (v * a - u * b) / k

# next turn
a, b = c, d
k = q

return [abs(a), abs(b)]


実際にやってみましょう。


for i in xrange(5, 10000, 4):
if is_prime(i):
x, y = two_square_sum(i)
print i, x, y

5 2 1
13 3 2
17 4 1
29 5 2
37 6 1
41 5 4
53 7 2
61 6 5
73 8 3
89 5 8
97 9 4
101 10 1
109 10 3
113 8 7
137 11 4
.....


この本で読みました。

Sunday, January 04, 2009

組み合わせを作るアルゴリズム#3

以前、k-subsetのことを書きました。Pythonで書きましたがC言語に直訳できそうな書き方でした。私のPythonの書き方はいつもそれで、たぶんPythonらしくないのでしょう。ダサイ書き方なんですよね。昨日某所にて同じことを以下のように書いているのをみてそう感じました。


def ksubset2(n, k, pos):
if k > 0:
for i in xrange(pos, n - k + 1):
for s in ksubset2(n, k - 1, i + 1):
yield [i] + s
else:
yield []


すごいなあ。格好いい。私がやってもこういう風には絶対ならないですね。これはジェネレータなのでksubset2(5,3,0)に前と同じ組み合わせが入っていることになります。

しかし、自分でやるとしたら、やっぱりCに置き換えやすい書き方をしておきますかね。Pythonで遅ければCに直すという理由で。つまりジェネレータは使わないのです。

ちなみに、某所でみたのは上のようにnで上限を渡すのではなくて、同じことではありますが以下のようにseqに元の集合を入れて渡すものでした。

def ksubset3(seq, k):
if k > 0:
for i in xrange(0, len(seq) - k + 1):
for s in ksubset3(seq[i+1:], k - 1):
yield [seq[i]] + s
else:
yield []