2021SC@SDUSC

本章综述

  大整数ZZ类主要实现了任意长度大整数表示、最大公因数、Jacobi符号和素性检验。笔者将通过逐个分析ZZ.cpp源代码中函数的形式来一步步向读者展示NTL是如何实现上述功能的。

计算最大公因数（gcd）

（1）数学基础：（广义）欧几里得除法

知识储备（定理，公立，公式）

·如果 b|a ,则(a,b) = b;

·如果a,b为两整数，则(a,b) = (b,a)

·如果 p为素数，a为整数，且p ∤ a，则a和p互素

   证明：设（a,p） = d ,则有d|p，且d|a 。因为p是素数，所以d = 1或者d = p

        对于d = p ,和p ∤ a矛盾，所以d = 1，即，（a，p） = 1.结论成立

·设b为任意正整数，（0，b） = b

·设a,b,c≠0且为整数。若c|a,c|b,则c|(sa+tb)。（若c|a,c|b，则c整除a,b的任意线性组合）

  证明：s*a+t*b  = s*θ*c+t*β*c = c*(s* θ+t* β)

·设a,b,c为三个不全为零的整数，如果a = q*b+c，其中q为整数，则（a,b） = (b,c)

  证明：d = (a,b) ，=> d|a，d|b  =>d|(a+(-q)b)  =>d|c  ,d|b   =>d为b,c的公因数 =>d≤d1

        d1 = (b,c),同理得：d1≤d。=>d1=d  =>(a,b) = (b,c)

广义欧几里得除法：

欧几里得算法又称辗转相除法，是指用于计算两个非负整数a，b的最大公约数。应用领域有数学和计算机两个方面。计算公式gcd(a,b) = gcd(b,a mod b)。欧几里得算法是用来求两个正整数最大公约数的算法。古希腊数学家欧几里得在其著作《The Elements》中最早描述了这种算法,所以被命名为欧几里得算法。扩展欧几里得算法可用于RSA加密等领域。

假如需要求 1997 和 615 两个正整数的最大公约数,用欧几里得算法，是这样进行的：

1997 / 615 = 3 (余 152)

615 / 152 = 4(余7)

152 / 7 = 21(余5)

7 / 5 = 1 (余2)

5 / 2 = 2 (余1)

2 / 1 = 2 (余0)

至此，最大公约数为1

以除数和余数反复做除法运算，当余数为 0 时，取当前算式除数为最大公约数，所以就得出了 1997 和 615 的最大公约数 1。

（2）代码分析：

广义欧几里得除法求最大公因数

long GCD(long a, long b)

{

   long u, v, t, x;



   if (a < 0) {

      if (a < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");

      a = -a;                    //判断输入的长整型是否溢出

   }



   if (b < 0) {

      if (b < -NTL_MAX_LONG) ResourceError("GCD: integer overflow");

      b = -b;                   //判断输入的长整型是否溢出



   }





   if (b==0)

      x = a;

   else {

      u = a;

      v = b;

      do {

         t = u % v;       //a = qb+r欧几里得核心算法（详见上知识储备-广义欧几里得除法）

         u = v;

         v = t;

      } while (v != 0);



      x = u;

   }



   return x;

}

贝祖公式（广义欧几里得的逆方法）

利用广义欧几里得的算法步骤一步步回溯，就可以找到这样一组（x,y）（详细代码见@元解~殇怀）

·贝祖公式：

求s,t的一种方法也是简单方法：是利用广义欧几里得除法先得到（a,b) 即a,b的最大公因数，然后回代得到整数s,t。

2.2.4关于模运算

数学基础：乘法逆元求解

如下例：

乘法逆元：已知A和N互素，则存在一个整数B，使得A*B=1(mod N)

利用广义欧几里得除法（辗转相除法）对余数N进行辗转相除，然后将商逆序排列（如上图第一行黑色字），然后利用贝祖公式的变形求出上图第三行的相关数据（其中第三行第一的数字永远为1，第二个数字和第一行的第一个数字一样），如：5 = 1+1*4 ；9 = 4+5*1 ；14 = 5+9*1 ；23 = 9+14*1 ；37 = 14+23*1 ……

最后求出550即为550关于模1769的乘法逆元，即（550*550）mod 1769 = 1。

long InvModStatus(long& x, long a, long n)  //求a关于n的乘法逆元，并将其赋给X。同时返回a和ns是否互素。

{

   long d, s, t;



   XGCD(d, s, t, a, n);

   if (d != 1) {

      x = d;

      return 1;

   }

   else {

      if (s < 0)

         x = s + n;                      //规定乘法逆元要是正的，如果求出来不是正数，则要加上一倍的n转成正数输出

      else

         x = s;



      return 0;

   }

}

long InvMod(long a, long n)

{

   long d, s, t;



   XGCD(d, s, t, a, n);

   if (d != 1) {

      InvModError("InvMod: inverse undefined");

   }

   if (s < 0)

      return s + n;

   else

      return s;

}

long PowerMod(long a, long ee, long n)

{

   long x, y;



   unsigned long e;



   if (ee < 0)

      e = - ((unsigned long) ee);

   else

      e = ee;



   x = 1;

   y = a;

   while (e) {

      if (e & 1) x = MulMod(x, y, n);   //算法加速，一次循环，两次计算乘积，通过右移指令和按位取与指令做到快速的循环迭代，减少循环次数，加快运算速度。

                                       //利用公式：((a{x} mod n)*(a{y} mod n))mod n = a{x+y} mod n  ---证明见后

      y = MulMod(y, y, n);

      e = e >> 1;

   }



   if (ee < 0) x = InvMod(x, n);



   return x;

}

证明：((a{x} mod n)*(a{y} mod n))mod n = a{x+y} mod n

原式左侧 = ((a{x}-k1*n)*(a{y}-k2*n))mod n

         = (a{x+y}-n(k2* a{x}+k1* a{y})+n{2}*k1*k2)mod n

        = (a{x+y}) mod n