Saturday, October 15, 2011

2平方数の和(2)

nがあたえられたときに、x^2+y^2=nを解くにはどうしたらよいか。nが素数pの場合については以前書いた。では合成数の場合はどのようにすればよいか。

n=pqのとき、以前の方法でp=u^2+v^2, q = a^2+b^2ということが分かっているものとして、つぎの恒等式をつかってnを平方数の和に分ける。

n = pq = (u^2+v^2)(a^2+b^2) = (ua+vb)^2 + (va-ub)^2

これでnはua+vbとabs(va-ub)に分けられた(どちらも正数だけ採用することにする)。これを繰り返してゆけばnが2個以上の素数で分解出来る場合も対処できる。

元々、pが4で割ったとき1余る場合、x^2+y^2=pの解はあり、pの余りが3の場合に解はない。nにそのようなpが含まれている場合、pの個数が奇数のとき、x^2+y^2=nに解がなくなる。偶数の場合には解ありである。例えば、n=9の時、p=3だが3が2個あるのでx^2+y^2=1の解 (1,0)を3倍した(3, 0)が解となる。また、n=3やn=27には解がない。

コードにしたのが以下。素因数分解は状況によって違うものを使いたい(試し割りのものと、フルイをつかったもの)ので、初期化で渡すようにした。以下ではfacs()がそれで、これに期待しているのは、n=p^r * q^sのとき、[p,r],[q,s]を渡してくれること。fscsの中で使っているfactors()が素因数分解のコア部分。

名前をprocess_primeに変えたが、これが以前書いた、素数のx^2+y^2=pを解くところ。comb()が上に書いた恒等式で展開するもの。

class square_sum(object):
  def __init__(self, func):
    self.func = func
  def square_mod(self, p):
    from random import randint
    while True:
      a = randint(2, p - 1)
      b = pow(a, (p - 1) / 4) % p
      if (b*b) % p == p - 1:
        return b

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

    a = self.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)] 
  def comb(self, X, Y):
    u, v = X
    a, b = Y

    res = set()
    r1 = [u * a + v * b, abs(v * a - u * b)]
    r1.sort()
    res.add(tuple(r1))

    r2 = [u * b + v * a, abs(v * b - u * a)]
    r2.sort()
    res.add(tuple(r2))

    return res

  def times(self, M, seq):
    return map(lambda x: x*M, seq)

  def solv(self, n):
    fac = self.func(n)

    # p%4==3 prime
    M = 1
    for p,r in fac:
      if p % 4 == 3:
        if r % 2 == 0:
          M *= p**(r/2)
        else:
          return []

    # p==2 or p%4==1 prime
    seq = []
    for p,r in fac:
      if p == 2 or p % 4 == 1:
        sq = self.process_prime(p)
        seq += [sq]*r

    if len(seq) ==  0:
      return [self.times(M, [0,1])]

    res = set()
    res.add(tuple(seq[0]))
    for k in xrange(1, len(seq)):
      new = set()
      for e in res:
         new |= self.comb(e, seq[k])
      res = new

    return map(lambda x: self.times(M,x), res) 
# if n=p1^r1 x ... x pm^rm, returns [p1,r1]...[pm,rm]
def facs(n):
  seq = factors(n)
  s = set(seq)
  r = []
  for p in s:
    r.append([p, seq.count(p)])
  return r
使い方は、
sq = square_sum(facs)
sq.solv(1234)
のような感じ。
比べてみると、しらみつぶしでx^2+y^2=nを解くよりはかなりはやい。




1 comment:

Anonymous said...

合成数の平方数の和を探していて、ここにたどり着きました。こんなに簡単にできるのですね。参考になりました(Chouei)