ストラッセン-ショーンハーゲ法 (Modified:2005-04-20)

多倍長整数の高速乗算法であるストラッセン-ショーンハーゲ法の実装法について検討する。作譜に必要な個々の算法について見渡した後、32ビット語長の場合について、その主要部分を算譜としてまとめる。

　整数n,m>0及び、R=Z/mZ上の1の原始n乗根ζを適当に定めることによって、R上の離散Fourier変換を定義することができる。([F1]参照) ストラッセン-ショーンハーゲ法においては、K=2^kとして、

(n,m,ζ)

＝

(K,

2^2K+1,

2⁴)

という離散Fourier変換を使用する。この変換はK=2^kであることから、Cooley-Tukey型の高速算法が適用できることと、1の原始K乗根ζ=2⁴を乗ずることが、乗算命令ではなくて、2進演算では非常に効率がいいビットシフト命令で実行することができることから効率の良い最適化が期待できる。まず、これが離散Fourier変換の必要条件を満たしていることを証明する。

補題1.1
　K=2^k(k>0)として、

(n,m,ζ)

＝

(K,

2^2K+1,

2⁴)

とするとき、整数の剰余環R=Z/mZは、次の離散Fourier変換の必要条件を満たす。

ζは1の原始K乗根
∃K^-1∈R=Z/mZ
K>l>0 → ∑_i=0^K-1ζ^il＝0

証明　1.は次の式による。

ζ^K/2＝2^2K≡-1

(mod 2^2K+1),

　ζ^K≡1

(mod 2^2K+1)

2.は(K, 2^2K+1)=1による。3.については、まず、

≠

(mod K)

⇒

(ζ^l-1,

2^2K+1)

となることを示す。

ζ^l

≡

(mod p),

2^2K

≡

-1

(mod p)

なる素数pが存在したと仮定すると

ζ^l

＝

2^4l

≡

(mod p)

より、4K|4lからK|lとなりl≡0(mod K)が導かれる。

ζ^l

K-1
Σ
i=0

ζ^il

K-1
Σ
i=0

ζ^(i+1)l

≡

K-1
Σ
i=0

ζ^il

ζ^Kl

ζ^0l

(mod 2^2K+1)

≡

K-1
Σ
i=0

ζ^il

(mod 2^2K+1)

より、

(ζ^l-1)

K-1
Σ
i=0

ζ^il

≡

0	(mod 2^2K+1)

したがって(ζ^l-1, 2^2K+1)=1より、

K-1
Σ
i=0

ζ^il

≡

0	(mod 2^2K+1)

■

仮定1.2 (多倍長整数の高速乗算法)
　

＝

2^k

0),

＝

K²

とし、正の整数a,bが、基数P=2^Kにより次のようにP進表現されるものとする。

＝

a_K-1P^K-1

…

a₁P

a₀,

0≦a_i＜P

＝

b_K-1P^K-1

…

b₁P

b₀,

0≦b_i＜P

この場合、0≦a,b≦2^N-1であるが、さらに

0≦ab≦2^N-1

と仮定する。

　仮定より、abは2^N-1を法として求めればよい。P^K≡1(mod 2^N-1)より、

＝

≡

K-1
Σ
r=0

(

Σ
s+t≡r(mod K)

a_sb_t

)P^r

(mod 2^N-1)

となる。この係数c_r=∑_{s+t≡r(mod K)}a_sb_tの大きさを評価すると、

≦

c_r

≦

K(2^K-1)²

＜

K(2^2K+1)

となる。したがってc_rはK(2^2K+1)を法として求めればよい。(K, 2^2K+1)=1であるから、孫子の剰余定理から次のことがわかる。

補題1.3
　c_rは、c'_r及びc''_rが既知であるとして、連立合同式

c_r

≡

c'_r

(mod K)

c_r

≡

c''_r

(mod 2^2K+1)

を解くことによって求めることができる。

算法1.4
　c'_rは、次の手順で計算する。

a'_r≡a_r (mod K)
b'_r≡b_r (mod K)
c'_r≡∑_{s+t≡r(mod K)}a'_sb'_t (mod K)

算法1.5
　c''_rは次のように計算する。

F(a)_i

≡

K-1
Σ
s=0

a_sζ^si

(mod 2^2K+1)

F(b)_i

≡

K-1
Σ
t=0

b_tζ^ti

(mod 2^2K+1)

c''_r

≡

K^-1

K-1
Σ
i=0

F(a)_iF(b)_iζ^-ir

(mod 2^2K+1)

　孫子の剰余定理は、３世紀後半に著された「孫子算経」に連立合同式を満たす整数を求める問題があることから付けられた現代中国で使われている名称である。欧米では中国の剰余定理(Chinese remainder thorem)と呼ばれている。([A2]参照)

定理2.1 (孫子の剰余定理)
　正の整数m₁,m₂,…,m_kはどの２つも互いに素であると仮定する。このとき、任意の整数u₁,u₂,…,u_kに対して、連立合同式：

u≡u₁

(mod m₁),

u≡u₂

(mod m₂),…,

u≡u_k

(mod m_k)

を満たす整数uが存在し、m₁m₂…m_kを法として一意的に定まる。したがって、特に

≦

＜

m₁m₂…m_k

を満たすものはただ一つ存在する。

証明　一意性をまず示す。整数u,vが連立合同式を満たすとすると

≡

(mod m_j),

1≦j≦k

となり、u-vはすべてのm_jの倍数であり、したがってm=m₁m₂…m_kの倍数となる。すなわち、

≡

(mod m₁m₂…m_k)

である。存在することを構成的に証明する。

M_j

＝

(m/m_j)^φ(m_j),

1≦j≦k

とおく。ここでφ(m_j)はオイラーの関数である。そうするとオイラーの定理により

M_j

≡

(mod m_j)

M_j

≡

(mod m_k),

k≠j

が成り立つ。このとき、

k
Σ
i=1

u_iM_i

はすべての合同式を満足する。 ■

算法2.2 (孫子の剰余定理の解の高速算法)
　まず、ユークリッドの互除法により、

c_ijm_i

≡

(mod m_j),

1≦i＜j≦k

なるc_ijを計算しておく。これは

am_i

bm_j

＝

においてc_ij=aとすればよい。このとき、v_r, 1≦r≦kを

v₁

＝

u₁

m₁

v₂

＝

(u₂-v₁)c₁₂

m₂

：

v_r

＝

(…((u_r-v₁)c_1r-v₂)c_2r-

…

-v_r-1)c_(r-1)r

m_r

とする。このとき、

＝

v_km_k-1…m₁

…

v₃m₂m₁

v₂m₁

v₁

が解である。

証明　

v_jm_j-1…m₁

≡

m_j-1…m₁(…((u_j-v₁)c_1j-v₂)c_2j-

…

-v_j-1)c_(j-1)j

≡

m_j-2…m₁(…((u_j-v₁)c_1j-v₂)

c_2j-

…

-v_j-2)c_(j-2)j

v_j-1m_j-2…m₁

：

≡

u_j-v₁-v₂m₁-

…

-v_j-1m_j-2…m₁

(mod m_j)

■

算譜2.3 (孫子の剰余定理の解法)
　

integer c[k][k], m[k], u[k], v[k];
integer i, r, u;

// v[k]の計算
for (r = 0; r < k; r++)
{ i = 0;
  tmp = u[r];
  while (i < r)
  { tmp = (tmp - v[i])*c[i][r];
    i++;
  }
  v[r] = tmp % m[r];
}

// 解uの計算
u = 0;
for (r = 0; r < k; r++)
{ i = 0;
  tmp = v[r];
  while (i < r)
  { tmp *= m[i];
    i++;
  }
  u += tmp;
}

　ユークリッドの互除法は最大公約数を求める算法であり、「原論」のⅦ巻に定理として述べられている。近代の多くの初等整数論の教科書においても、そのまま最大公約数を求める算法として解説されているが、これが必ずしも最良の算法であるとはいえない。ここでは最大公約数を求める算法をある種の有限列を計算する方法として形式化し、その形式でユークリッドの互除法を表現した後、J. Steinにより発見された除算命令を使わない算法について述べる。([S2]参照)

定義2.4 (最大公約数の算法のベクトル列表現)
　２つの負でない整数a,bについて、

最大公約数d=(a,b)を求める
不定方程式d=ma+nbを満たすm,nを求める

算法について考える。

を満たす整数ベクトル(x,y,z)の組の有限列{(x_i,y_i,z_i),(s_i,t_i,u_i)}_i=0ⁿが、次の条件を満たすものとする。

z₀＝a, u₀＝b
z_i＞u_i≧0, i=1,…,n
max(z_i,u_i)＞max(z_i+1,u_i+1), i=0,…,n-1
(z_i,u_i)＝(z_i+1,u_i+1), i=0,…,n-1
z_n＝d, u_n＝0

もしこのような有限列を何らかの手段で生成することができれば最大公約数d=z_n及び不定方程式の解m=x_n,y=z_nが求まる。これを最大公約数の算法のベクトル列表現と呼ぶ。

算法2.5 (拡張されたユークリッドの互除法)
　２つの負でない整数a,bが与えられたとき、

(x₀,y₀,z₀)＝(1,0,a),(s₀,t₀,u₀)＝(0,1,b)
u_i≠0ならば
(x_i+1,y_i+1,z_i+1)＝(s_i,t_i,u_i)
(s_i+1,t_i+1,u_i+1)＝(x_i,y_i,z_i)-(s_i,t_i,u_i)z_i/u_i

によって生成されるベクトル列は、最大公約数の算法のベクトル列表現である。

算譜2.6 (拡張されたユークリッドの互除法)
　

integer a, b;
integer x, y, z, s, t, u;
integer q, d, e, f;

x = 1; y = 0; z = a;
s = 0; t = 1; u = b;
while (u != 0)
{ q = z/u;
  d = s; e = t; f = u;
  s = x - s*q; t = y - t*q; u = z - u*q;
  x = d; y = e; z = f;
}

算法2.7 (J. Steinによる最大公約数の算法)
　二つの負でない整数a,bが与えられたとき、もしともに2を素因数としてもつならば、必要なだけ2のべき乗で割ることによって、2を共通の約数としてもたないようにしておく。このとき、

(x₀,y₀,z₀)＝(1,0,a),(s₀,t₀,u₀)＝(b,1-a,b)
z_iが偶数のとき、
- x_i,y_iがともに偶数のとき、
  (x_i+1,y_i+1,z_i+1)＝(x_i,y_i,z_i)/2
- x_i,y_iのいずれかが奇数のとき、
  (x_i+1,y_i+1,z_i+1)＝(x_i+b,y_i-a,z_i)/2
u_iが偶数のとき、
- s_i,t_iがともに偶数のとき、
  (s_i+1,t_i+1,u_i+1)＝(s_i,t_i,u_i)/2
- s_i,t_iのいずれかが奇数のとき、
  (s_i+1,t_i+1,u_i+1)＝(s_i+b,t_i-a,u_i)/2
z_i,u_iがともに奇数のとき、
(x_i+1,y_i+1,z_i+1)＝(x_i,y_i,z_i)-(s_i,t_i,u_i)
(s_i+1,t_i+1,u_i+1)＝(s_i,t_i,u_i)
もし2.3.4.の操作で、z_i+1＜u_i+1となればベクトルの順序を入れ替える。

によって生成されるベクトル列は、最大公約数の算法のベクトル列表現である。

証明　最大公約数が求まることは、aが偶数でbが奇数のとき(a,b)=(a/2,b)であることと(a,b)=(a-b,b)からわかる。2.と3.の操作で、ax+by=zが偶数で、x,yがともに偶数でないときx+bとy-aがともに偶数になることを示す。xが奇数であるとする。aが奇数のとき、bとyはともに奇数でなければならない。aが偶数のとき、仮定からbは奇数である。したがってyは偶数である。yが奇数である場合も同様である。 ■

算譜2.8 (J. Steinによる最大公約数の算法)
　

integer a, b;
integer x, y, z, s, t, u;
integer d, e, f;
integer k;

k = 0;
while (a,bはともに偶数)
{ a >>= 1; b >>= 1;
  k++;
}

x = 1; y = 0; z = a;
s = 0; t = 1; u = b;
if (u > z)
{ d = s; e = t; f = u;
  s = x; t = y; u = z;
  x = d; y = e; z = f;
}
while (u != 0)
{ if (zは偶数)
  { if (x,yはともに偶数)
    { x >>= 1; y >>= 1; z >>= 1;
    } else {
      x += b; y -= a;
      x >>= 1; y >>= 1; z >>= 1;
    }
  }
  else if (uは偶数)
  { if (s,tはともに偶数)
    { s >>= 1; t >>= 1; u >>= 1;
    } else {
      s += b; t -= a;
      s >>= 1; t >>= 1; u >>= 1;
    }
  }
  else
  { x -= s; y -= t; z -= u;
  }
  if (u > z)
  { d = s; e = t; f = u;
    s = x; t = y; u = z;
    x = d; y = e; z = f;
  }
}
z <<= k;

　ストラッセン-ショーンハーゲ法の中に2^2K+1を法とする演算がでてくるが、ここでは一般に2^e+1を法とする加法と乗法を2^eを法とする加法と乗法を利用して計算する方法についてまとめておく。証明は自明なので省略する。

命題2.9
　0≦u,v＜2^e+1とするとき、

(u+v)

(2^e+1)

＝

u+v,

　u+v

＜

2^e+1

＝

(u+v)

2^e

　otherwise

(2^e+1)

＝

2^e

2^e,

2^e

≧

＝

2^e

　otherwise

仮定3.1 (多倍長整数の高速乗算法(32ビット語長))
　

K＝2⁵＝32,

N＝K²＝1024

とし、正の整数a,bが、基数P=2³²により次のようにP進表現されるものとする。

＝

a_K-1P^K-1

…

a₁P

a₀,

0≦a_i＜P

＝

b_K-1P^K-1

…

b₁P

b₀,

0≦b_i＜P

この場合、0≦a,b≦2^N-1であるが、さらに

0≦ab≦2^N-1

と仮定する。

　a_u,b_vを0でない最大の係数としたとき、0≦ab≦2^N-1の判定を次の条件によって行なう。

(u+1)

(v+1)

≦

K-1

この条件はきびしすぎるが、正確な判定を行なうにはかなりの計算量が必要になるので簡単のためにこうしておく。

　c'_r,c''_rが与えられたときに、補題1.3の連立合同式を解いてc_rを求める部分をまとめておく。算法2.2を適用するとm₁=2⁵,m₂=2⁶⁴+1である。まず、

c₁₂m₁≡1

(mod m₂)

なるc₁₂を求める。

2⁵⁹2⁵≡-1

(mod 2⁶⁴+1)

の両辺を２乗することによって、

2¹²³2⁵≡1

(mod 2⁶⁴+1)

c₁₂は2⁶⁴+1を法として求めればよいので、

c₁₂

＝

2¹²³

＝

2⁵⁹(2⁶⁴+1)

2⁵⁹

≡

2⁶⁴

2⁵⁹

(mod 2⁶⁴+1)

≡

2⁵⁹31

(mod 2⁶⁴+1)

となる。

ストラッセン-ショーンハーゲ法

梅谷武

目　　次

ストラッセン-ショーンハーゲ法

梅谷 武

目 次

梅谷武

目　　次