[빠른곱셈-2] 쇤하게-슈트라센 (Schönhage–Strassen) 알고리즘 - 1

원리

빠른곱셈-1에서 했던 FFT연산과 방법은 동일합니다

단, DFT를 정수범위에서 수행합니다

 

숫자 변환

길이가 n인 두 숫자 u,v를 생각해봅시다

$u = \sum_{i=0}^{n-1}u'_i2^i$

$v = \sum_{i=0}^{n-1}v'_i2^i$

이를 l의 길이를 가진 블록 b개로 쪼개줍니다

$u = \sum_{i=0}^{b-1}u_i2^{il}$

$v = \sum_{i=0}^{b-1}v_i2^{il}$

여기서

$n=2^k$

$b=\lfloor 2^{k/2} \rfloor$

$l=\lceil 2^{k/2} \rceil$

이다 (b*l=n임)

 

곱셈

${y=u*v = \sum_{j=0}^{2b-2}y_j2^{lj} \\
y_j =\begin{Bmatrix}
\sum_{i=0}^{j}u_iv_{j-i} &&if,\, (0\leq j <b-1) \\ 
\sum_{i=j-b+1}^{b-1}u_iv_{j-i} &&if,\, (b-1\leq j <2b-1) \\ 
\end{Bmatrix}}$

 

Modulo

곱셈의 결과 y에 $2^n+1$을 mod해주자

$2^n=-1(mod\,2^n+1)$이므로

${\begin{align*} y=u*v &= \sum_{j=0}^{2b-2}y_j2^{lj} (mod\,2^n+1)\\
&=\sum_{j=0}^{b-1}y_j2^{lj}+\sum_{j=0}^{b-1}y_{j+b}2^{lj+lb} (mod\,2^n+1) \\
&=\sum_{j=0}^{b-1}y_j2^{lj}-\sum_{j=0}^{b-1}y_{j+b}2^{lj} (mod\,2^n+1) \\
&=\sum_{j=0}^{b-1}(y_j-y_{j+b})2^{lj} (mod\,2^n+1)\\
&=\sum_{j=0}^{b-1}(z_j)2^{lj} (mod\,2^n+1)\\
(z_j=y_j-y_{j+b}, lb=n)\end{align*}}$

참고로, 숫자를 맞추기 위해 $y_{2b-1}$까지 고려한것이며 반드시 $y_{2b-1}=0$ 입니다

 

z는 u와v의 negative wrapped convolution입니다

$z_j=y_j-y_{j+b}=\sum_{i=0}^{j}u_iv_{j-i}-\sum_{i=j+1}^{b-1}u_iv_{j+b-i}$

 

영상에서는 z를 쓰는데, 왜 z를 도입하는 지 잘 모르겠습니다.

실제로는 z를 거치지 않고 바로 y를 구하면 될  것 같습니다.

 

z 범위

z값의 범위를 알아보자

$z_j=y_j-y_{j+b}$

$z_j$의 최대값는 $y_j$의 최대값에 의존한다

$z_j$의 최소값는 $y_{j+b}$의 최소값에 의존한다

 

${y=u*v = \sum_{j=0}^{2b-2}y_j2^{lj} \\
y_j =\begin{Bmatrix}
\sum_{i=0}^{j}u_iv_{j-i} &&if,\, (0\leq j 
\sum_{i=j-b+1}^{b-1}u_iv_{j-i} &&if,\, (b-1\leq j <2b-1) \\ 
\end{Bmatrix}}$

이고,

$u_av_b<2^{2l}$이다

따라서,

$-[(b-1)-(b+j-b+1)+1]2^{2l}<z_j<(j-0+1)2^{2l}$

최소값을 구할 때는 $j->b+j$로 넣어 구해야 함을 주의하자

 

$z_j$가 될 수 있는 값의 범위는 $b2^{2l}$이다 (최대값 빼기 최소값)

즉, $z_j(mod\,b(2^{2l}+1))$를 해주어도 정보의 손실이 없다

b는 2의 지수, $(2^{2l}+1)$는 홀수이기 때문에 서로소이므로 chinese remainder theorem을 적용할 수 있다

 

Mod $(2^{2l}+1)$와 Mod b를 따로따로 수행해준 뒤 합쳐주면 됩니다

 

Mod 2^(2l)+1

$z_j(mod\,(2^{2l}+1))$를 계산해봅시다

$2^{2l}+1=w^{b/2}+1$ if $w=2^{4l/b}$

$z_j(mod\,(2^{2l}+1))=z_j(mod\,(w^{b/2}+1))$

이고, 이는 primitive $b$-th root of unity와 같습니다

(앞서 벡터의 크기를 b 라고 정의했습니다)

 

FFT를 적용하여 $z_j(mod\,(2^{2l}+1))$를 구할 수 있습니다

 

Mod b

z는 u와v의 negative wrapped convolution입니다

$u= (u_0,u_1,...,u_{b-1}) $

$u_i'=u_i(mod\,b)$ 라고 합시다

$u_i$를 2진수로 표현했을 때 크기는 최대$logb$ 입니다 (정확하게는 $log_2b$)

$u''=u_0'000u_1'000...000u_{b-1}'$

$u_i'$의 앞 뒤로 zero padding을 추가한 숫자를 $u''$ 라고 합시다

zero padding의 크기는 $2logb$로 정합니다

그러면,

$u''*v''$를 계산할 때 정보의 손실이 일어나지 않습니다

 

왜냐하면,

$u''*v''=(u_0'*v_0')000...000\sum{u_i*v_{i-j}}000 ...000(u_{b-1}'*v_{b-1}')$

$\sum{u_i*v_{i-j}}$에서 

$u_i<b$ 이고 $v_i<b$ 인데, 이것이 $\sum$에 의해 최대 b번 더해집니다

따라서,

$\sum{u_i*v_{i-j}}<b^3$이고 이것은 최대 $3logb$의 크기를 가질 수 있습니다

따라서, zero padding을 $3logb$빼기 $logb$($u_i$의 크기)로 정하면 정보의 손실이 없습니다

$u''*v''$를 구한 후 비트 이동만 시켜주면 $y_i(mod\,b)$를 구할 수 있습니다

 

이때의 시간복잡도는 u'',v''의 크기가 3blogb이므로,

M(3blogb)입니다. 카라추바 곱셈을 활용하면, (3blogb)^(1.585)가 되고 b는 n^(1/2)이기 때문에,

이때의 시간복잡도는 n보다는 작습니다

 

Mod b(2^(2l)+1)

chinese remainder theorem으로 $z_j(mod\,b(2^{2l}+1))$을 계산하겠습니다

앞선 포스팅에서,

$u(mod\,p)=\sum_{j=0}^{n-1}u_jc_jd_j(mod\,p)$

$u(mod \,p_i) = u_i$

이었습니다

 

$p_0=b$

$p_1=2^(2l)+1$

$z_i'=z_i(mod\,b)$

$z_i''=z_i(mod\,2^{2l}+1)$

라고 합시다

c의 정의에 의해 $c_1 = p_2$, $c_2= p_1$입니다

 

d의 정의에 의해

$p_1*d_0(mod\,b)= 1$

$p_0*d_1(mod\,2^{2l}+1)= 1$

이때, $d_0=1$, $d_1=2^{2l}+1-2^{2l}/b$를 대입하면 식이 성립합니다.

 

${\begin{align}
z_i&=z_i'*(2^{2l}+1)*1+z_i''*(b)*(2^{2l}+1-2^{2l}/b) (mod\,b(2^{2l}+1)) \\
&=z_i'*(2^{2l}+1)*1+z_i''*(b*(2^{2l}+1)-(2^{2l}+1)+1)(mod\,b(2^{2l}+1))  \\
&=(2^{2l}+1)(z_i'-z_i'')+z_i''b(2^{2l}+1)+z_i''(mod\,b(2^{2l}+1)) \\
&=(2^{2l}+1)(z_i'-z_i'')+z_i'' (mod\,b(2^{2l}+1))\\
\end{align}}$

 

다음할 내용

- 시간복잡도 계산 (FFT의 시간복잡도를 구하는 것을 이해 못함)

- python 및 C로 코드 구현하여 built-in 곱셈(카라추바)와 시간복잡도 비교

 

 

 

 

출처

 

https://www.youtube.com/watch?v=OGUMsBkZqkc&list=PLbMVogVj5nJTA6ZeVkqHEiUup_3ccdg_C&index=26