C++ RSA

首先要说证明不全是我自己证明的，有些是从各个地方整理出来的，证明包含了每个细节，可以说是最详细的一套从头到脚的证明，有助于对RSA算法的绝对理解，主要是关于a,b互质的不定方程用辗转相除求解，其中还有一些相关的定理证明，还有欧拉定理用于RSA公式的证明。代码是自己写的非高精RSA实现，高精看着太难受了，写着也难受，就暂时先放掉，也就是说任何一个地方都不能超过long long 的取值范围

#include<iostream>
#include<cstdlib>
#include<cmath>
#define ull long long
using namespace std;

ull Euclidean(ull a,ull b){
	if(a==0){
		return 0;
	}
	if(b==0){
		return a;
	}
	Euclidean(b,a%b);
}

//这个方法用于a,b互质，不互质的数可以转化为互质的两个数，可以在外面包装一层处理
void EuclideanEx(ull a,ull b,ull& x,ull& y){
	if(a==0){
		return;
	}
	if(b==0){
		x=1;
		y=0;
		return;
	}
	EuclideanEx(b,a%b,x,y);
	ull tx=x;
	x=y;
	y=tx-a/b*y;
}

//b>0
ull Positive(ull x0,ull b){
	if(x0<0){
		ull nx0=x0*-1;
		ull t=1;
		if(b<nx0){
			t=ceil(nx0/b);
		}
		x0=x0+b*t;
	}
	return x0;
}

ull MontgomeryMod(ull base,ull p,ull n){
	ull ret=1;
	ull t=base;
	while(p>0){
		if(p&1){
			ret=ret*t%n;
		}
		t=t*t%n;
		p>>=1;
	}
	return ret;
}

ull RSAEncrypt(ull data,ull n,ull e){
	return MontgomeryMod(data,e,n);
}
ull RSADecrypt(ull endata,ull n,ull d){
	return MontgomeryMod(endata,d,n);
}

int main(){
	srand(time(NULL));
	ull e,p,q,x,y,n,fn;
	cin>>p>>q;
	cout<<Euclidean(p,q)<<endl;
	n=p*q;
	fn=(p-1)*(q-1);
	do{
		e=rand()%fn+1;
	}while(Euclidean(fn,e)!=1);
	EuclideanEx(e,fn,x,y);
	x=Positive(x,fn);
	cout<<x<<endl;
	ull m;
	cin>>m;
	m=RSAEncrypt(m,n,e);
	cout<<m<<endl;
	m=RSADecrypt(m,n,x);
	cout<<m<<endl;
	return 0;
}

a*x+b*y==gcd(a, b) 与 a/gcd(a,b)*x+b/(gcd)*y=1 在有整数解的时候是等价的 
所以只要证明：若a,b互质,a*x+b*y=1即证明了上式


a,b互质，a*x+b*y=1，一定存在整数x,y使得等式成立。关于这个的证明：
因为a,b为正整数，所以a+b>=2，x与y异号，令m=-y => a*x-b*m=1; m=(a*x-1)/b


a,b互质，a,2a,3a,...,(b-1)a,ba 除以b的余数包含 0,1,2,3,...,b-2,b-1
[关于这个的证明，主要是证明其间每个数的余数都不同：若存在i,j 0<i<j<b,使得ia与ja的余数相同，那么(j-i)*a % b=0，因为 i,j<b 且 a,b互质，所以(j-i)*a % b=0不可能成立，所以不存在这样的i,j]
顺带提一下，欧拉定理的证明第一步，也会用到这个地方，就先说下关于欧拉定理的证明：a^φ(b) % b =1
1*a*2*a*3*a*...*(b-1)*a % b = (1*a)%b *(2*a)%b *(3*a)%b *... *(b-1)*a%b 这一步是模乘法运算公式，证明很简单，则排序后是 =1*2*3*...*(b-1) % b
a^φ(b)*(b-1)!=(b-1)!
a^φ(b)=1


在其中选取k,x=k使得a*k % b=1, m=(a*k-1)/b =>m一定为整数，那么x=k,y=-m时，a*x+b*y=1成立
那么从这也可以看出，通过枚举是可以求出解的


接下来再说下辗转相除：
a=m*gcd(a,b);b=n*gcd(a,b),r=a%b=?*gcd(a,b);r2=b%r=??*gcd(a,b);r3=r%r2=???*gcd(a,b).....r?%r??=0 
可见余数里始终含有gcd(a,b)这个因子，那么辗转的过程就是再不断缩小余数里另外个因子，因为gcd(a,b)是最大公因数，1又是任何整数的公因数，即当r??=1*gcd(a,b)时，结束这个过程


扩展欧几里德定律正是利用了这个性质，我们都认为a,b互质，gcd(a, b)==gcd(b, a%b)==1
所以a*x1+b*y1=b*x2+(a%b)*y2;   　　因为r=a%b,   r =a-k*b所以==>
a*x1+b*y1=b*x2+(a-k*b)*y2;       　　因为k=a/b;所以　==>
a*x1+b*y1=b*x2+(a-(a/b)*b)*y2;   　　展开得到　　==>　　　　
a*x1+b*y1=b*x2+a*y2-b*(a/b)*y2;　　转换得到      ==>
a*x1+b*y1=a*y2+b*(x2-(a/b)*y2);
观察上式可知 x1=y2, y1=x2-a/b*y2;
通过不断的递归，yn=0时结束，此时xn=gcd(a,b)=1。


如果d小于0，那么还需要得出不定方程的通解，使用d大于0的一个解。
如果已知x0,y0为不定方程a*x+b*y=c的解，那么
x=x0+b*t; 因为b一定是大于0的，我这就
y=y0-a*t;


接下来就是RSA公钥和私钥的生成算法：
1：随机选择两个不相等的质数m、n
2：计算欧拉函数φ(n)=(m-1)*(n-1) 表示小于n与n互质的数有多少个
3：选取随机数e，1<e<φ(n)，且e与φ(n)互质 （一般选取65537）
4：计算e对于φ(n)的模反元素d
e*d % φ(n)=1
e*d-kφ(n)=1
e和φ(n)已知，相当于解 a*x+b*y=1，并且解得的x 要使得 e*x % φ(n) =1 所以x必须是正整数
5：公钥 (n,e) 私钥 (n,d)


加密过程：
加密信息m必须小于n
m^e % n =c


解密过程：
c^d % n =m


证明：
m^e%n=c=m^e-kn
已知ed%φ(n)=1,c^d%n=(m^e-kn)^d%n=m^(ed)%n=m^(h*φ(n)+1)%n
1.若m,n互质，欧拉定理，m^(φ(n))%n=1,上式=((m^(φ(n)))^h)*m=m
2.m、n不是互质关系,n=p*q, p,q为质数，所以m=kp或m=kq
因为m<n，m=kp，k<q，kp与q互质
(kp)^(q-1) % q=1
((kp)^(q-1))^(h*(p-1))*kp % q =kp
(kp)^(e*d) % q=kp
(kp)^(e*d) %p=0，易证(kp)^(e*d) % (p*q) =kp =>m^(e*d)%n=m
2这个证明真不容易想到




这再顺带说下模乘法运算公式，高次幂暴力取模算法也会用到这个公式，证明：
设 x1%c=r1 ,x2%c=r2 ,x3%c=r2 ,..., xn%c=rn
x1=(r1+k1*c),x2=(r2+k2*c),x3=(r3+k3*c),... ,xn=(rn+kn*c);
x1*x2*x3*x4*...*xn % c = (r1+k1*c)(r2+k2*c)(r3+k3*c)...(rn+kn*c) % c=r1*r2*r3*...*rn + ?*c - ?*c =r1*r2*r3*...*rn=(x1%c)(x2%c)(x3%c)..(xn%c) % c
那么高次幂取模可以转换为 a%c*a%c*a%c... 这么一个循环的过程，这样可以在一定程度上减少溢出的可能
但这种算法时间复杂度会随着幂的增加而增加


蒙哥马利高次幂取模算法：
int get_mod(int a, int b, int c)
{
    long long res = 1;//声明为long long类型防止数据溢出
    int temp = a;
    while(b > 0)
    {
        if( b & 1)//取幂指数二进制最后一位
        {
            res = (res * temp) % c;
        }
        temp = (temp * temp) % c;
        b >>= 1;//幂指数二进制向右移动一位
    }
    return static_cast<int> (res);//以int类型返回最终结果
}
没想到蒙哥马利这算法还真和我最开始想到一块去了，就是尽可能将乘积合到一块取一次模，然而数学家是刁刁的
任何一个整数都可以写成 2^k 相加，同底幂函数相乘等于同底指数相加，那么：
? 是 0或1 假设都是1 ，代码里面是可以判断的，我这不好写
c^k % n=c^(?*2^m1+?*2^m2+?*2^m3+...?*2^0) % n= (c^(2^m1) %n *c^(2^m2) %n*....*c^(2^0) %n)%n

那么从c^(2^0)右往左开始迭代即可