开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
_integer.c
资源名称：leda.tar.gz　[点击查看]
上传用户：gzelex
上传日期：2007-01-07
资源大小：707k
文件大小：39k
源码类别：

数值算法/人工智能

开发平台：

MultiPlatform

_integer.c：源码内容
  1. /*******************************************************************************
  2. +
  3. +  LEDA-R  3.2.3
  4. +
  5. +  _integer.c
  6. +
  7. +  Copyright (c) 1995  by  Max-Planck-Institut fuer Informatik
  8. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  9. +  All rights reserved.
  10. + 
  11. *******************************************************************************/
  14. //------------------------------------------------------------------------------
  15. //
  16. //  integer:  big integers
  17. //
  18. //  by Christian Uhrig and Stefan Naeher   (1994/1995)
  19. //
  20. //  - 32 bit (unsigned long) vector representation
  21. //  - use of a handle class concept (to avoid copy overhead)
  22. //  - use of the LEDA memory management
  23. //  - sparc assembler code for 32 bit multiplication and division
  24. //  - Karatsuba multiplication for long integers
  25. //  - loop unrolling
  26. //
  27. //  - sparc assembler for add and sub (S. Naeher, 1995)
  28. //
  29. //------------------------------------------------------------------------------
  32. #include <LEDA/integer.h>
  33. #include <math.h>
  34. #include <ctype.h>
  37. typedef int_word_type word;
  39. typedef unsigned int sz_t;
  42. #define WORD_LENGTH  32 
  43. #define HALF_LENGTH  16 
  44. #define KARA_LIMIT    8 
  46. #define MAX_WORD     0xFFFFFFFF
  47. #define RIGHT_MASK   0x0000FFFF
  48. #define LEFT_MASK    0xFFFF0000
  50. #define LOW_WORD(w)  (w & RIGHT_MASK)
  51. #define HIGH_WORD(w) (w >> HALF_LENGTH)
  53. enum {ZERO=0, NEGATIVE=-1, POSITIVE=1 };
  55. const double pow32 = ldexp(1,WORD_LENGTH);
  59. inline int next_power(int s)
  60. { int p = 1;
  61.   while (p < s) p <<= 1;
  62.   return p;
  63. }
  66. void delete_integer_rep(integer_rep* p)
  67. { int s = sizeof(integer_rep) + (p->size-1)*sizeof(unsigned long);
  68.   if (s < 256)
  69.      deallocate_bytes(p,s);
  70.   else
  71.      delete[] p;
  72. }
  75. integer::integer()
  76. { PTR = new_integer_rep(1);
  77.   PTR->vec[0] = 0;
  78.   PTR->sign = ZERO;
  79.   PTR->used = 0;
  80.  }
  84. integer::integer(unsigned int n)
  85. {
  86.   PTR = new_integer_rep(1);
  88.   PTR->vec[0] = n;
  90.   if (n == 0)
  91.      { PTR->sign = ZERO;
  92.        PTR->used = 0; }
  93.   else
  94.      { PTR->sign = POSITIVE;
  95.        PTR->used = 1; }
  96. }
  100. integer::integer(int n)
  101. {
  102.   PTR = new_integer_rep(1);
  104.   PTR->used = 1;
  105.   PTR->vec[0] = n;
  106.   PTR->sign = POSITIVE;
  108.   if (n == 0)
  109.   { PTR->sign = ZERO;
  110.     PTR->used = 0;
  111.    }
  113.   if (n < 0) 
  114.   { PTR->sign = NEGATIVE;
  115.     PTR->vec[0] = -n;
  116.    }
  117. }
  120. integer::integer(long n)
  121. {
  122.   PTR = new_integer_rep(1);
  124.   PTR->used = 1;
  125.   PTR->vec[0] = n;
  126.   PTR->sign = POSITIVE;
  128.   if (n == 0)
  129.   { PTR->sign = ZERO;
  130.     PTR->used = 0;
  131.    }
  133.   if (n < 0) 
  134.   { PTR->sign = NEGATIVE;
  135.     PTR->vec[0] = -n;
  136.    }
  137. }
  140. integer::integer(unsigned long n)
  141. {
  142.   PTR = new_integer_rep(1);
  144.   if (n == 0)
  145.   { PTR->sign = ZERO;
  146.     PTR->used = 0;
  147.    }
  148.   else
  149.   { PTR->sign = POSITIVE;
  150.     PTR->used = 1;
  151.     PTR->vec[0] = n;
  152.    }
  154. }
  158. integer::integer(double x) 
  159. { 
  160.   int sig = POSITIVE;
  162.   if (x < 0) 
  163.   { sig = NEGATIVE;
  164.     x = -x;
  165.    }
  167.   if (x < 1)
  168.   { PTR = new_integer_rep(1);
  169.     PTR->used = 0;
  170.     PTR->sign = ZERO;
  171.     return;
  172.    }
  174.   if (x < double(MAX_WORD))
  175.   { PTR = new_integer_rep(1);
  176.     PTR->vec[0] = word(x);
  177.     PTR->used = 1;
  178.     PTR->sign = sig;
  179.     return;
  180.    }
  183.   int expt;
  184.   double mantissa = frexp(x, &expt);
  185.   double int_mant;
  187.   int r = expt % WORD_LENGTH;
  188.   int q = expt / WORD_LENGTH;
  190.   int used = q;
  191.   if (r) used++;
  193.   int sz = next_power(used);
  195.   PTR = new_integer_rep(sz);
  196.   PTR->used = used;
  197.   PTR->sign = sig;
  199.   word* p = PTR->vec+used-1;
  201.   if (r)
  202.   { mantissa *= ldexp(1,r);
  203.     mantissa = modf(mantissa, &int_mant);
  204.     *p--= (unsigned long)int_mant; // first r bits
  205.    }
  207.   while (q--)
  208.     if (mantissa != 0)
  209.     { mantissa *= pow32;
  210.       mantissa = modf(mantissa, &int_mant);
  211.       *p-- = (unsigned long)int_mant; // next 32 bits
  212.      }
  213.     else *p-- = 0;
  215.   if (used == 0) PTR->sign = ZERO;
  217. }
  221. //------------------------------------------------------------------------------
  222. // local functions
  223. //------------------------------------------------------------------------------
  225. inline word DivTwoWordsByOne(word high, word low, word D, word* q)
  226. { 
  227.   // precondition:  high < D  
  228.   // *q = [high,low] / D, returns remainder
  229.   // the result has at most 32 bits  
  230.   // we can compute it with double precision floating point division
  232.   double H = high;
  233.   double L = low;
  234.   word   Q = word((H/D)*pow32 + L/D);
  236.   D *= Q;
  238.   *q = Q;
  239.   return low-D;
  240. }
  243. static int absolute_cmp(word* a, int a_used, word* b, int b_used)
  244. { 
  245.   if (a_used > b_used) return  1;
  246.   if (a_used < b_used) return -1;
  247.   word* ap = a + a_used;
  248.   word* bp = b + b_used;
  249.   while(ap > a)
  250.   { word A = *--ap;
  251.     word B = *--bp;
  252.     if (A > B) return  1;
  253.     if (A < B) return -1;
  254.    }
  255.   return 0;
  256. }
  259. static void Init_Vector(word* a, sz_t a_used)
  260. { 
  261.   switch (a_used % 16)
  262.   { case 15:*a++ = 0;
  263.     case 14:*a++ = 0;
  264.     case 13:*a++ = 0;
  265.     case 12:*a++ = 0;
  266.     case 11:*a++ = 0;
  267.     case 10:*a++ = 0;
  268.     case  9:*a++ = 0;
  269.     case  8:*a++ = 0;
  270.     case  7:*a++ = 0;
  271.     case  6:*a++ = 0;
  272.     case  5:*a++ = 0;
  273.     case  4:*a++ = 0;
  274.     case  3:*a++ = 0;
  275.     case  2:*a++ = 0;
  276.     case  1:*a++ = 0;
  277.    }
  279.   int n = a_used / 16;
  281.   while (n--)
  282.     { *a++ = 0;  *a++ = 0;  *a++ = 0;  *a++ = 0;
  283.       *a++ = 0;  *a++ = 0;  *a++ = 0;  *a++ = 0;
  284.       *a++ = 0;  *a++ = 0;  *a++ = 0;  *a++ = 0;
  285.       *a++ = 0;  *a++ = 0;  *a++ = 0;  *a++ = 0; }
  286. }
  290. static word* Copy_Vector(word* a, word* b, sz_t b_used)
  291. { 
  292.   // copy b[0..b_used-1] to a[]
  294.   switch (b_used % 16)
  295.   { case 15:*a++ = *b++;
  296.     case 14:*a++ = *b++;
  297.     case 13:*a++ = *b++;
  298.     case 12:*a++ = *b++;
  299.     case 11:*a++ = *b++;
  300.     case 10:*a++ = *b++;
  301.     case  9:*a++ = *b++;
  302.     case  8:*a++ = *b++;
  303.     case  7:*a++ = *b++;
  304.     case  6:*a++ = *b++;
  305.     case  5:*a++ = *b++;
  306.     case  4:*a++ = *b++;
  307.     case  3:*a++ = *b++;
  308.     case  2:*a++ = *b++;
  309.     case  1:*a++ = *b++; }
  311.   int n = b_used / 16;
  313.   while (n--)
  314.     { *a++ = *b++;  *a++ = *b++;  *a++ = *b++;  *a++ = *b++;
  315.       *a++ = *b++;  *a++ = *b++;  *a++ = *b++;  *a++ = *b++;
  316.       *a++ = *b++;  *a++ = *b++;  *a++ = *b++;  *a++ = *b++;
  317.       *a++ = *b++;  *a++ = *b++;  *a++ = *b++;  *a++ = *b++; }
  319.   return a;
  320. }
  324. static void Print_Vector(word* a, int a_used, ostream& out) 
  325. {
  326.   char* digits   = new char[10*a_used];
  327.   word* tmp      = new word[a_used];
  328.   word* tmp_high = Copy_Vector(tmp,a,a_used) - 1;
  330.   int i = 0;
  331.   while (tmp_high >= tmp)
  332.   { word r = 0;
  333.     for (word* p = tmp_high; p >= tmp; p--) 
  334.        r = DivTwoWordsByOne(r,*p,10,p);
  335.     if (*tmp_high == 0) tmp_high--;
  336.     digits[i++] = '0' + char(r);
  337.    }
  339.   while (i--) out << digits[i];
  341.   delete[] digits;
  342.   delete[] tmp;
  343. }
  347. #if defined(sparc)
  349. //------------------------------------------------------------------------------
  350. // sparc assembler code
  351. //------------------------------------------------------------------------------
  352. // _sparc_mult.s: Multiply_Words & Mult_Inner_Loop
  353. // _sparc_div.s : Div_Inner_Loop
  354. // _sparc_add.s : School_Add
  355. // _sparc_sub.s : School_Sub
  356. //------------------------------------------------------------------------------
  358. #if defined(__svr4__)
  359. #define Multiply_Words  _Multiply_Words
  360. #define Mult_Inner_Loop _Mult_Inner_Loop
  361. #define Div_Inner_Loop  _Div_Inner_Loop
  362. #define School_Add      _School_Add
  363. #define School_Sub      _School_Sub
  364. #endif
  366. extern "C" sz_t School_Add(word*,sz_t,word*,sz_t,word*);
  367. extern "C" sz_t School_Sub(word*,sz_t,word*,sz_t,word*);
  368. extern "C" word Multiply_Words(word a, word b, word* high);
  369. extern "C" word Mult_Inner_Loop(word* p, word* a, word* a_stop, word B);
  370. extern "C" word Div_Inner_Loop(word* p, word* a, word* a_stop, word B);
  372. #else
  373. /* if not sparc */
  375. inline word Multiply_Words(word a, word b, word* high)
  376. { word al = LOW_WORD(a);
  377.   word ah = HIGH_WORD(a);
  378.   word bl = LOW_WORD(b);
  379.   word bh = HIGH_WORD(b);
  380.   word c,L,H;
  381.   c = bl*al;           
  382.   L = LOW_WORD(c);
  383.   c = HIGH_WORD(c) + bl*ah;
  384.   H = HIGH_WORD(c);
  385.   c = LOW_WORD(c)  + bh*al;
  386.   L |= (LOW_WORD(c) << HALF_LENGTH);
  387.   *high =  H + HIGH_WORD(c) + bh*ah;
  388.   return L;
  389. }
  392. inline word Mult_Inner_Loop(word* p, word* a, word* a_stop, word B)
  393. {
  394.    // p[] += a[]*B, return carry
  396.    word low,high;
  397.    word carry = 0;
  398.    while (a < a_stop)
  399.    { word P = *p;
  400.      low = Multiply_Words(*a++,B,&high);
  401.      P += low;
  402.      if (P < low )  high++; // carry in addition
  403.      P += carry;
  404.      if (P < carry) high++; // carry in addition
  405.      *p++ = P;
  406.      carry = high;
  407.    }
  408.    *p = carry;
  409.    return carry;
  410. }
  413. inline word Div_Inner_Loop(word* p, word* a, word* a_stop, word B)
  414. { 
  415.    // p[] -= a[]*B, return carry
  417.    word low,high;
  418.    word carry = 0;
  419.    word P;
  420.    while (a < a_stop)
  421.    { P = *p;
  422.      low = Multiply_Words(*a++,B,&high);
  423.      if (P < low) high++;  // carry
  424.      P -= low;
  425.      if (P < carry) high++;  // carry
  426.      P -= carry;
  427.      *p++ = P;
  428.      carry = high;
  429.    }
  430.    P = *p;
  431.    *p = P - carry;
  432.   
  433.    return (P < carry);
  434. }
  438. static sz_t School_Add(word *a, sz_t a_used, word *b, sz_t b_used, word* sum)
  439. {
  440.   // compute sum = a + b   (a_used >= b_used)
  442.   int carry = 0;
  443.   word aa;
  444.   word bb;
  445.   word* s_last = sum + a_used;
  447. #define ADD_LOOP_BODY {
  448.   aa = *a++;
  449.   bb = *b++;
  450.   aa += bb;
  451.   if (carry) { if (++aa > bb) carry = false; }
  452.   else if (aa < bb) carry = true;
  453.   *sum++ = aa; }
  455.   int r = b_used % 16;
  457.   switch (r) { 
  458.       case 15: ADD_LOOP_BODY;
  459.       case 14: ADD_LOOP_BODY;
  460.       case 13: ADD_LOOP_BODY;
  461.       case 12: ADD_LOOP_BODY;
  462.       case 11: ADD_LOOP_BODY;
  463.       case 10: ADD_LOOP_BODY;
  464.       case  9: ADD_LOOP_BODY;
  465.       case  8: ADD_LOOP_BODY;
  466.       case  7: ADD_LOOP_BODY;
  467.       case  6: ADD_LOOP_BODY;
  468.       case  5: ADD_LOOP_BODY;
  469.       case  4: ADD_LOOP_BODY;
  470.       case  3: ADD_LOOP_BODY;
  471.       case  2: ADD_LOOP_BODY;
  472.       case  1: ADD_LOOP_BODY;
  473.      }
  475.   int n = b_used / 16;
  477.   while (n--)
  478.   { ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; 
  479.     ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; 
  480.     ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; 
  481.     ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY; ADD_LOOP_BODY;
  482.    }
  484.   if (sum != a) // copy rest of a to sum
  485.   { Copy_Vector(sum,a,a_used-b_used);
  486.     *s_last = 0;
  487.    }
  489.   if (carry) // propagate carry
  490.   { while (++*sum == 0) sum++;
  491.     if (sum == s_last) a_used++;
  492.    }
  494.   return a_used;
  495. }
  498. static sz_t School_Sub(word *a, sz_t a_used, word *b, sz_t b_used, word* diff)
  499. {
  500.   // compute diff = a - b    (a > b)
  502. #define SUB_LOOP_BODY {
  503.   aa = *a++;
  504.   bb = *b++;
  505.   if (borrow && ++bb) borrow = false;
  506.   if (aa < bb)  borrow = true;
  507.   *diff++ = aa - bb; }
  509.   word  aa;
  510.   word  bb;
  511.   bool  borrow = false;
  513.   word* d_stop = diff + a_used;
  515.   switch (b_used % 16) 
  516.    { case 15: SUB_LOOP_BODY;
  517.      case 14: SUB_LOOP_BODY;
  518.      case 13: SUB_LOOP_BODY;
  519.      case 12: SUB_LOOP_BODY;
  520.      case 11: SUB_LOOP_BODY;
  521.      case 10: SUB_LOOP_BODY;
  522.      case  9: SUB_LOOP_BODY;
  523.      case  8: SUB_LOOP_BODY;
  524.      case  7: SUB_LOOP_BODY;
  525.      case  6: SUB_LOOP_BODY;
  526.      case  5: SUB_LOOP_BODY;
  527.      case  4: SUB_LOOP_BODY;
  528.      case  3: SUB_LOOP_BODY;
  529.      case  2: SUB_LOOP_BODY;
  530.      case  1: SUB_LOOP_BODY;
  531.     }
  533.   int n = b_used / 16;
  535.   while (n--)
  536.   { SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; 
  537.     SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; 
  538.     SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; 
  539.     SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY; SUB_LOOP_BODY;
  540.    }
  542.   if (diff != a) Copy_Vector(diff,a,a_used-b_used);
  544.   if (borrow)
  545.     while (--*diff == MAX_WORD) diff++;
  547.   while (*--d_stop == 0) a_used--;
  549.   return a_used;
  550. }
  552. /* end of not-sparc section */
  553. #endif
  556. static sz_t School_Mult(word *a, sz_t a_used, word *b,int b_used, word* prod) 
  557. { 
  558.   sz_t p_used = a_used + b_used;
  559.   word* a_stop = a + a_used;
  560.   word* b_last = b + b_used - 1;
  562.   Init_Vector(prod,a_used);
  564.   while (b < b_last)
  565.       Mult_Inner_Loop(prod++,a,a_stop,*b++);
  567.   if (Mult_Inner_Loop(prod,a,a_stop,*b) == 0) p_used--;   // no carry
  569.   return p_used;
  570. }
  574. static void Karatsuba(word* prod, word* a, word* b, word* tmp, int n)
  575. {       
  577.  /*  Input: 
  578.                 4n      3n      2n       n        0
  580.                                          A1       A0
  581.                                  +-------+--------+
  582.    a[0..2n-1]                    |  a1   |  a0    |
  583.                                  +-------+--------+
  585.                                          B1       B0
  586.                                  +-------+--------+
  587.    b[0..2n-1]                    |  b1   |  b0    |
  588.                                  +-------+--------+
  590.                                  T2      T1       T0
  591.                  +-------+-------+-------+--------+
  592.    tmp[0..4n-1]  |       |       |       |        |    space for temporaries
  593.                  +-------+-------+-------+--------+
  596.   Output:                P3      P2      P1       P0
  597.                  +-------+-------+-------+--------+
  598.   prod[0..4n-1]  |           a  *  b              | 
  599.                  +-------+-------+-------+--------+
  603.  Let  
  604.          a = a1*B^n + a0 
  605.          b = b1*B^n + b0
  607.  then
  608.        a*b = (a1*b1)*B^2n + (a1*b0)*B^n + (a0*b1)*B^n + (a0*b0)
  610.            = (a0*b0) + (a1*b1)*B^2n + (a1*b1 + a0*b0 + (a1-a0)(b0-b1))*B^n 
  613.        4n      3n      2n       n        0
  614.         +-------+-------+-------+--------+
  615.         |    a1 * b1    |    a0 * b0     |
  616.         +-------+-------+-------+--------+
  617.         +-------+-------+-------+--------+
  618.       + |   0   |    a1 * b1    |   0    |
  619.         +-------+-------+-------+--------+
  620.         +-------+-------+-------+--------+
  621.       + |   0   |    a0 * b0    |   0    |
  622.         +-------+-------+-------+--------+
  623.         +-------+-------+-------+--------+
  624.       + |   0   |(a1-a0)*(b0-b1)|   0    |
  625.         +-------+-------+-------+--------+
  626.      ----------------------------------------
  627.         +--------------------------------+
  628.       = |         a     *     b          |
  629.         +--------------------------------+
  632.    We compute the three partial products a1*b1, a0*b0, and (a1-a0)*(b0-b1)
  633.    recursively as long as n is even and greater or equal to KARA_LIMIT
  634.    and by a call of School_Mult otherwise.
  636.   */
  639.   if (n & 1 || n < KARA_LIMIT)
  640.   { School_Mult(a,2*n,b,2*n,prod);
  641.     return;
  642.    }
  645.   word* A0 = a;
  646.   word* A1 = a + n;
  648.   word* B0 = b;
  649.   word* B1 = b + n;
  651.   word* P0 = prod;
  652.   word* P1 = prod + n;
  653.   word* P2 = prod + 2*n;
  654.   word* P3 = prod + 3*n;
  656.   word* T0 = tmp;
  657.   // word* T1 = tmp + n; (never used)
  658.   word* T2 = tmp + 2*n;
  661.   int a_sign = absolute_cmp(A1,n,A0,n);
  662.   int b_sign = absolute_cmp(B0,n,B1,n);
  664.   Karatsuba(P0,A0,B0,T0,n/2);   // prod[0..2n]  := a0 * b0 
  666.   if (a_sign >= 0)
  667.      School_Sub(A1,n,A0,n,P2);  // prod[2n..3n] := a1 - a0
  668.   else 
  669.      School_Sub(A0,n,A1,n,P2);  // prod[2n..3n] := a0 - a1
  671.   if (b_sign >= 0)
  672.      School_Sub(B0,n,B1,n,P3);  // prod[3n..4n] := b0 - b1
  673.   else
  674.      School_Sub(B1,n,B0,n,P3);  // prod[3n..4n] := b1 - b0
  677.   Karatsuba(T2,P2,P3,T0,n/2);   // tmp[2n..4n]  := |a1-a0|*|b0-b1|
  678.   Karatsuba(P2,A1,B1,T0,n/2);   // prod[2n..4n] := a1 * b1
  680.   Copy_Vector(T0,P0,2*n);
  682.   /* now we have 
  683.                       P3      P2      P1      P0
  684.               +-------+-------+-------+-------+
  685.      prod =   |     a1*b1     |     a0*b0     |
  686.               +-------+-------+-------+-------+
  687.               
  688.                               T2      T1      T0
  689.               +-------+-------+-------+-------+
  690.       tmp =   ||a1-a0|*|b0-b1||     a0*b0     |
  691.               +-------+-------+-------+-------+
  692.   */
  695.   // add up partial products
  697.   School_Add(P1,3*n,P2,2*n,P1);   // prod[1n..3n] += prod[2n..4n]
  698.   School_Add(P1,3*n,T0,2*n,P1);   // prod[1n..3n] +=  tmp[0n..2n]
  700.   if (a_sign == b_sign)  
  701.     School_Add(P1,3*n,T2,2*n,P1); // prod[1n..3n] += tmp[2n..4n]
  702.   else
  703.     School_Sub(P1,3*n,T2,2*n,P1); // prod[1n..3n] -= tmp[2n..4n]
  705. }
  710. static sz_t Kara_Mult(word* a, sz_t a_used, word* b, sz_t b_used, word* prod)
  711. { 
  712.   // a_used >= b_used
  714.   if (b_used < 4*KARA_LIMIT)
  715.      return School_Mult(a, a_used, b, b_used, prod);
  717.   sz_t p_used = a_used + b_used;
  719.   int m = next_power(b_used)/KARA_LIMIT;
  720.   int n = 0;
  722.   while (n < b_used) n += m; // n = i*2^k with i < KARA_LIMIT
  724.   word* ap = a;
  725.   word* bp = b;
  727.   // extend length of b to n (append zeroes) 
  728.   if (n > b_used)
  729.   { bp = new word[n];
  730.     word* p = Copy_Vector(bp,b,b_used);
  731.     Init_Vector(p,n-b_used);
  732.    }
  734.   // partition a into t intervals of size n each
  735.   // r is size of remainder
  737.   int t = a_used / n;
  738.   int r = a_used % n;
  740.   if (r > n/2) // make length(a) a muliple of n
  741.   { t++;
  742.     r = 0;
  743.     ap  = new word[t*n];
  744.     word* p = Copy_Vector(ap,a,a_used);
  745.     Init_Vector(p,t*n-a_used);
  746.     a_used = t*n;
  747.    }
  749.   word* res = new word[a_used + n]; // result
  750.   Init_Vector(res,a_used+n);
  752.   word* tpr = new word[2*n];        // temporary product
  753.   word* tmp = new word[2*n];        // temporary space used by Karatsuba
  755.   word* q = ap;   // current position in ap
  756.   word* p = res;  // current position in result
  758.   while (t--)
  759.   { Karatsuba(tpr, q, bp, tmp, n/2);  // tpr := q[0..n-1] * bp[0..n-1]
  760.     School_Add(p,2*n,tpr,2*n,p);      // add tpr to result
  761.     q += n;
  762.     p += n;
  763.    }
  764.      
  765.   // use School_Mult to multiply rest of a with bp
  766.   while (r--) Mult_Inner_Loop(p++,bp,bp+b_used, *q++);
  769.   p = Copy_Vector(prod,res,p_used);  // copy result to pred
  770.   while (*--p == 0) p_used--;        // and adjust p_used
  772.   if (ap != a) delete[] ap;
  773.   if (bp != b) delete[] bp;
  774.   delete[] tmp;
  775.   delete[] res;
  777.   return p_used; 
  778. }
  779.   
  782. static word Shift_Left(word *a, word *b, sz_t length, sz_t shift)
  783. { 
  784.   // auxiliary function used in School_Div
  785.   // a := (b << shift) and return carry   (0 < shift < 32)
  786.   
  787.   word* b_stop = b + length;
  789.   if (shift == 0)
  790.   { while (b < b_stop) *a++ = *b++;
  791.     return 0;
  792.    }
  794.   int  r = WORD_LENGTH-shift;
  795.   word carry = 0;
  796.   while (b < b_stop)
  797.   { word B = *b++;
  798.     carry |= (B << shift);
  799.     *a++ = carry;
  800.     carry = (B >> r);
  801.    }
  802.   return carry; 
  803. }
  807. static sz_t School_Div(word* A, sz_t a_used, word* B, sz_t b_used,
  808.                         word* Q, word* R, sz_t* R_used)
  809. {
  810.   word* a      = new word[a_used + 1];
  811.   word* b      = new word[b_used];
  812.   word* b_stop = b + b_used;
  814.   sz_t quot_used = a_used - b_used + 1;
  816.   // we first normalize A and B, i.e. ...
  818.   sz_t shift = 0;
  819.   word sh = B[b_used-1];
  820.   while ((sh & 0x80000000) == 0)
  821.    { sh <<= 1;
  822.      shift++;
  823.    }
  825.   a[a_used] = Shift_Left(a, A, a_used, shift);
  826.   Shift_Left(b, B, b_used, shift);
  828.   // two often used values
  830.   word b1 = b[b_used-1];               
  831.   word b2 = (b_used > 1) ? b[b_used-2] : 0;
  833.   /* division main loop */
  835.   word* a_ptr = a + a_used;        // current position in a
  836.   word* Q_ptr = Q + quot_used - 1; // current position in Q
  837.  
  838.   while (Q_ptr >= Q)
  839.   {
  840.     word q = MAX_WORD;
  842.     if (*a_ptr != b1) 
  843.     { word r = DivTwoWordsByOne(*a_ptr, *(a_ptr-1), b1, &q);
  844.       word h11;
  845.       word h10 = Multiply_Words(q, b2, &h11);
  846.       word h21 = r;
  847.       word h20 = *(a_ptr-2);
  849.       while( h11 > h21 || (h11 == h21 && h10 > h20) ) // [h11,h10] > [h21,h20]
  850.       { q--;
  851.         if (h10 < b2) h11--;
  852.         h10 -= b2;
  853.         h21 += b1;
  854.         if (h21 < b1) break;
  855.        }
  856.      }
  859.      if (Div_Inner_Loop(a_ptr-b_used,b,b_stop,q)) // a -= b*q
  860.      { // if carry add one b back to a and decrement q  (never executed ?)
  861.        word* ap;
  862.        word* bp;
  863.        bool  carry = false;
  864.        for(bp = b, ap = a_ptr-b_used; bp < b_stop; bp++, ap++) 
  865.        { word aa = *ap;
  866.          word bb = *bp;
  867.          bb += aa;
  868.          if (carry)
  869.            { if (++bb > aa) carry = false; }
  870.          else
  871.              if (bb < aa) carry = true;
  872.          *ap = bb;
  873.         }
  874.        if (carry) (*ap)++; 
  875.        q--;
  876.      }
  878.     *Q_ptr-- = q;
  879.      a_ptr--;
  880.   }
  882.   
  883.   if (R)  // copy remainder to R[]
  884.   { word* p = a + b_used;
  885.     word* q = R + b_used;
  886.     // copy and shift a to R[]
  887.     if (shift) 
  888.       { int  r = WORD_LENGTH-shift;
  889.         word carry = 0;
  890.         while (p > a)
  891.         { word A = *--p;
  892.           carry |= (A >> shift);
  893.           *--q = carry;
  894.           carry = (A << r);
  895.          }
  896.        }
  897.     else
  898.        while (p > a) *--q = *--p;
  900.     for(q = R+b_used-1; b_used && *q==0; q--)  b_used--; 
  902.     *R_used = b_used;
  903.    }
  906.   delete[] a;
  907.   delete[] b;
  909.   if (Q[quot_used-1] == 0) quot_used--;
  910.   return quot_used;
  911. }
  916. static sz_t DivByOneWord(word *a, sz_t a_used, word B, word* quot, word* rem)
  917. {
  918.   word  r = 0;
  919.   word* a_ptr = a + a_used;
  920.   word* q_ptr = quot + a_used;
  922.   while (a_ptr > a) r = DivTwoWordsByOne(r,*--a_ptr,B,--q_ptr);
  924.   if (quot[a_used-1] == 0) a_used--;
  925.   if (rem) *rem = r;
  926.   return a_used;
  927. }
  932. //------------------------------------------------------------------------------
  933. // friend operators  + - * / % == < > ...
  934. //------------------------------------------------------------------------------
  938. integer operator+(const integer & a, const integer & b) 
  939. {
  940.   sz_t a_sign = a.PTR->sign;
  941.   sz_t b_sign = b.PTR->sign;
  943.   sz_t sum_size;
  945.   if (a_sign == ZERO) return b;
  946.   if (b_sign == ZERO) return a;
  948.   integer_rep* sum_ptr;
  950.   sz_t a_used = a.PTR->used;
  951.   sz_t b_used = b.PTR->used;
  952.   sz_t a_size = a.PTR->size;
  953.   sz_t b_size = b.PTR->size;
  955.   word* av = a.PTR->vec;
  956.   word* bv = b.PTR->vec;
  958.   // we distinguish the cases of equal and unequal sign
  960.   if (a_sign == b_sign)
  961.   { if (a_used >= b_used)
  962.      { sum_size = (a_used+1 <= a_size) ? a_size : 2*a_size;
  963.        sum_ptr = new_integer_rep(sum_size);
  964.        sum_ptr->used = School_Add(av, a_used, bv, b_used, sum_ptr->vec); 
  965.       }
  966.     else
  967.      { sum_size = (b_used+1 <= b_size) ? b_size : 2*b_size;
  968.        sum_ptr = new_integer_rep(sum_size);
  969.        sum_ptr->used = School_Add(bv, b_used, av, a_used, sum_ptr->vec); 
  970.       }
  972.     sum_ptr->sign = a_sign;
  973.     return sum_ptr;
  974.    }
  977. /* if |a| and |b| have different signs, we must subtract the smaller absolute
  978.  * value from the greater one and set the sign appropriately
  979.  */
  981.   int rel = a_used - b_used;
  983.   if (rel == 0)
  984.   { word* p = av + a_used;
  985.     word* q = bv + b_used;
  986.     while (a_used && *--p == *--q) 
  987.     { a_used--;
  988.       b_used--;
  989.      }
  990.     if (a_used == 0) return 0;
  991.     rel  = (*p > *q) ? 1 : -1;
  992.    }
  994.   if (rel > 0) 
  995.     { // |a| > |b|
  996.       sum_ptr = new_integer_rep(a_size);
  997.       sum_ptr->used = School_Sub(av, a_used, bv, b_used, sum_ptr->vec);
  998.       sum_ptr->sign = a_sign;
  999.      }
  1000.   else
  1001.     { // |a| < |b|
  1002.       sum_ptr = new_integer_rep(b_size);
  1003.       sum_ptr->used = School_Sub(bv, b_used, av, a_used, sum_ptr->vec);
  1004.       sum_ptr->sign = b_sign;
  1005.      }
  1007.   return sum_ptr;
  1008. }
  1013. integer operator-(const integer & a, const integer & b)
  1014. {
  1015.   if (a.PTR == b.PTR) return 0;
  1017.   sz_t a_sign = a.PTR->sign;
  1018.   sz_t b_sign = b.PTR->sign;
  1019.   sz_t a_used = a.PTR->used;
  1020.   sz_t b_used = b.PTR->used;
  1021.   sz_t a_size = a.PTR->size;
  1022.   sz_t b_size = b.PTR->size;
  1025.   word* av = a.PTR->vec;
  1026.   word* bv = b.PTR->vec;
  1028.   sz_t diff_size;
  1030.   if (b_sign == ZERO) return a;
  1031.   if (a_sign == ZERO) return -b;
  1033.   integer_rep* diff_ptr;
  1035.   /* we distinguish the cases of equal and unequal signs */
  1037.   if (a_sign == b_sign)
  1038.   { 
  1039.     int rel = a_used - b_used;
  1041.     if (rel == 0)
  1042.     { word* p = av + a_used;
  1043.       word* q = bv + b_used;
  1044.       while (a_used && *--p == *--q) 
  1045.       { a_used--;
  1046.         b_used--;
  1047.        }
  1048.       if (a_used == 0) return 0;
  1049.       rel  = (*p > *q) ? 1 : -1;
  1050.      }
  1051.     
  1052.     if (rel > 0)
  1053.     { // |a| > |b| 
  1054.       diff_ptr = new_integer_rep(a_size);
  1055.       diff_ptr->used = School_Sub(av, a_used, bv, b_used, diff_ptr->vec);
  1056.       diff_ptr->sign = a_sign;  
  1057.       return diff_ptr;
  1058.      }
  1059.     else 
  1060.     {  // |b| > |a|
  1061.       diff_ptr = new_integer_rep(b_size);
  1062.       diff_ptr->used = School_Sub(bv, b_used, av, a_used, diff_ptr->vec);
  1063.       diff_ptr->sign = -a_sign;
  1064.       return diff_ptr;
  1065.     }
  1067.   }
  1069. /*
  1070.  * if |a| and |b| have different signs, we must carry out an addition and set
  1071.  * the sign appropriately
  1072.  */
  1074.   if (a_used >= b_used)
  1075.     { diff_size = (a_used+1 <= a_size) ? a_size : 2*a_size;
  1076.       diff_ptr  = new_integer_rep(diff_size);
  1077.       diff_ptr->used = School_Add(av, a_used, bv, b_used, diff_ptr->vec);
  1078.      }
  1079.   else
  1080.     { diff_size = (b_used+1 <= b_size) ? b_size : 2*b_size;
  1081.       diff_ptr  = new_integer_rep(diff_size);
  1082.       diff_ptr->used = School_Add(bv, b_used, av, a_used, diff_ptr->vec);
  1083.      }
  1085.   diff_ptr->sign = a_sign;
  1087.   return diff_ptr;
  1088. }
  1093. integer operator*(const integer& a, const integer& b) 
  1094. {
  1095.   sz_t a_sign = a.PTR->sign;
  1096.   sz_t b_sign = b.PTR->sign;
  1098.   if (a_sign == ZERO) return a;
  1099.   if (b_sign == ZERO) return b;
  1101.   sz_t a_used = a.PTR->used;
  1102.   sz_t b_used = b.PTR->used;
  1104.   integer_rep* prod_ptr;
  1106.   if (a_used == 1 && b_used == 1)
  1107.   { prod_ptr = new_integer_rep(2);
  1108.     word* p  = prod_ptr->vec;
  1109.     prod_ptr->sign = (a_sign == b_sign ? POSITIVE : NEGATIVE);
  1110.     *p = Multiply_Words(a.PTR->vec[0], b.PTR->vec[0], p+1);
  1111.     prod_ptr->used = (p[1]) ? 2 : 1;
  1112.     return prod_ptr;
  1113.    }
  1115.   int sz = next_power(a_used+b_used);
  1117.   prod_ptr = new_integer_rep(sz);
  1118.   prod_ptr->sign = (a_sign == b_sign ? POSITIVE : NEGATIVE);
  1121.   if (a_used >= b_used)
  1122.       prod_ptr->used = Kara_Mult(a.PTR->vec, a_used, b.PTR->vec, b_used,
  1123.                                  prod_ptr->vec);
  1124.   else
  1125.       prod_ptr->used = Kara_Mult(b.PTR->vec, b_used, a.PTR->vec, a_used,
  1126.                                  prod_ptr->vec);
  1127.   return prod_ptr;
  1128. }
  1133. integer operator/(const integer & a, const integer & b) 
  1134. { 
  1135.   sz_t a_sign = a.PTR->sign;
  1136.   sz_t b_sign = b.PTR->sign;
  1137.   sz_t a_used = a.PTR->used;
  1138.   sz_t b_used = b.PTR->used;
  1140.   if (b_sign == ZERO) error_handler(1,"division by zero");
  1142.   if (a_sign == ZERO) return a;
  1144.   int rel = absolute_cmp(a.PTR->vec,a_used, b.PTR->vec, b_used);
  1146.   if (rel == 0) return (a_sign == b_sign) ? 1 : -1;  // |a| == |b|
  1147.   if (rel  < 0) return 0;                            // |a| <  |b|
  1149.   // |a| > |b|
  1151.   int sz = next_power(a_used - b_used + 1);
  1153.   integer_rep* quot_ptr = new_integer_rep(sz);
  1155.   quot_ptr->sign = (a_sign == b_sign) ? POSITIVE : NEGATIVE;
  1157.   if (b_used == 1) 
  1158.       quot_ptr->used = DivByOneWord(a.PTR->vec, a_used, b.PTR->vec[0],
  1159.                                     quot_ptr->vec,0);
  1160.   else 
  1161.       // |a| and |b| both have at least 2 digits
  1162.       quot_ptr->used = School_Div(a.PTR->vec, a_used, b.PTR->vec, b_used,
  1163.                                   quot_ptr->vec,0,0);
  1165.   return quot_ptr;
  1166. }
  1171. integer operator%(const integer & a, const integer & b) 
  1172. { 
  1173.   sz_t a_sign = a.PTR->sign;
  1174.   sz_t b_sign = b.PTR->sign;
  1175.   sz_t a_used = a.PTR->used;
  1176.   sz_t b_used = b.PTR->used;
  1178.   if (b_sign == ZERO) error_handler(1,"mod by zero");
  1180.   if (a_sign == ZERO) return a;
  1182.   int rel = absolute_cmp(a.PTR->vec,a_used, b.PTR->vec, b_used);
  1184.   if (rel == 0) return 0;  // |a| == |b|
  1185.   if (rel  < 0) return a;  // |a| <  |b|
  1187.   // |a| > |b|
  1189.   int qsz = next_power(a_used - b_used + 1);
  1190.   int rsz = b.PTR->size;
  1192.   integer_rep* quot_ptr = new_integer_rep(qsz);
  1193.   integer_rep* rem_ptr = new_integer_rep(rsz);
  1196.   if (b_used == 1)
  1197.     { quot_ptr->used = DivByOneWord(a.PTR->vec, a_used, b.PTR->vec[0],
  1198.                                     quot_ptr->vec, rem_ptr->vec);
  1199.       rem_ptr->used = (rem_ptr->vec[0] == 0) ? 0 : 1;
  1200.      }
  1201.   else 
  1202.     { sz_t r_used;
  1203.       School_Div(a.PTR->vec, a_used, b.PTR->vec, b_used,
  1204.                  quot_ptr->vec,rem_ptr->vec,&r_used);
  1205.       rem_ptr->used = r_used;
  1206.      }
  1208.   delete_integer_rep(quot_ptr);
  1210.   if (rem_ptr->used == 0) 
  1211.      rem_ptr->sign = ZERO;
  1212.   else
  1213.      rem_ptr->sign = a_sign;
  1215.   return rem_ptr;
  1216. }
  1221. integer operator & (const integer & a, const integer & b)   // bitwise and
  1222. {
  1223.   sz_t a_sign = a.PTR->sign;
  1224.   sz_t b_sign = b.PTR->sign;
  1226.   if (a_sign == 0) return a;
  1227.   if (b_sign == 0) return b;
  1229.   sz_t a_used = a.PTR->used;
  1230.   sz_t b_used = b.PTR->used;
  1232.   integer_rep* and_ptr;
  1234.   word* ap = a.PTR->vec;
  1235.   word* bp = b.PTR->vec;
  1236.   word* a_stop = ap + a_used; 
  1237.   word* b_stop = bp + b_used; 
  1239.   int used;
  1240.   word* p;
  1242.   if (a_used <= b_used)
  1243.      { and_ptr = new_integer_rep(a.PTR->size);
  1244.        used = a_used;
  1245.        p = and_ptr->vec;
  1246.        while (ap < a_stop) *p++ = *ap++ & *bp++;
  1247.      }
  1248.   else
  1249.      { and_ptr = new_integer_rep(b.PTR->size);
  1250.        used = b_used;
  1251.        p = and_ptr->vec;
  1252.        while (bp < b_stop) *p++ = *ap++ & *bp++;
  1253.       }
  1255.   while (*--p == 0 && used > 0) used--;
  1258.   and_ptr->used = used;
  1259.   and_ptr->sign = (used > 0) ? a_sign : ZERO;
  1262.   return and_ptr;
  1263. }
  1264.  
  1266. integer operator | (const integer & a, const integer & b)   // bitwise or
  1267. {
  1268.   sz_t a_sign = a.PTR->sign;
  1269.   sz_t b_sign = b.PTR->sign;
  1271.   if (a_sign == 0) return a;
  1272.   if (b_sign == 0) return b;
  1274.   sz_t a_used = a.PTR->used;
  1275.   sz_t b_used = b.PTR->used;
  1277.   integer_rep* and_ptr;
  1279.   word* ap = a.PTR->vec;
  1280.   word* bp = b.PTR->vec;
  1281.   word* a_stop = ap + a_used; 
  1282.   word* b_stop = bp + b_used; 
  1284.   if (a_used >= b_used)
  1285.      { and_ptr = new_integer_rep(a.PTR->size);
  1286.        and_ptr->used = a_used;
  1287.        and_ptr->sign = a_sign;
  1288.        word* p = and_ptr->vec;
  1289.        while (bp < b_stop) *p++ = *ap++ | *bp++;
  1290.        while (ap < a_stop) *p++ = *ap++;
  1291.      }
  1292.   else
  1293.      { and_ptr = new_integer_rep(b.PTR->size);
  1294.        and_ptr->used = b_used;
  1295.        and_ptr->sign = a_sign;
  1296.        word* p = and_ptr->vec;
  1297.        while (ap < a_stop) *p++ = *ap++ | *bp++;
  1298.        while (bp < b_stop) *p++ = *bp++;
  1299.       }
  1301.   and_ptr->sign = a_sign;
  1303.   return and_ptr;
  1304. }
  1305.  
  1308. bool operator > (const integer & a, const integer & b)
  1309. {
  1310.   if (a.PTR == b.PTR) return false;
  1312.   int a_sign = a.PTR->sign;
  1313.   int b_sign = b.PTR->sign;
  1315.   if (a_sign > b_sign) return true;
  1316.   if (a_sign < b_sign) return false;
  1318.   /* the signs are equal */
  1320.   if (a_sign == ZERO) return false;
  1322.   if (a_sign == POSITIVE)
  1323.     return (absolute_cmp(a.PTR->vec,a.PTR->used,
  1324.                          b.PTR->vec,b.PTR->used) > 0);
  1325.   else
  1326.     return (absolute_cmp(a.PTR->vec,a.PTR->used,
  1327.                          b.PTR->vec,b.PTR->used) < 0);
  1328. }
  1331. bool operator < (const integer & a, const integer & b)
  1332. {
  1333.   if (a.PTR == b.PTR) return false;
  1335.   int a_sign = a.PTR->sign;
  1336.   int b_sign = b.PTR->sign;
  1338.   if (a_sign < b_sign) return true;
  1339.   if (a_sign > b_sign) return false;
  1341.   /* the signs are equal */
  1343.   if (a_sign == ZERO) return false;
  1345.   if (a_sign == POSITIVE)
  1346.     return (absolute_cmp(a.PTR->vec,a.PTR->used,
  1347.                          b.PTR->vec,b.PTR->used) < 0);
  1348.   else
  1349.     return (absolute_cmp(a.PTR->vec,a.PTR->used,
  1350.                          b.PTR->vec,b.PTR->used) > 0);
  1351. }
  1354. bool operator == (const integer & a, const integer & b)
  1355. {
  1356.   if (a.PTR == b.PTR) return true;
  1358.   int a_sign = a.PTR->sign;
  1359.   int b_sign = b.PTR->sign;
  1361.   if (a_sign != b_sign) return false;
  1363.   if (a_sign == ZERO) return true;
  1365.   return (absolute_cmp(a.PTR->vec,a.PTR->used,
  1366.                        b.PTR->vec,b.PTR->used) == 0);
  1367. }
  1372. integer gcd(const integer& a, const integer& b)
  1373. { 
  1374.   // gcd is not efficient. See Knuth 2 for improvements. 
  1376.   int a_sign = a.PTR->sign;
  1377.   int b_sign = b.PTR->sign;
  1378.   
  1379.   if (a_sign == ZERO)
  1380.     if (b_sign == ZERO)
  1381.         return 1;
  1382.     else
  1383.         return abs(b);
  1386.   // a is non-zero
  1388.   if (b_sign == 0) return abs(a);
  1389.   
  1391.   // both a and b are non-zero
  1393.   integer u = abs(a);
  1394.   integer v = abs(b);
  1396.   if (u < v) v = v%u;
  1397.    
  1398.   while (sign(v) != 0)
  1399.   { integer tmp = u % v; 
  1400.     u = v;
  1401.     v = tmp;
  1402.    }
  1404.   return u;
  1405. }
  1409. //------------------------------------------------------------------------------
  1410. // member functions
  1411. //------------------------------------------------------------------------------
  1413. void integer::absolute() // *this = abs(*this), nur effizienter
  1414. {
  1415.   if (PTR->count == 1) PTR->sign = POSITIVE;
  1416.   else 
  1417.     *this = abs(*this); 
  1418. }
  1420. int integer::zeros() const
  1421. { // gives the number of zeros at the end of the integer
  1422.   if (PTR->sign == ZERO) return 0;
  1423.   int k=0;
  1424.   while (!(PTR->vec[k])) k++;
  1425.   int len = k * WORD_LENGTH;
  1426.   word low = PTR->vec[k];
  1427.   if (! (low & 65535)) { low >>= 16; len += 16; }
  1428.   if (! (low & 255))  { low >>= 8; len += 8; }
  1429.   if (! (low & 15))  { low >>= 4; len += 4; }
  1430.   if (! (low & 3))  { low >>= 2; len += 2; }
  1431.   if (! (low & 1))  { low >>= 1; len += 1; }
  1432.   return len;
  1433. }
  1435. integer integer::div(const integer & b, integer& r) 
  1436. { 
  1437.   // returns quotient and assings remainder to r
  1439.   sz_t a_sign = PTR->sign;
  1440.   sz_t b_sign = b.PTR->sign;
  1441.   sz_t a_used = PTR->used;
  1442.   sz_t b_used = b.PTR->used;
  1444.   if (b_sign == ZERO) error_handler(1,"division by zero");
  1446.   if (a_sign == ZERO) 
  1447.   { r = 0;
  1448.     return *this;
  1449.    }
  1451.   int rel = absolute_cmp(PTR->vec,a_used, b.PTR->vec, b_used);
  1453.   if (rel == 0) 
  1454.   { r = 0;
  1455.     return (a_sign == b_sign) ? 1 : -1; 
  1456.    }
  1458.   if (rel  < 0) 
  1459.   { r = *this;
  1460.     return 0;                     
  1461.    }
  1463.   // |a| > |b|
  1465.   int sz = next_power(a_used - b_used + 1);
  1467.   integer_rep* quot_ptr = new_integer_rep(sz);
  1469.   quot_ptr->sign = (a_sign == b_sign) ? POSITIVE : NEGATIVE;
  1471.   integer_rep* rem_ptr = new_integer_rep(b.PTR->size);
  1473.   sz_t r_used;
  1475.   if (b_used == 1)
  1476.     { quot_ptr->used = DivByOneWord(PTR->vec, a_used, b.PTR->vec[0],
  1477.                                     quot_ptr->vec, rem_ptr->vec);
  1478.       r_used = (rem_ptr->vec[0] == 0) ? 0 : 1;
  1479.      }
  1480.   else 
  1481.      School_Div(PTR->vec, a_used, b.PTR->vec, b_used,
  1482.                 quot_ptr->vec,rem_ptr->vec,&r_used);
  1484.   rem_ptr->used = r_used;
  1485.   rem_ptr->sign = (r_used == 0) ? ZERO : a_sign;
  1487.   r = rem_ptr;
  1489.   return quot_ptr;
  1490. }
  1494. int integer::length()  const   // number of bits
  1495. { if (PTR->sign == ZERO) return 0;
  1496.   int len = PTR->used - 1;
  1497.   word hi= PTR->vec[len];
  1498.   len *= WORD_LENGTH;
  1499.   while(hi)
  1500.   { hi >>= 1;
  1501.     len++;
  1502.    }
  1503.   return len;
  1504.  }
  1505.   
  1507. integer integer::operator-()  const
  1508. { // unary minus
  1509.   if (PTR->sign == ZERO) return *this;
  1510.   integer_rep* p = copy_integer_rep(PTR);
  1511.   p->sign = -PTR->sign;
  1512.   return p;
  1513. }
  1516. integer integer::operator~()  const
  1517. { // negation
  1518.   integer_rep* p = copy_integer_rep(PTR);
  1519.   word* wp = p->vec + p->used;
  1520.   while (--wp >= p->vec) *wp = ~*wp;
  1521.   return p;
  1522. }
  1526. bool integer::operator==(int l) const
  1527. { int sig = PTR->sign;
  1528.   if (l==0) return sig == ZERO;
  1530.   if (PTR->used > 1) return false;
  1532.   if (l > 0 )
  1533.   { if (sig <= 0) return false;
  1534.     return PTR->vec[0] == word(l);
  1535.    }
  1537.   // l < 0
  1538.   if (sig >= 0) return false;
  1539.   return PTR->vec[0] == word(-l);
  1540. }
  1543. bool integer::operator<(int l) const
  1544. { int sig = PTR->sign;
  1545.   if (l==0) return sig == NEGATIVE;
  1547.   if (l > 0 )
  1548.   { if (sig <= 0) return true;
  1549.     if (PTR->used > 1) return false;
  1550.     return PTR->vec[0] < word(l);
  1551.    }
  1552.   
  1553.   // l < 0
  1554.   if (sig >= 0) return false;
  1555.   if (PTR->used > 1) return true;
  1556.   return PTR->vec[0] > word(-l);
  1557. }
  1560. bool integer::operator>(int l) const
  1561. { int sig = PTR->sign;
  1562.   if (l==0) return sig == POSITIVE;
  1564.   if (l < 0 )
  1565.   { if (sig >= 0) return true;
  1566.     if (PTR->used > 1) return false;
  1567.     return PTR->vec[0] < word(-l);
  1568.    }
  1569.   
  1570.   // l > 0
  1571.   if (sig <= 0) return false;
  1572.   if (PTR->used > 1) return true;
  1573.   return PTR->vec[0] > word(l);
  1574. }
  1577. double integer::todouble() const
  1578. { 
  1579.   int i = PTR->used;
  1580.   if (i == 0) return 0;
  1581.   if (i > 32) return 2*MAXDOUBLE; //Infinity
  1583.   word*  v = PTR->vec;
  1584.   double d = 0;
  1586.   while (i > 0) d = pow32*d + v[--i];
  1588.   return (PTR->sign == NEGATIVE) ? -d : d;
  1589. }
  1592. integer integer::sqrt() const
  1593. {
  1594.   if (PTR->sign == NEGATIVE)
  1595.     error_handler(1, "negative argument in integer::sqrt");
  1597.   if (PTR->sign == ZERO) return *this;
  1599.   integer root = (*this) >> (length()/2);  // first approx.
  1601.   integer fix  = (root*root + *this)/(root << 1);
  1603.   while (fix != root)
  1604.   { fix  = (root*root + *this)/(root << 1);
  1605.     root = (fix*fix   + *this)/(fix << 1);
  1606.   }
  1608.   return root;  
  1609. }
  1612. integer integer::operator<<(long n)  const
  1613. {
  1614.   if (PTR->sign == ZERO || n == 0) return *this;
  1616.   int w_shift = int(n / WORD_LENGTH);
  1617.   int b_shift = int(n % WORD_LENGTH);
  1618.   int used = PTR->used + w_shift;
  1619.   int r = WORD_LENGTH - b_shift;
  1621.   if (b_shift && (PTR->vec[PTR->used-1] >> r)) used++;
  1623.   int sz = next_power(used);
  1625.   integer_rep* result;
  1626.   result = new_integer_rep(sz);
  1627.   result->used = used;
  1628.   result->sign = PTR->sign;
  1630.   word* q = PTR->vec;
  1631.   word* q_stop = q + PTR->used;
  1632.   word* p = result->vec;
  1634.   while (w_shift--) *p++ = 0;
  1636.   if (b_shift) 
  1637.     { word carry = 0;
  1638.       while (q < q_stop)
  1639.       { word B = *q++;
  1640.         carry |= (B << b_shift);
  1641.         *p++ = carry;
  1642.         carry = (B >> r);
  1643.        }
  1644.       if (carry) *p = carry;
  1645.      }
  1646.   else
  1647.      while (q < q_stop) *p++ = *q++;
  1649.   return result;
  1650. }
  1654. integer integer::operator>>(long n)  const
  1655. {
  1656.   if (n <= 0 || PTR->sign == ZERO) return *this;
  1658.   int sticks = int(n / WORD_LENGTH);
  1659.   int shifts = int(n % WORD_LENGTH);
  1660.   int used   = PTR->used - sticks;
  1661.  
  1662.   int sz = next_power(used);
  1663.  
  1664.   integer_rep* result;
  1665.  
  1666.   result = new_integer_rep(sz);
  1667.   result->sign = PTR->sign;
  1668.  
  1669.   word* p = result->vec + used;
  1670.   word* q = PTR->vec + PTR->used;
  1671.  
  1672.   if (shifts != 0)
  1673.     { word carry = 0;
  1674.       int r = WORD_LENGTH - shifts;
  1675.       while (p > result->vec)
  1676.       { word B = *--q;
  1677.         carry |= (B >> shifts);
  1678.         *--p = carry;
  1679.         carry = (B << r);
  1680.        }
  1681.      }
  1682.   else
  1683.      while (p > result->vec) *--p = *--q;
  1684.  
  1685.   if (result->vec[used-1] == 0) used--;
  1686.   if (used == 0) result->sign = ZERO;
  1687.  
  1688.   result->used = used;
  1689.  
  1690.   return result;
  1691. }
  1692.  
  1695. //------------------------------------------------------------------------------
  1696. // stream i/o
  1697. //------------------------------------------------------------------------------
  1700. void integer::hex_print(ostream& out)
  1701. { if (PTR->sign == ZERO)
  1702.   { out << "0";
  1703.     return;
  1704.    }
  1705.   word* p = PTR->vec + PTR->used;
  1706.   out << string("%x",*--p);
  1707.   while (p > PTR->vec) out << string("%0x",*--p);
  1708. }
  1712. ostream & operator << (ostream & out, const integer & x) 
  1713. {
  1714.   if (x.PTR->sign == ZERO) 
  1715.   { out << "0";
  1716.     return out;
  1717.    }
  1719.   if (x.PTR->sign == NEGATIVE) out << "-";
  1721.   Print_Vector(x.PTR->vec,x.PTR->used,out);
  1723.   return out;
  1724. }
  1728. istream & operator >> (istream & in, integer & a) 
  1729. {
  1730.   bool negative = false;
  1731.   const int null = '0';
  1732.   char c;
  1734.   while (in.get(c) && isspace(c));
  1736.   if (c == '-') 
  1737.   { negative = true;
  1738.     while (in.get(c) && isspace(c));
  1739.    }
  1741.   if (isdigit(c)) 
  1742.   { a = c - '0';
  1743.     while (in.get(c) && isdigit(c)) a = 10*a + (c-null);
  1744.    }
  1746.   if (sign(a) != 0 && negative) a = -a;
  1748.   return in;
  1750. }
  1754. //------------------------------------------------------------------------------
  1755. // miscellaneous 
  1756. //------------------------------------------------------------------------------
  1759. integer_rep* copy_integer_rep(integer_rep* x)
  1760. { integer_rep* result = new_integer_rep(x->size);
  1761.   result->sign = x->sign;
  1762.   result->used = x->used;
  1763.   Copy_Vector(result->vec,x->vec,x->size);
  1764.   return result;
  1765. }
  1768. integer abs(const integer & a)
  1769. { if (a.PTR->sign >= 0) return a;
  1770.          integer_rep* ptr = copy_integer_rep(a.PTR);
  1771.          ptr->sign = POSITIVE;
  1772.          return ptr;
  1773. }
  1776. int log(const integer & a)
  1777. {
  1778.   if (a.PTR->sign != POSITIVE)
  1779.         error_handler(999,"domain error in log(integer)");
  1781.   int l  = a.PTR->used - 1;
  1782.   int lg = l * WORD_LENGTH;
  1784.   word hi = a.PTR->vec[l];
  1786.   while (hi >>= 1) lg++;
  1788.   return lg;
  1789. }
  1792. int integer::cmp(const integer & a, const integer & b)
  1793. { int a_sign = a.PTR->sign;
  1794.   int b_sign = b.PTR->sign;
  1796.   if (a_sign < b_sign) return -1;
  1797.   if (a_sign > b_sign) return  1;
  1799.   // the signs are equal 
  1801.   if (a_sign == ZERO) return 0;
  1803.   if (a_sign == POSITIVE)
  1804.       return  absolute_cmp(a.PTR->vec,a.PTR->used,
  1805.                            b.PTR->vec,b.PTR->used);
  1806.   else
  1807.       return -absolute_cmp(a.PTR->vec,a.PTR->used,
  1808.                            b.PTR->vec,b.PTR->used);
  1809. }
  1813. integer integer::random(int n)  // return n-bit random integer
  1814. { 
  1815.   sz_t w = n / WORD_LENGTH;
  1816.   sz_t r = n % WORD_LENGTH;
  1817.   if (r) w++;
  1819.   int sz = next_power(w);
  1821.   integer_rep* ptr = new_integer_rep(sz);
  1822.   ptr->sign = POSITIVE;
  1823.   ptr->used = w;
  1825.   for (int i = 0; i < w; i++) rand_int >> ptr->vec[i];
  1827.   if (r) ptr->vec[w-1] >>= (WORD_LENGTH-r);
  1829.   return ptr;
  1830.  }