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

• 本章综述
• 代码分析
• 贝祖公式

本章综述

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

代码分析

（1）判断输入的I/O是否合法：

static NTL_CHEAP_THREAD_LOCAL long iodigits = 0;
static NTL_CHEAP_THREAD_LOCAL long ioradix = 0;
  static void InitZZIO()
{
   long x;

   x = (NTL_WSP_BOUND-1)/10;          
   iodigits = 0;
   ioradix = 1;

   while (x) {
      x = x / 10;
      iodigits++;
      ioradix = ioradix * 10;
   }

   if (iodigits <= 0) TerminalError("problem with I/O");
}

注释：
这里的NTL_WSP_BOUND的值为(1L<<30)。因为长整型是一个有符号数，其二进制数的最高位是是符号位，所以一个长整型数最多表示范围是[-2147483647, 2147483647]，所以最多有30位二进制数代表数字大小。1L代表1的数据类型是long类型的，也就是长整型。“<<30”代表将二进制的1左移30位。原来1在第一位，左移30位后，1在第31位。此时代表的数字是2147483648，超出了长整型表示范围，所以需要减一。此时代表的数字是2147483647，代表最大的长整型数的大小和位数。

Iodigits 变量记录最大长整型数的十进制位数（9位）
Ioradix 变量记录最大长整型数的十进制数+1（10位，1000000000）
Iodigits，Ioradix均为静态变量。

（2）接收大整数

istream& operator>>(istream& s, ZZ& x)
{
   long c; long cval; long sign; long ndigits; long acc;
   NTL_ZZRegister(a);
   if (!s) NTL_INPUT_ERROR(s, "bad ZZ input");        //若无输入，则直接报错

   if (!iodigits) InitZZIO();     //一个长整型最大接收的十进制位数是9位，该数据存在
    a = 0;                        // iodigits中。故先要将iodigits中数据进行初始化

   SkipWhiteSpace(s);            //跳过输入时候有效数据前的空格
   c = s.peek();                  //peek()探测指针，提前检测下一位，并不改变指针指向
   if (c == '-') {
      sign = -1;
      s.get();
      c = s.peek();
   }
   else
      sign = 1;                   //确认该大整数的正负号
   cval = CharToIntVal(c);         //接受大整数的时候只接受十进制，数字在0-9之间才合法
   if (cval < 0 || cval > 9) NTL_INPUT_ERROR(s, "bad ZZ input");
   ndigits = 0;
   acc = 0;
   while (cval >= 0 && cval <= 9) {
      acc = acc*10 + cval;        
      ndigits++;
      if (ndigits == iodigits) {                //用于九位九位的存储大整数，大整数超过九位会采用截断存储处理方法
         mul(a, a, ioradix);
         add(a, a, acc);
         ndigits = 0;
         acc = 0;
      }
      s.get();                              //从输入流中提取已经存好的一位
      c = s.peek();                         //探测下一位
      cval = CharToIntVal(c);               //检测下一位输入是否合法
   }
   if (ndigits != 0) {
      long mpy = 1;
      while (ndigits > 0) {
         mpy = mpy * 10;
         ndigits--;
      }
      mul(a, a, mpy);
      add(a, a, acc);
   }
   if (sign == -1)
      negate(a, a);
   x = a;
   return s;
}

注释：
NTL_ZZRegister(a)这个函数相当于一种注册，生成了一个叫a的大整数类对象。和我们平时写的代码ZZ a;功能类似。
关于截断存储大整数的方法：NTL中定义了一个无符号长整型的数组（*cc），以十进制九位九位的将大整数进行拆分然后存入数组中。方便以后取用。
2.2.2输出大整数
（1）定义_ZZ_local_stack堆栈（采用结构体的方式定义）

	struct _ZZ_local_stack {
   long top;
   Vec<long> data;

   _ZZ_local_stack() { top = -1; }           //构造函数

   long pop() { return data[top--]; }       //弹栈，出栈
   long empty() { return (top == -1); }     //堆栈判空
   void push(long x);                       //压栈
};

void _ZZ_local_stack::push(long x)          
{
   if (top+1 >= data.length())              //判断堆栈的是否已满。若满，则扩容
      data.SetLength(max(32, long(1.414*data.length())));

   top++;
   data[top] = x;
}

注释：
堆栈的作用。因为大整数在运算的时候是从低位开始运算的，但是输出的时候是从高为开始逐渐输出的。所以低位做好运算之后适合进行压栈处理，不断地压栈，直到大整数运算结束再从堆栈中一一出栈，得到正确的输出结果。
（2）定义长整型输出函数

	static
void PrintDigits(ostream& s, long d, long justify)
{
   NTL_TLS_LOCAL_INIT(Vec<char>, buf, (INIT_SIZE, iodigits));   //初始化数组buf
   long i = 0;
   while (d) {
      buf[i] = IntValToChar(d % 10);         //将要输出的长整型数的每一位存到buf数组中
      d = d / 10;
      i++;
   }
   if (justify) {
      long j = iodigits - i;
      while (j > 0) {
         s << "0";
         j--;
      }
   }
   while (i > 0) {
      i--;
      s << buf[i];             //逐位输出长整型数据（只是*cc数组中的一个元素）
   }
}

(3)对于输出流运算符的重载

ostream& operator<<(ostream& s, const ZZ& a)
{
   ZZ b;
   _ZZ_local_stack S;
   long r;
   long k;

   if (!iodigits) InitZZIO();

   b = a;

   k = sign(b);

   if (k == 0) {
      s << "0";                         //判断输出符号，若数值为零，则直接输出“0”
      return s;
   }

   if (k < 0) {
      s << "-";                      //判断输出符号，若为负号，则先输出“-”
      negate(b, b);
   }

   do {
      r = DivRem(b, b, ioradix);
      S.push(r);
   } while (!IsZero(b));

   r = S.pop();
   PrintDigits(s, r, 0);

   while (!S.empty()) {
      r = S.pop();
      PrintDigits(s, r, 1);
   }
      
   return s;
}

贝祖公式

void XGCD(long& d, long& s, long& t, long a, long b)
{
   long  u, v, u0, v0, u1, v1, u2, v2, q, r;

   long aneg = 0, bneg = 0;

   if (a < 0) {
      if (a < -NTL_MAX_LONG) ResourceError("XGCD: integer overflow");
      a = -a;
      aneg = 1;
   }

   if (b < 0) {
      if (b < -NTL_MAX_LONG) ResourceError("XGCD: integer overflow");
      b = -b;
      bneg = 1;
   }

   u1=1; v1=0;
   u2=0; v2=1;
   u = a; v = b;

   while (v != 0) {
      q = u / v;
      r = u % v;
      u = v;
      v = r;                 //到此为止是广义欧几里得算法
      u0 = u2;               //从此开始是贝祖公式的求解
      v0 = v2;
      u2 =  u1 - q*u2;
      v2 = v1- q*v2;
      u1 = u0;
      v1 = v0;
   }

   if (aneg)
      u1 = -u1;

   if (bneg)
      v1 = -v1;

   d = u;
   s = u1;
   t = v1;                  //最后返回u = sa+tb
}

贝祖公式的数学证明详见@回首，阑珊