NTL密码算法开源库-大整数ZZ类（三）

2021SC@SDUSC

中国剩余定理

//中国剩余定理模板代码
typedef long long ll;
ll china(ll a[],ll b[],int n)//a[]为除数，b[]为余数
{
    ll M=1,y,x=0;
    for(int i=0;i<n;++i)  //算出它们累乘的结果
        M*=a[i];
    for(int i=0;i<n;++i)
    {
        ll w=M/a[i];
        ll tx=0;
        int t=exgcd(w,a[i],tx,y);  //计算逆元
        x=(x+w*(b[i]/t)*x)%M; 
    }
    return (x+M)%M;
}

具体数学原理证明见@回首~阑珊

一次同余式

设m是正整数,a是满足ma的整数.则一次同余式
ax≡1(mod m)
有解的充分必要条件是(a,m)=1.而且,当同余式(3)有解时,解是惟一的
证充分性.(存在性)因为(a,m)=1,根据广义欧几里得除法
可找到整数s,t使得
s · a + t · m = (a,m) = 1.
因此,x=s(modm)是同余式(3)的解(惟一性)若还有解x1,即ax1=1(modm),则有
a(x – x1 ) ≡ 0 (mod m)
因为(a,m)=1,所以x1≡x(modm).解是惟一的
再证必要性：若同余式(3)有解x≡xo(modm),则存在整数q,使得a·x0=1+q·m.根据定理有(a,m)=1.定理成立.证毕

乘法逆元

定义3.1.2设m是一个正整数,a是一个整数.如果存在整数a′使得
a · a 0 ≡ a 0 · a ≡ 1 (mod m)
成立,则a叫做模m可逆元
根据定理，在模m的意义下,a′是惟一存在的.这时a叫做a
的模m逆元,记作
a 0 = a −1 (mod m).
因此,在定理的条件下,同余式(3)即
ax≡1(modm)的解可写成x ≡ a ^(−1) (mod m).

简化剩余的等价描述

设m是一个正整数.则整数a是模m简化剩余的充要条件是整数α是模m逆元
证：必要性.如果整数α是模m简化剩余,则(a,m)=1.根据定
理，存在整数a′使得
a · a 0 ≡ a 0 · a ≡ 1 (mod m).
因此，a是模m逆元
充分性.如果α是模m逆元,则存在整数α′使得
a · a 0 ≡ 1 (mod m).
即同余式
ax≡1(modm)
有解x≡a(modm).又有(a,m)=1.因此,整数a是
模m简化剩余.

最后,考虑通常的一次同余式的求解
设m是一个正整数,a是满足m不整除a的整数，则一次同余式
ax ≡ b (mod m)
有解的充分必要条件是(a,m)|b.而且,当同余式有解时,其解

二次同余式

求解形如 x^2 = a (mod p) 这样的同余式
代码：

/*
模P平方根：
    求 X ^2 = a (mod p)
    定理：当P为奇素数的时候
    先判断（a / p )的勒让德符号， 若为-1则无解，若为1则有解
    分解P-1,然后求B，然后求出X(t-1),和a的逆元
    然后开始求 ans = （ a的逆元 * 上一个X的平方（t-k））的(t-k-1)次方对P取模
    如果  ans == -1 则 J = 1；
    如果  ans ==  1 则 J = 0；
    然后开始求现在的 X = （上一个X * B的（J*2的k次方）次方
    直到求出X0，也就是最后的解
*/
#include<bits/stdc++.h>
using namespace std;
 
void Divide(int p, int &t, int &s)
{
    t = 0, s= 0;
    while(p % 2 == 0)
    {
        t++;
        p /= 2;
    }
    s = p;
    return ;
}
 
int Pow(int a, int b, int mod)
{
    int ans = 1, base = a;
    while(b != 0)
    {
        if(b & 1)
            ans = (ans * base) % mod;
        base = (base * base) % mod;
        b >>= 1;
    }
    return ans;
}
 
int Legendre(int a, int p)
{
    if(a == 2)
    {
        int x = (p+1)*(p-1)/8;
        if(x % 2 == 0)
            return 1;
        else
            return -1;
    }
    else
    {
        int ans = Pow(a, (p-1)/2, p);
        if(ans == p-1)
            return -1;
        else
            return 1;
    }
}
 
int FindN(int p)
{
    for(int i = 1; i < p; i++)
    {
        if(Legendre(i, p) == -1)
            return i;
    }
}
 
int e_gcd(int a, int b, int &x, int &y)
{
    if(b == 0)
    {
        x = 1; y = 0;
        return a;
    }
    int q = e_gcd(b, a%b, y, x);
    y -= a/b*x;
    return q;
}
 
int Inverse(int a, int m)
{
    int x, y;
    int gcd = e_gcd(a, m, x, y);
    if(1 % gcd != 0)
        return -1;
    x *= 1/gcd;
    m = abs(m);
    int ans = x % m;
    if(ans <= 0)
        ans += m;
    return ans;
}
 
int JudgeJ(int A, int x, int t, int p)
{
    cout << A << " " << x << " " << t << " " << p << endl;
    x = ((x % p) * (x % p) % p);
    int xx = (A * x) % p;
    int ans = Pow(xx, pow(2, t), p);
    if(ans == p-1)
        return 1;
    else
        return 0;
}
 
int main()
{
    int a, p;
    printf("请输入所求算式的 a 和 p：\n");
    while( scanf("%d %d", &a, &p) != EOF)
    {
        if(Legendre(a, p) == -1)
        {
            cout << "无解" << endl;
            continue;
        }
        int t, s;
        Divide(p-1, t, s); //求t和s
        int n = FindN(p);  //找到那个不符合条件的n
        int b = Pow(n, s, p);
        int *X;
        X = (int *) malloc(sizeof(int) * (t+5));
        X[t-1] = Pow(a, (s+1)/2, p);
        t--;
        int A = Inverse(a, p);          //求A的逆元
        int k = 0;
        while(t > 0)
        {
            int j = JudgeJ(A, X[t], t-1, p);
            int B = Pow(b, j * pow(2, k), p);
            X[t-1] = ((X[t] % p) * (B % p)) % p;
            t--; k++;
        }
        printf("所求的解为：");
        cout << X[0] << endl;
    }
    return 0;
}

雅可比符号

long Jacobi(const ZZ& aa, const ZZ& nn)
{
   ZZ a, n;
   long t, k;
   long d;

   a = aa;
   n = nn;
   t = 1;

   while (a != 0) {
      k = MakeOdd(a);
      d = trunc_long(n, 3);
      if ((k & 1) && (d == 3 || d == 5)) t = -t;

      if (trunc_long(a, 2) == 3 && (d & 3) == 3) t = -t;
      swap(a, n);
      rem(a, a, n);
   }

   if (n == 1)
      return t;
   else
      return 0;
}