从Mentgomery到CIOS

从Mentgomery到CIOS

1.重新回忆Montgomery算法

​ Montgomery算法的核心思想是将正常数域上的取模操作，变换到适合于计算机操作的redix-2数域上，再进行取模操作。

step1：数域转换，设原始数据为a，b,映射到Montgomery域

$$\begin{gathered} A=Mont(a,(R^{2}modN),N)=a\ast R^2\ast R^{-1}modN=a\ast RmodN \\ B=Mont(b,(R^{2}modN),N)=b\ast R^2\ast R^{-1}modN=b\ast RmodN \end{gathered}$$

step2：进行Montgomery算法

/*
	Input:A,B,N,R
	Output:A*B*R^(-1)modN
*/
Mont(A,B,N,R):
1.	N' = -N^(-1)modR;
2.	T = A*B;
3.	M = T*N'modR;
4.	T = T + M*N;
5.  T = T/R;
6   if(T>N)return T-N;//约减
7.	else return T;

step3：改进M

​ 上面的计算是要一次性计算出M，这样做的缺点是，M可能会很大，大到难以存储，那么假设在b进制下，就可以将M看作(假设M在2进制下有s位)：
$$\begin{gathered} T=(t_{s-1},t_{s-2},...t_0{)}_2 \\ M=(M_{s-1},M_{s-2},...M_0{)}_2 \\ M=\sum M_i\ast 2^i（注：i\in [0,s-1]） \end{gathered}$$
​ 那么上述第三第四步骤就可以进行合并：
f o r ( i = 0 ; i < = s − 1 ; i + + ) {                                                  m i = t i N ′ m o d R ; T = T + m i 2 i N ; }                                                  for(i=0;i<=s-1;i++)\\ \{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ m_i = t_iN'modR;\\ T = T + m_i2^iN;\\ \}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ for(i=0;i<=s−1;i++){                                                mi​=ti​N′modR;T=T+mi​2iN;}                                                

step4：代换N’

​ 与M同样，N’也存在过宽而导致无法存储的问题，下面规定，一次存储的位宽w=32bit。使用下式替换调N’:
$$\begin{gathered} 令R=2^w; \\ {N}_0^{{}^{\prime }}=N^{\prime }modR或者N_0^{{}^{\prime }}=-N_0^{-1}modR \end{gathered}$$
​ 则Montgomery算法可以变换为：
N 0 ′ = N ′ m o d R                          T = A ∗ B                                f o r ( i = 0 ; i < = s − 1 ; i + + ) {                                                  m i = t i N 0 ′ m o d R ; T = T + m i 2 i N ; }                                                   T = T / R ;                                    i f ( T > N ) r e t u r n   T − N ;        e l s e   r e t u r n T ;                             N_0^{'}=N'modR\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ T=A*B\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ for(i=0;i<=s-1;i++)\\ \{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ m_i = t_iN_0'modR;\\ T = T + m_i2^iN;\\ \} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ T = T/R;\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ if(T>N)return \ T-N;\ \ \ \ \ \ \\ else \ return T;\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ N0′​=N′modR                        T=A∗B                              for(i=0;i<=s−1;i++){                                                mi​=ti​N0′​modR;T=T+mi​2iN;}                                                 T=T/R;                                  if(T>N)return T−N;      else returnT;                           
一致性说明：可以把上述式子中for循环部分展开：
$$\begin{gathered} m_0=t_0{N}_0^{\prime }modR \\ {m}_1=t_1{N}_0^{\prime }modR \\ {m}_2=t_2{N}_0^{\prime }modR \\ ... \\ {m}_{s-1}=t_{s-1}{N}_0^{\prime }modR \\ {T}_0=T+m_0\ast 2^{0}N=(t_0{N}_0^{\prime }modR)2^{0}N=(N_0^{\prime }N{t}_0{2}^0)modR \\ {T}_1=T+(N_0^{\prime }N{t}_0{2}^0)modR+(N_0^{\prime }N{t}_1{2}^1)modR \\ ... \\ {T}_{s-1}=T+(N_0^{\prime }N{t}_0{2}^0)modR+...+(N_0^{\prime }N{t}_{s-1}{2}^{s-1})modR \\ =T+N_0^{\prime }N\sum (t_i{2}^i)(modR) \end{gathered}$$
带入：
$$\begin{gathered} N_0^{{}^{\prime }}=N^{\prime }modR \\ {T}_{s-1}=T+N_0^{\prime }N\sum (t_i{2}^i)(modR) \\ =T+N^{\prime }N(modR)\sum (t_i{2}^i)(modR) \\ =T+N^{\prime }N\sum (t_i{2}^i)(modR) \\ 又T=\sum (t_i{2}^i) \\ 原式=T+N(T\ast N^{\prime }modR)=T+MN \\ 可见与"1"中Mont(A,B,N,R)步骤完全一致 \end{gathered}$$

step5：代换T

​ 上面步骤都需要实现计算T，同样的T=AB也是一个很大的数，同样可以将AB的过程进行拆分，对A*B的拆分比较容易：
KaTeX parse error: Undefined control sequence: \ at position 248: … \ \ \ \ \ \ \ \̲ ̲
最终算法更新为：
N 0 ′ = N ′ m o d R ;                          T = 0 ; / / 初 值 0                                f o r ( i = 0 ; i < = s − 1 ; i + + ) {                                                  T = T + a i b i B ; m i = t i N 0 ′ m o d R ; T = T + m i 2 i N ; }                                                   T = T / R ;                                    i f ( T > N ) r e t u r n   T − N ;        e l s e   r e t u r n T ;                             N_0^{'}=N'modR;\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ T=0;//初值0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ for(i=0;i<=s-1;i++)\\ \{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ T = T + a_ib^iB;\\ m_i = t_iN_0'modR;\\ T = T + m_i2^iN;\\ \} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ T = T/R;\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ if(T>N)return \ T-N;\ \ \ \ \ \ \\ else \ return T;\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ N0′​=N′modR;                        T=0;//初值0                              for(i=0;i<=s−1;i++){                                                T=T+ai​biB;mi​=ti​N0′​modR;T=T+mi​2iN;}                                                 T=T/R;                                  if(T>N)return T−N;      else returnT;                           

2.从Mentgomery到CIOS

​ 按照上述五个步骤，可以对Mentgomery算法进行改正，那么建立在上述五个改进步骤之上，产生了CIOS算法：

输入：A,B,n,n'[0]
输出：T = A*B*R^(-1)(modN)
T = 0;   
(C,S) = 0;
for(i=0;i<=s-1;i++)
{
	C = 0;
	for (j=0;j<=s-1;j++)
	{
		(C,S) = t[j] + a[j]*b[i] + C;
		t[j] = S;
	}
	(C,S) = t[s] + C;
	t[s] = S;
	t[s+1] = C ;
	C = 0;
	m = t[0]*n'[0] mod W;//W = 2^w
	for (j=0;j<=s-1;j++)
	{
		(C,S) = t[j] + m*n[j] + C;
		t[j] = S;
	}
	(C,S) = t[s] + C;
	t[s] = S;
	t[s+1] = t[s+1] + C;
	for (j=0;j<=s;j++)//※※
		t[j] = t[j+1];
}
T = {ts,t(s-1),...,t1,t0};
if(T>N)return T-N;//约减
else return T;
//(version2020.6.3:这里还是存在问题，T和N仍可能很占空间，不可能一次做完)

下面逐块分析CIOS算法，首先明确两个字符C和S：
$$\begin{gathered} C:Carry，进位字符，表示上次计算的进位值。 \\ S:Sum,表示本次加法计算的结果。 \end{gathered}$$

上面提到的CIOS算法需要遵循下面假设：
	1.每次操作使用的变量位宽w = 32bit。
	2.每个数组单元：例如t[i]，指的是当前操作的32bit存储空间，换句话来说"//※※"处的for循环，会导致t[0]被t[1]覆盖，
	等价于对T = {ts,t(s-1),...,t1,t0};一次性右移32bit，
	等价于T/W;其中W = 2^w
	由于"//※※"在内循环j中，而外循环（即变量i的循环）会执行S次，就相当于T会被移动S个W次，即右移2^(sw)个bit，而恰好R = 2^(sw)。完成了T/R的操作。

计算A*B:

for(i=0;i<=s-1;i++)
{
	C = 0;
	for (j=0;j<=s-1;j++)
	{
		(C,S) = t[j] + a[j]*b[i] + C;
		t[j] = S;
	}
	(C,S) = t[s] + C;
	t[s] = S;
	t[s+1] = C ;
}

​ 如上式，如果删去后面的一部分，就可以明显的看出内循环中，第一个for循环是用来计算A*B,即完成Step5所对应的修正。

计算M修正：

	C = 0;
	m = t[0]*n'[0] mod W;//W = 2^w
	for (j=0;j<=s-1;j++)
	{
		(C,S) = t[j] + m*n[j] + C;
		t[j] = S;
	}
	(C,S) = t[s] + C;
	t[s] = S;
	t[s+1] = t[s+1] + C;
	for (j=0;j<=s;j++)//※※
		t[j] = t[j+1];

我将最后一个循环也放在了这个里面，主要是因为，想说明在做最后一个for循环（实质上式在计算T>>32或者说T/(2^32)）时，是如何能保证整除的。
$$\begin{gathered} 定理：假设n和r的gcd(n,r)=1。且满足 \\ r\ast r^{-1}-n\ast n^{\prime }=1 \\ {n}^{\prime }=-n^{-1}modr， \\ 则对所有整数t，当m=t\ast n^{\prime }modr时 \\ (t+m\ast n)/r是一个整数，且满足(t+m\ast n)/r\equiv t\ast r^{-1}modn。 \\ 证明：因为r\ast r^{-1}-n\ast n^{\prime }=1,则 \\ (t+m\ast n)/r=(t+t\ast n^{\prime }\ast n)/r=(t+t\ast (r\ast r^{-1}-1))/r=t\ast r^{-1}modn \end{gathered}$$
​ 则可以看到，每次内循环当j=0时，均有：

(C,S) = t[0] + m*n[0] + C;
	==>带入m
	==>t[0] + t[0]*n'[0]*n[0] + C;
	==>将t[0]+C整体看做定理中所述的整数t，而定理中的r = 2^w = 2^32
	==>显然t[0]是可以被2^32整出的

又由于外循环(i循环的存在)，最终每次增加的C会被加入到T中去，进而完成了：

4.	T = T + M*N;
5.  T = T/R;

这两个步骤。