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


この本で読みました。

2 comments:

mamiko said...

突然お邪魔します。整数2009を平方数a+平方数bで表すことができるでしょうか。
できるならば、その二つの平方数を応えなさい。できないならば、その理由を説明しなさい。・・・これは攻玉社中の算数単科入試で出た今年の問題です。
息子は平方数を適当に組み合わせて勘があたったようで20秒くらいで解いたそうです。
しかし、こんな解法があったのですね。

Tsuchiya Yoshihiro said...

はじめまして。
28*28 + 35*35 = 2009
ですね。
しかし、入試の問題としては、かなりキツいですねえ。