home *** CD-ROM | disk | FTP | other *** search
/ ftp.cs.arizona.edu / ftp.cs.arizona.edu.tar / ftp.cs.arizona.edu / icon / historic / v941.tgz / icon.v941src.tar / icon.v941src / src / runtime / rlrgint.r < prev    next >
Text File  |  2002-01-18  |  53KB  |  2,318 lines
  1. /*
  2.  * File: rlrgint.r
  3.  *  Large integer arithmetic
  4.  */
  5.  
  6. #ifdef LargeInts
  7.  
  8. extern int over_flow;
  9.  
  10. /*
  11.  *  Conventions:
  12.  *
  13.  *  Lrgints entering this module and leaving it are too large to
  14.  *  be represented with T_Integer.  So, externally, a given value
  15.  *  is always T_Integer or always T_Lrgint.
  16.  *
  17.  *  Routines outside this module operate on bignums by calling
  18.  *  a routine like
  19.  *
  20.  *      bigadd(da, db, dx)
  21.  *
  22.  *  where da, db, and dx are pointers to tended descriptors.
  23.  *  For the common case where one argument is a T_Integer, these
  24.  *  call routines like
  25.  *
  26.  *      bigaddi(da, IntVal(*db), dx).
  27.  *
  28.  *  The bigxxxi routines can convert an integer to bignum form;
  29.  *  they use itobig.
  30.  *
  31.  *  The routines that actually do the work take (length, address)
  32.  *  pairs specifying unsigned base-B digit strings.  The sign handling
  33.  *  is done in the bigxxx routines.
  34.  */
  35.  
  36. /*
  37.  * Type for doing arithmetic on (2 * NB)-bit nonnegative numbers.
  38.  *  Normally unsigned but may be signed (with NB reduced appropriately)
  39.  *  if unsigned arithmetic is slow.
  40.  */
  41.  
  42. /* The bignum radix, B */
  43.  
  44. #define B            ((word)1 << NB)
  45.  
  46. /* Lrgint digits in a word */
  47.  
  48. #define WORDLEN  (WordBits / NB + (WordBits % NB != 0))
  49.  
  50. /* size of a bignum block that will hold an integer */
  51.  
  52. #define INTBIGBLK  sizeof(struct b_bignum) + sizeof(DIGIT) * WORDLEN
  53.  
  54. /* lo(uword d) :            the low digit of a uword
  55.    hi(uword d) :            the rest, d is unsigned
  56.    signed_hi(uword d) :     the rest, d is signed
  57.    dbl(DIGIT a, DIGIT b) : the two-digit uword [a,b] */
  58.  
  59. #define lo(d)        ((d) & (B - 1))
  60. #define hi(d)        ((uword)(d) >> NB)
  61. #define dbl(a,b)     (((uword)(a) << NB) + (b))
  62.  
  63. #if ((-1) >> 1) < 0
  64. #define signed_hi(d) ((word)(d) >> NB)
  65. #else
  66. #define signbit      ((uword)1 << (WordBits - NB - 1))
  67. #define signed_hi(d) ((word)((((uword)(d) >> NB) ^ signbit) - signbit))
  68. #endif
  69.  
  70. /* LrgInt(dptr dp) : the struct b_bignum pointed to by dp */
  71.  
  72. #define LrgInt(dp)   ((struct b_bignum *)&BlkLoc(*dp)->bignumblk)
  73.  
  74. /* LEN(struct b_bignum *b) : number of significant digits */
  75.  
  76. #define LEN(b)       ((b)->lsd - (b)->msd + 1)
  77.  
  78. /* DIG(struct b_bignum *b, word i): pointer to ith most significant digit */
  79. /*  (NOTE: This macro expansion often results in a very long string,
  80.  *   that has been known to cause problems with the C compiler under VMS.
  81.  *   So when DIG is used, keep it to one use per line.)
  82.  */
  83.  
  84. #define DIG(b,i)     (&(b)->digits[(b)->msd+(i)])
  85.  
  86. /* ceil, ln: ceil may be 1 too high in case ln is inaccurate */
  87.  
  88. #undef ceil
  89. #define ceil(x)      ((word)((x) + 1.01))
  90. #define ln(n)        (log((double)n))
  91.  
  92. /* determine the number of words needed for a bignum block with n digits */
  93.  
  94. #define LrgNeed(n)   ( ((sizeof(struct b_bignum) + ((n) - 1) * sizeof(DIGIT)) \
  95.                + WordSize - 1) & -WordSize )
  96.  
  97. /* copied from rconv.c */
  98.  
  99. #define tonum(c)     (isdigit(c) ? (c)-'0' : 10+(((c)|(040))-'a'))
  100.  
  101. /* copied from oref.c */
  102.  
  103. #define RandVal (RanScale*(k_random=(RandA*(long)k_random+RandC)&0x7fffffffL))
  104.  
  105. /*
  106.  * Prototypes.
  107.  */
  108.  
  109. static int mkdesc    (struct b_bignum *x, dptr dx);
  110. static void itobig    (word i, struct b_bignum *x, dptr dx);
  111.  
  112. static void decout    (FILE *f, DIGIT *n, word l);
  113.  
  114. static int bigaddi    (dptr da, word i, dptr dx);
  115. static int bigsubi    (dptr da, word i, dptr dx);
  116. static int bigmuli    (dptr da, word i, dptr dx);
  117. static int bigdivi    (dptr da, word i, dptr dx);
  118. static int bigmodi    (dptr da, word i, dptr dx);
  119. static int bigpowi    (dptr da, word i, dptr dx);
  120. static int bigpowii    (word a, word i, dptr dx);
  121. static word bigcmpi    (dptr da, word i);
  122.  
  123. static DIGIT add1    (DIGIT *u, DIGIT *v, DIGIT *w, word n);
  124. static word sub1    (DIGIT *u, DIGIT *v, DIGIT *w, word n);
  125. static void mul1    (DIGIT *u, DIGIT *v, DIGIT *w, word n, word m);
  126. static int div1
  127.    (DIGIT *a, DIGIT *b, DIGIT *q, DIGIT *r, word m, word n, struct b_bignum *b1, struct b_bignum *b2);
  128. static void compl1    (DIGIT *u, DIGIT *w, word n);
  129. static word cmp1    (DIGIT *u, DIGIT *v, word n);
  130. static DIGIT addi1    (DIGIT *u, word k, DIGIT *w, word n);
  131. static void subi1    (DIGIT *u, word k, DIGIT *w, word n);
  132. static DIGIT muli1    (DIGIT *u, word k, int c, DIGIT *w, word n);
  133. static DIGIT divi1    (DIGIT *u, word k, DIGIT *w, word n);
  134. static DIGIT shifti1    (DIGIT *u, word k, DIGIT c, DIGIT *w, word n);
  135. static word cmpi1    (DIGIT *u, word k, word n);
  136.  
  137. #define bdzero(dest,l)  memset(dest, '\0', (l) * sizeof(DIGIT))
  138. #define bdcopy(src, dest, l)  memcpy(dest, src, (l) * sizeof(DIGIT))
  139.  
  140. /*
  141.  * mkdesc -- put value into a descriptor
  142.  */
  143.  
  144. static int mkdesc(x, dx)
  145. struct b_bignum *x;
  146. dptr dx;
  147. {
  148.    word xlen, cmp;
  149.    static DIGIT maxword[WORDLEN] = { 1 << ((WordBits - 1) % NB) };
  150.  
  151.    /* suppress leading zero digits */
  152.  
  153.    while (x->msd != x->lsd &&
  154.           *DIG(x,0) == 0)
  155.       x->msd++;
  156.  
  157.    /* put it into a word if it fits, otherwise return the bignum */
  158.  
  159.    xlen = LEN(x);
  160.  
  161.    if (xlen < WORDLEN ||
  162.        (xlen == WORDLEN &&
  163.         ((cmp = cmp1(DIG(x,0), maxword, (word)WORDLEN)) < 0 ||
  164.         (cmp == (word)0 && x->sign)))) {
  165.       word val = -(word)*DIG(x,0);
  166.       word i;
  167.  
  168.       for (i = x->msd; ++i <= x->lsd; )
  169.          val = (val << NB) - x->digits[i];
  170.       if (!x->sign)
  171.      val = -val;
  172.       dx->dword = D_Integer;
  173.       IntVal(*dx) = val;
  174.       }
  175.    else {
  176.       dx->dword = D_Lrgint;
  177.       BlkLoc(*dx) = (union block *)x;
  178.       }
  179.    return Succeeded;
  180. }
  181.  
  182. /*
  183.  *  i -> big
  184.  */
  185.  
  186. static void itobig(i, x, dx)
  187. word i;
  188. struct b_bignum *x;
  189. dptr dx;
  190. {
  191.    x->lsd = WORDLEN - 1;
  192.    x->msd = WORDLEN;
  193.    x->sign = 0;
  194.  
  195.    if (i == 0) {
  196.       x->msd--;
  197.       *DIG(x,0) = 0;
  198.       }
  199.    else if (i < 0) {
  200.       word d = lo(i);
  201.  
  202.       if (d != 0) {
  203.          d = B - d;
  204.          i += B;
  205.          }
  206.       i = - signed_hi(i);
  207.       x->msd--;
  208.       *DIG(x,0) = d;
  209.       x->sign = 1;
  210.       }
  211.  
  212.    while (i != 0) {
  213.       x->msd--;
  214.       *DIG(x,0) = lo(i);
  215.       i = hi(i);
  216.       }
  217.  
  218.    dx->dword = D_Lrgint;
  219.    BlkLoc(*dx) = (union block *)x;
  220. }
  221.  
  222. /*
  223.  *  string -> bignum
  224.  */
  225.  
  226. word bigradix(sign, r, s, end_s, result)
  227. int sign;                      /* '-' or not */
  228. int r;                          /* radix 2 .. 36 */
  229. char *s, *end_s;                        /* input string */
  230. union numeric *result;          /* output T_Integer or T_Lrgint */
  231. {
  232.    struct b_bignum *b;
  233.    DIGIT *bd;
  234.    word len;
  235.    int c;
  236.  
  237.    if (r == 0)
  238.       return CvtFail;
  239.    len = ceil((end_s - s) * ln(r) / ln(B));
  240.    Protect(b = alcbignum(len), return Error);
  241.    bd = DIG(b,0);
  242.  
  243.    bdzero(bd, len);
  244.  
  245.    if (r < 2 || r > 36)
  246.       return CvtFail;
  247.  
  248.    for (c = ((s < end_s) ? *s++ : ' '); isalnum(c);
  249.         c = ((s < end_s) ? *s++ : ' ')) {
  250.       c = tonum(c);
  251.       if (c >= r)
  252.      return CvtFail;
  253.       muli1(bd, (word)r, c, bd, len);
  254.       }
  255.  
  256.    /*
  257.     * Skip trailing white space and make sure there is nothing else left
  258.     *  in the string. Note, if we have already reached end-of-string,
  259.     *  c has been set to a space.
  260.     */
  261.    while (isspace(c) && s < end_s)
  262.       c = *s++;
  263.    if (!isspace(c))
  264.       return CvtFail;
  265.  
  266.    if (sign == '-')
  267.       b->sign = 1;
  268.  
  269.    /* put value into dx and return the type */
  270.  
  271.    { struct descrip dx;
  272.    (void)mkdesc(b, &dx);
  273.    if (Type(dx) == T_Lrgint)
  274.      result->big = (struct b_bignum *)BlkLoc(dx);
  275.    else
  276.      result->integer = IntVal(dx);
  277.    return Type(dx);
  278.    }
  279. }
  280.  
  281. /*
  282.  *  bignum -> real
  283.  */
  284.  
  285. double bigtoreal(da)
  286. dptr da;
  287. {
  288.    word i;
  289.    double r = 0;
  290.    struct b_bignum *b = &BlkLoc(*da)->bignumblk;
  291.  
  292.    for (i = b->msd; i <= b->lsd; i++)
  293.       r = r * B + b->digits[i];
  294.  
  295.    return (b->sign ? -r : r);
  296. }
  297.  
  298. /*
  299.  *  real -> bignum
  300.  */
  301.  
  302. int realtobig(da, dx)
  303. dptr da, dx;
  304. {
  305.  
  306. #ifdef Double
  307.    double x;
  308. #else                    /* Double */
  309.    double x = BlkLoc(*da)->realblk.realval;
  310. #endif                    /* Double */
  311.  
  312.    struct b_bignum *b;
  313.    word i, blen;
  314.    word d;
  315.    int sgn;
  316.  
  317. #ifdef Double
  318.     {
  319.         int    *rp, *rq;
  320.         rp = (int *) &(BlkLoc(*da)->realblk.realval);
  321.         rq = (int *) &x;
  322.         *rq++ = *rp++;
  323.         *rq = *rp;
  324.     }
  325. #endif                    /* Double */
  326.  
  327.    if (x > 0.9999 * MinLong && x < 0.9999 * MaxLong) {
  328.       MakeInt((word)x, dx);
  329.       return Succeeded;        /* got lucky; a simple integer suffices */
  330.       }
  331.  
  332.    if (sgn = x < 0)
  333.       x = -x;
  334.    blen = ln(x) / ln(B) + 0.99;
  335.    for (i = 0; i < blen; i++)
  336.       x /= B;
  337.    if (x >= 1.0) {
  338.       x /= B;
  339.       blen += 1;
  340.       }
  341.  
  342.    Protect(b = alcbignum(blen), return Error);
  343.    for (i = 0; i < blen; i++) {
  344.       d = (x *= B);
  345.       *DIG(b,i) = d;
  346.       x -= d;
  347.       }
  348.  
  349.    b->sign = sgn;
  350.    return mkdesc(b, dx);
  351. }
  352.  
  353. /*
  354.  *  bignum -> string
  355.  */
  356.  
  357. int bigtos(da, dx)
  358. dptr da, dx;
  359. {
  360.    tended struct b_bignum *a, *temp;
  361.    word alen = LEN(LrgInt(da));
  362.    word slen = ceil(alen * ln(B) / ln(10));
  363.    char *p, *q;
  364.  
  365.    a = LrgInt(da);
  366.    Protect(temp = alcbignum(alen), fatalerr(0,NULL));
  367.    if (a->sign)
  368.       slen++;
  369.    Protect(q = alcstr(NULL,slen), fatalerr(0,NULL));
  370.    bdcopy(DIG(a,0),
  371.           DIG(temp,0),
  372.           alen);
  373.    p = q += slen;
  374. #ifdef VMS
  375.    {
  376.       DIGIT *tempdg;
  377.       for (;;) {
  378.          tempdg = DIG(temp,0);
  379.          if (!cmpi1(tempdg, (word)0, alen))
  380.             break;
  381.          *--p = '0' + divi1(tempdg, (word)10, tempdg, alen);
  382.       }
  383.    }
  384. #else
  385.    while (cmpi1(DIG(temp,0),
  386.                 (word)0, alen))
  387.       *--p = '0' + divi1(DIG(temp,0),
  388.                          (word)10,
  389.                          DIG(temp,0),
  390.                          alen);
  391. #endif            /* VMS */
  392.    if (a->sign)
  393.       *--p = '-';
  394.    StrLen(*dx) = q - p;
  395.    StrLoc(*dx) = p;
  396.    return NoCvt;    /* The mnemonic is wrong, but the signal means */
  397.             /* that the string is allocated and not null- */
  398.             /* terminated. */
  399. }
  400.  
  401. /*
  402.  *  bignum -> file
  403.  */
  404.  
  405. void bigprint(f, da)
  406. FILE *f;
  407. dptr da;
  408. {
  409.    struct b_bignum *a, *temp;
  410.    word alen = LEN(LrgInt(da));
  411.    word slen, dlen;
  412.    struct b_bignum *blk = &BlkLoc(*da)->bignumblk;
  413.  
  414.    slen = blk->lsd - blk->msd;
  415.    dlen = slen * NB * 0.3010299956639812    /* 1 / log2(10) */
  416.       + log((double)blk->digits[blk->msd]) * 0.4342944819032518 + 0.5;
  417.                         /* 1 / ln(10) */
  418.    if (dlen >= MaxDigits) {
  419.       fprintf(f, "integer(~10^%ld)",(long)dlen);
  420.       return;
  421.       }
  422.  
  423.    /* not worth passing this one back */
  424.    Protect(temp = alcbignum(alen), fatalerr(0, NULL));
  425.  
  426.    a = LrgInt(da);
  427.    bdcopy(DIG(a,0),
  428.           DIG(temp,0),
  429.           alen);
  430.    if (a->sign)
  431.       putc('-', f);
  432.    decout(f,
  433.           DIG(temp,0),
  434.           alen);
  435. }
  436.  
  437. /*
  438.  * decout - given a base B digit string, print the number in base 10.
  439.  */
  440. static void decout(f, n, l)
  441. FILE *f;
  442. DIGIT *n;
  443. word l;
  444. {
  445.    DIGIT i = divi1(n, (word)10, n, l);
  446.  
  447.    if (cmpi1(n, (word)0, l))
  448.       decout(f, n, l);
  449.    putc('0' + i, f);
  450. }
  451.  
  452. /*
  453.  *  da -> dx
  454.  */
  455.  
  456. int cpbignum(da, dx)
  457. dptr da, dx;
  458. {
  459.    struct b_bignum *a, *x;
  460.    word alen = LEN(LrgInt(da));
  461.  
  462.    Protect(x = alcbignum(alen), return Error);
  463.    a = LrgInt(da);
  464.    bdcopy(DIG(a,0),
  465.           DIG(x,0),
  466.           alen);
  467.    x->sign = a->sign;
  468.    return mkdesc(x, dx);
  469. }
  470.  
  471. /*
  472.  *  da + db -> dx
  473.  */
  474.  
  475. int bigadd(da, db, dx)
  476. dptr da, db;
  477. dptr dx;
  478. {
  479.    tended struct b_bignum *a, *b;
  480.    struct b_bignum *x;
  481.    word alen, blen;
  482.    word c;
  483.  
  484.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  485.       alen = LEN(LrgInt(da));
  486.       blen = LEN(LrgInt(db));
  487.       a = LrgInt(da);
  488.       b = LrgInt(db);
  489.       if (a->sign == b->sign) {
  490.          if (alen > blen) {
  491.             Protect(x = alcbignum(alen + 1), return Error);
  492.             c = add1(DIG(a,alen-blen),
  493.                      DIG(b,0),
  494.                      DIG(x,alen-blen+1),
  495.                      blen);
  496.             *DIG(x,0) =
  497.                 addi1(DIG(a,0),
  498.                       c,
  499.                       DIG(x,1),
  500.                       alen-blen);
  501.             }
  502.          else if (alen == blen) {
  503.             Protect(x = alcbignum(alen + 1), return Error);
  504.             *DIG(x,0) =
  505.                add1(DIG(a,0),
  506.                     DIG(b,0),
  507.                     DIG(x,1),
  508.                     alen);
  509.             }
  510.          else {
  511.             Protect(x = alcbignum(blen + 1), return Error);
  512.             c = add1(DIG(b,blen-alen),
  513.                      DIG(a,0),
  514.                      DIG(x,blen-alen+1),
  515.                      alen);
  516.             *DIG(x,0) =
  517.                addi1(DIG(b,0),
  518.                      c,
  519.                      DIG(x,1),
  520.                      blen-alen);
  521.             }
  522.          x->sign = a->sign;
  523.          }
  524.       else {
  525.          if (alen > blen) {
  526.             Protect(x = alcbignum(alen), return Error);
  527.             c = sub1(DIG(a,alen-blen),
  528.                      DIG(b,0),
  529.                      DIG(x,alen-blen),
  530.                      blen);
  531.             subi1(DIG(a,0),
  532.                   -c,
  533.                   DIG(x,0),
  534.                   alen-blen);
  535.             x->sign = a->sign;
  536.             }
  537.          else if (alen == blen) {
  538.             Protect(x = alcbignum(alen), return Error);
  539.             if (cmp1(DIG(a,0),
  540.                      DIG(b,0),
  541.                      alen) > 0) {
  542.                (void)sub1(DIG(a,0),
  543.                           DIG(b,0),
  544.                           DIG(x,0),
  545.                           alen);
  546.                x->sign = a->sign;
  547.                }
  548.             else {
  549.                (void)sub1(DIG(b,0),
  550.                           DIG(a,0),
  551.                           DIG(x,0),
  552.                           alen);
  553.                x->sign = b->sign;
  554.                }
  555.             }
  556.          else {
  557.             Protect(x = alcbignum(blen), return Error);
  558.             c = sub1(DIG(b,blen-alen),
  559.                      DIG(a,0),
  560.                      DIG(x,blen-alen),
  561.                      alen);
  562.             subi1(DIG(b,0),
  563.                   -c,
  564.                   DIG(x,0),
  565.                   blen-alen);
  566.             x->sign = b->sign;
  567.             }
  568.          }
  569.       return mkdesc(x, dx);
  570.       }
  571.    else if (Type(*da) == T_Lrgint)    /* bignum + integer */
  572.       return bigaddi(da, IntVal(*db), dx);
  573.    else if (Type(*db) == T_Lrgint)    /* integer + bignum */
  574.       return bigaddi(db, IntVal(*da), dx);
  575.    else {                             /* integer + integer */
  576.       struct descrip td;
  577.       char tdigits[INTBIGBLK];
  578.  
  579.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  580.       return bigaddi(&td, IntVal(*db), dx);
  581.       }
  582. }
  583.  
  584. /*
  585.  *  da - db -> dx
  586.  */
  587.  
  588. int bigsub(da, db, dx)
  589. dptr da, db, dx;
  590. {
  591.    struct descrip td;
  592.    char tdigits[INTBIGBLK];
  593.    tended struct b_bignum *a, *b;
  594.    struct b_bignum *x;
  595.    word alen, blen;
  596.    word c;
  597.  
  598.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  599.       alen = LEN(LrgInt(da));
  600.       blen = LEN(LrgInt(db));
  601.       a = LrgInt(da);
  602.       b = LrgInt(db);
  603.       if (a->sign != b->sign) {
  604.          if (alen > blen) {
  605.             Protect(x = alcbignum(alen + 1), return Error);
  606.             c = add1(DIG(a,alen-blen),
  607.                      DIG(b,0),
  608.                      DIG(x,alen-blen+1),
  609.                      blen);
  610.             *DIG(x,0) =
  611.                addi1(DIG(a,0),
  612.                      c,
  613.                      DIG(x,1),
  614.                      alen-blen);
  615.             }
  616.          else if (alen == blen) {
  617.             Protect(x = alcbignum(alen + 1), return Error);
  618.             *DIG(x,0) =
  619.                add1(DIG(a,0),
  620.                     DIG(b,0),
  621.                     DIG(x,1),
  622.                     alen);
  623.             }
  624.          else {
  625.             Protect(x = alcbignum(blen + 1), return Error);
  626.             c = add1(DIG(b,blen-alen),
  627.                      DIG(a,0),
  628.                      DIG(x,blen-alen+1),
  629.                      alen);
  630.             *DIG(x,0) =
  631.                addi1(DIG(b,0),
  632.                      c,
  633.                      DIG(x,1),
  634.                      blen-alen);
  635.             }
  636.          x->sign = a->sign;
  637.          }
  638.       else {
  639.          if (alen > blen) {
  640.             Protect(x = alcbignum(alen), return Error);
  641.             c = sub1(DIG(a,alen-blen),
  642.                      DIG(b,0),
  643.                      DIG(x,alen-blen),
  644.                      blen);
  645.             subi1(DIG(a,0),
  646.                   -c,
  647.                   DIG(x,0),
  648.                   alen-blen);
  649.             x->sign = a->sign;
  650.             }
  651.          else if (alen == blen) {
  652.             Protect(x = alcbignum(alen), return Error);
  653.             if (cmp1(DIG(a,0),
  654.                      DIG(b,0),
  655.                      alen) > 0) {
  656.                (void)sub1(DIG(a,0),
  657.                           DIG(b,0),
  658.                           DIG(x,0),
  659.                           alen);
  660.                x->sign = a->sign;
  661.                }
  662.             else {
  663.                (void)sub1(DIG(b,0),
  664.                           DIG(a,0),
  665.                           DIG(x,0),
  666.                           alen);
  667.                x->sign = 1 ^ b->sign;
  668.                }
  669.             }
  670.          else {
  671.             Protect(x = alcbignum(blen), return Error);
  672.             c = sub1(DIG(b,blen-alen),
  673.                      DIG(a,0),
  674.                      DIG(x,blen-alen),
  675.                      alen);
  676.             subi1(DIG(b,0),
  677.                       -c,
  678.                       DIG(x,0),
  679.                       blen-alen);
  680.             x->sign = 1 ^ b->sign;
  681.             }
  682.          }
  683.       return mkdesc(x, dx);
  684.       }
  685.    else if (Type(*da) == T_Lrgint)     /* bignum - integer */
  686.       return bigsubi(da, IntVal(*db), dx);
  687.    else if (Type(*db) == T_Lrgint) {   /* integer - bignum */
  688.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  689.       alen = LEN(LrgInt(&td));
  690.       blen = LEN(LrgInt(db));
  691.       a = LrgInt(&td);
  692.       b = LrgInt(db);
  693.       if (a->sign != b->sign) {
  694.          if (alen == blen) {
  695.             Protect(x = alcbignum(alen + 1), return Error);
  696.             *DIG(x,0) =
  697.                add1(DIG(a,0),
  698.                     DIG(b,0),
  699.                     DIG(x,1),
  700.                     alen);
  701.             }
  702.          else {
  703.             Protect(x = alcbignum(blen + 1), return Error);
  704.             c = add1(DIG(b,blen-alen),
  705.                      DIG(a,0),
  706.                      DIG(x,blen-alen+1),
  707.                      alen);
  708.             *DIG(x,0) =
  709.                addi1(DIG(b,0),
  710.                      c,
  711.                      DIG(x,1),
  712.                      blen-alen);
  713.             }
  714.          x->sign = a->sign;
  715.          }
  716.       else {
  717.          if (alen == blen) {
  718.             Protect(x = alcbignum(alen), return Error);
  719.             if (cmp1(DIG(a,0),
  720.                      DIG(b,0),
  721.                      alen) > 0) {
  722.                (void)sub1(DIG(a,0),
  723.                           DIG(b,0),
  724.                           DIG(x,0),
  725.                           alen);
  726.                x->sign = a->sign;
  727.                }
  728.             else {
  729.                (void)sub1(DIG(b,0),
  730.                           DIG(a,0),
  731.                           DIG(x,0),
  732.                           alen);
  733.                x->sign = 1 ^ b->sign;
  734.                }
  735.             }
  736.          else {
  737.             Protect(x = alcbignum(blen), return Error);
  738.             c = sub1(DIG(b,blen-alen),
  739.                      DIG(a,0),
  740.                      DIG(x,blen-alen),
  741.                      alen);
  742.             subi1(DIG(b,0),
  743.                   -c,
  744.                   DIG(x,0),
  745.                   blen-alen);
  746.             x->sign = 1 ^ b->sign;
  747.             }
  748.          }
  749.       return mkdesc(x, dx);
  750.       }
  751.    else {                              /* integer - integer */
  752.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  753.       return bigsubi(&td, IntVal(*db), dx);
  754.       }
  755.  
  756. }
  757.  
  758. /*
  759.  *  da * db -> dx
  760.  */
  761.  
  762. int bigmul(da, db, dx)
  763. dptr da, db, dx;
  764. {
  765.    tended struct b_bignum *a, *b;
  766.    struct b_bignum *x;
  767.    word alen, blen;
  768.  
  769.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  770.       alen = LEN(LrgInt(da));
  771.       blen = LEN(LrgInt(db));
  772.       a = LrgInt(da);
  773.       b = LrgInt(db);
  774.       Protect(x = alcbignum(alen + blen), return Error);
  775.       mul1(DIG(a,0),
  776.            DIG(b,0),
  777.            DIG(x,0),
  778.            alen, blen);
  779.       x->sign = a->sign ^ b->sign;
  780.       return mkdesc(x, dx);
  781.       }
  782.    else if (Type(*da) == T_Lrgint)    /* bignum * integer */
  783.       return bigmuli(da, IntVal(*db), dx);
  784.    else if (Type(*db) == T_Lrgint)    /* integer * bignum */
  785.       return bigmuli(db, IntVal(*da), dx);
  786.    else {                             /* integer * integer */
  787.       struct descrip td;
  788.       char tdigits[INTBIGBLK];
  789.  
  790.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  791.       return bigmuli(&td, IntVal(*db), dx);
  792.       }
  793. }
  794.  
  795. /*
  796.  *  da / db -> dx
  797.  */
  798.  
  799. int bigdiv(da, db, dx)
  800. dptr da, db, dx;
  801. {
  802.    tended struct b_bignum *a, *b, *x, *tu, *tv;
  803.    word alen, blen;
  804.    struct descrip td;
  805.    char tdigits[INTBIGBLK];
  806.  
  807.    /* Put *da into large integer format. */
  808.    if (Type(*da) != T_Lrgint) {
  809.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  810.       da = &td;
  811.       }
  812.  
  813.    if (Type(*db) == T_Lrgint) {         /* bignum / bignum */
  814.       alen = LEN(LrgInt(da));
  815.       blen = LEN(LrgInt(db));
  816.       if (alen < blen) {
  817.          MakeInt(0, dx);
  818.          return Succeeded;
  819.          }
  820.       a = LrgInt(da);
  821.       b = LrgInt(db);
  822.       Protect(x = alcbignum(alen - blen + 1), return Error);
  823.       if (blen == 1)
  824.          divi1(DIG(a,0),
  825.                (word)*DIG(b,0),
  826.                DIG(x,0),
  827.                alen);
  828.       else {
  829.          Protect(tu = alcbignum(alen + 1), return Error);
  830.          Protect(tv = alcbignum(blen), return Error);
  831.          if (div1(DIG(a,0),
  832.                   DIG(b,0),
  833.                   DIG(x,0),
  834.                   NULL, alen-blen, blen, tu, tv) == Error)
  835.             return Error;
  836.             }
  837.       x->sign = a->sign ^ b->sign;
  838.       return mkdesc(x, dx);
  839.       }
  840.    else                                 /* bignum / integer */
  841.       return bigdivi(da, IntVal(*db), dx);
  842. }
  843.  
  844. /*
  845.  *  da % db -> dx
  846.  */
  847.  
  848. int bigmod(da, db, dx)
  849. dptr da, db, dx;
  850. {
  851.    tended struct b_bignum *a, *b, *x, *temp, *tu, *tv;
  852.    word alen, blen;
  853.    struct descrip td;
  854.    char tdigits[INTBIGBLK];
  855.  
  856.    /* Put *da into large integer format. */
  857.    if (Type(*da) != T_Lrgint) {
  858.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  859.       da = &td;
  860.       }
  861.  
  862.    if (Type(*db) == T_Lrgint) {        /* bignum % bignum */
  863.       alen = LEN(LrgInt(da));
  864.       blen = LEN(LrgInt(db));
  865.       if (alen < blen) {
  866.          cpbignum(da, dx);
  867.          return Succeeded;
  868.          }
  869.       a = LrgInt(da);
  870.       b = LrgInt(db);
  871.       Protect(x = alcbignum(blen), return Error);
  872.       if (blen == 1) {
  873.      Protect(temp = alcbignum(alen), return Error);
  874.          *DIG(x,0) =
  875.             divi1(DIG(a,0),
  876.                   (word)*DIG(b,0),
  877.                   DIG(temp,0),
  878.                   alen);
  879.          }
  880.       else {
  881.          Protect(tu = alcbignum(alen + 1), return Error);
  882.          Protect(tv = alcbignum(blen), return Error);
  883.          if (div1(DIG(a,0),
  884.                   DIG(b,0),
  885.                   NULL,
  886.                   DIG(x,0),
  887.                   alen-blen, blen, tu, tv) == Error)
  888.             return Error;
  889.             }
  890.       x->sign = a->sign;
  891.       return mkdesc(x, dx);
  892.       }
  893.    else                       /* bignum % integer */
  894.       return bigmodi(da, IntVal(*db), dx);
  895. }
  896.  
  897. /*
  898.  *  -i -> dx
  899.  */
  900.  
  901. int bigneg(da, dx)
  902. dptr da, dx;
  903. {
  904.    struct descrip td;
  905.    char tdigits[INTBIGBLK];
  906.    int cpstat;
  907.  
  908.    /* Put *da into large integer format. */
  909.    if (Type(*da) != T_Lrgint) {
  910.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  911.       da = &td;
  912.       }
  913.    LrgInt(da)->sign ^= 1;       /* Temporarily change the sign */
  914.    cpstat = cpbignum(da, dx);
  915.    LrgInt(da)->sign ^= 1;       /* Change it back */
  916.  
  917.    return cpstat;
  918. }
  919.  
  920. /*
  921.  *  da ^ db -> dx
  922.  */
  923.  
  924. int bigpow(da, db, dx)
  925. dptr da, db, dx;
  926. {
  927.  
  928.    if (Type(*db) == T_Lrgint) {
  929.       struct b_bignum *b;
  930.  
  931.       b = LrgInt ( db );
  932.  
  933.  
  934.       if (Type(*da) == T_Lrgint) {
  935.      if ( b->sign ) {
  936.         /* bignum ^ -bignum = 0 */
  937.         MakeInt ( 0, dx );
  938.         return Succeeded;
  939.         }
  940.      else
  941.         /* bignum ^ +bignum = guaranteed overflow */
  942.         ReturnErrNum(307, Error);
  943.      }
  944.       else if ( b->sign )
  945.      /* integer ^ -bignum */
  946.      switch ( IntVal ( *da ) ) {
  947.         case 1:
  948.            MakeInt ( 1, dx );
  949.            return Succeeded;
  950.         case -1:
  951.            /* Result is +1 / -1, depending on whether *b is even or odd. */
  952.            if ( ( b->digits[ b->lsd ] ) & 01 )
  953.           MakeInt ( -1, dx );
  954.            else
  955.           MakeInt ( 1, dx );
  956.            return Succeeded;
  957.         case 0:
  958.            ReturnErrNum(204,Error);
  959.         default:
  960.            /* da ^ (negative int) = 0 for all non-special cases */
  961.            MakeInt(0, dx);
  962.            return Succeeded;
  963.         }
  964.       else {
  965.      /* integer ^ +bignum */
  966.      word n, blen;
  967.      register DIGIT nth_dig, mask;
  968.  
  969.      b = LrgInt ( db );
  970.      blen = LEN ( b );
  971.  
  972.      /* We scan the bits of b from the most to least significant.
  973.       * The bit position in b is represented by the pair ( n, mask )
  974.       * where n is the DIGIT number (0 = most sig.) and mask is the
  975.       * the bit mask for the current bit.
  976.       *
  977.       * For each bit (most sig to least) in b,
  978.       *  for each zero, square the partial result;
  979.       *  for each one, square it and multiply it by a */
  980.      MakeInt ( 1, dx );
  981.      for ( n = 0; n < blen; ++n ) {
  982.         nth_dig = *DIG ( b, n );
  983.         for ( mask = 1 << ( NB - 1 ); mask; mask >>= 1 ) {
  984.            if ( bigmul ( dx, dx, dx ) == Error )
  985.           return Error;
  986.            if ( nth_dig & mask )
  987.           if ( bigmul ( dx, da, dx ) == Error )
  988.              return Error;
  989.            }
  990.         }
  991.      }
  992.       return Succeeded;
  993.       }
  994.    else if (Type(*da) == T_Lrgint)    /* bignum ^ integer */
  995.       return bigpowi(da, IntVal(*db), dx);
  996.    else                               /* integer ^ integer */
  997.       return bigpowii(IntVal(*da), IntVal(*db), dx);
  998. }
  999.  
  1000. int bigpowri( a, db, drslt )
  1001. double a;
  1002. dptr   db, drslt;
  1003. {
  1004.    register double retval;
  1005.    register word n;
  1006.    register DIGIT nth_dig, mask;
  1007.    struct b_bignum *b;
  1008.    word blen;
  1009.  
  1010.    b = LrgInt ( db );
  1011.    blen = LEN ( b );
  1012.    if ( b->sign ) {
  1013.       if ( a == 0.0 )
  1014.      ReturnErrNum(204, Error);
  1015.       else
  1016.      a = 1.0 / a;
  1017.       }
  1018.  
  1019.    /* We scan the bits of b from the most to least significant.
  1020.     * The bit position in b is represented by the pair ( n, mask )
  1021.     * where n is the DIGIT number (0 = most sig.) and mask is the
  1022.     * the bit mask for the current bit.
  1023.     *
  1024.     * For each bit (most sig to least) in b,
  1025.     *  for each zero, square the partial result;
  1026.     *  for each one, square it and multiply it by a */
  1027.    retval = 1.0;
  1028.    for ( n = 0; n < blen; ++n ) {
  1029.       nth_dig = *DIG ( b, n );
  1030.       for ( mask = 1 << ( NB - 1 ); mask; mask >>= 1 ) {
  1031.      retval *= retval;
  1032.      if ( nth_dig & mask )
  1033.         retval *= a;
  1034.      }
  1035.       }
  1036.  
  1037.    Protect(BlkLoc(*drslt) = (union block *)alcreal(retval), return Error);
  1038.    drslt->dword = D_Real;
  1039.    return Succeeded;
  1040. }
  1041.  
  1042. /*
  1043.  *  iand(da, db) -> dx
  1044.  */
  1045.  
  1046. int bigand(da, db, dx)
  1047. dptr da, db, dx;
  1048. {
  1049.    tended struct b_bignum *a, *b, *x, *tad, *tbd;
  1050.    word alen, blen, xlen;
  1051.    word i;
  1052.    DIGIT *ad, *bd;
  1053.    struct descrip td;
  1054.    char tdigits[INTBIGBLK];
  1055.  
  1056.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  1057.       alen = LEN(LrgInt(da));
  1058.       blen = LEN(LrgInt(db));
  1059.       xlen = alen > blen ? alen : blen;
  1060.       a = LrgInt(da);
  1061.       b = LrgInt(db);
  1062.       Protect(x = alcbignum(xlen), return Error);
  1063.  
  1064.       if (alen == xlen && !a->sign)
  1065.          ad = DIG(a,0);
  1066.       else {
  1067.          Protect(tad = alcbignum(xlen), return Error);
  1068.          ad = DIG(tad,0);
  1069.          bdzero(ad, xlen - alen);
  1070.          bdcopy(DIG(a,0),
  1071.                 &ad[xlen-alen], alen);
  1072.          if (a->sign)
  1073.         compl1(ad, ad, xlen);
  1074.          }
  1075.  
  1076.       if (blen == xlen && !b->sign)
  1077.          bd = DIG(b,0);
  1078.       else {
  1079.          Protect(tbd = alcbignum(xlen), return Error);
  1080.          bd = DIG(tbd,0);
  1081.          bdzero(bd, xlen - blen);
  1082.          bdcopy(DIG(b,0),
  1083.                 &bd[xlen-blen], blen);
  1084.          if (b->sign)
  1085.         compl1(bd, bd, xlen);
  1086.          }
  1087.  
  1088.       for (i = 0; i < xlen; i++)
  1089.          *DIG(x,i) =
  1090.             ad[i] & bd[i];
  1091.  
  1092.       if (a->sign & b->sign) {
  1093.          x->sign = 1;
  1094.          compl1(DIG(x,0),
  1095.                 DIG(x,0),
  1096.                 xlen);
  1097.          }
  1098.       }
  1099.    else if (Type(*da) == T_Lrgint) {   /* iand(bignum,integer) */
  1100.       itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
  1101.       alen = LEN(LrgInt(da));
  1102.       blen = LEN(LrgInt(&td));
  1103.       xlen = alen > blen ? alen : blen;
  1104.       a = LrgInt(da);
  1105.       b = LrgInt(&td);
  1106.       Protect(x = alcbignum(alen), return Error);
  1107.  
  1108.       if (alen == xlen && !a->sign)
  1109.          ad = DIG(a,0);
  1110.       else {
  1111.          Protect(tad = alcbignum(xlen), return Error);
  1112.          ad = DIG(tad,0);
  1113.          bdzero(ad, xlen - alen);
  1114.          bdcopy(DIG(a,0),
  1115.                 &ad[xlen-alen], alen);
  1116.          if (a->sign)
  1117.         compl1(ad, ad, xlen);
  1118.          }
  1119.  
  1120.       if (blen == xlen && !b->sign)
  1121.          bd = DIG(b,0);
  1122.       else {
  1123.          Protect(tbd = alcbignum(xlen), return Error);
  1124.          bd = DIG(tbd,0);
  1125.          bdzero(bd, xlen - blen);
  1126.          bdcopy(DIG(b,0),
  1127.                 &bd[xlen-blen], blen);
  1128.          if (b->sign)
  1129.         compl1(bd, bd, xlen);
  1130.          }
  1131.  
  1132.       for (i = 0; i < xlen; i++)
  1133.          *DIG(x,i) =
  1134.             ad[i] & bd[i];
  1135.  
  1136.       if (a->sign & b->sign) {
  1137.          x->sign = 1;
  1138.          compl1(DIG(x,0),
  1139.                 DIG(x,0),
  1140.                 xlen);
  1141.          }
  1142.       }
  1143.    else if (Type(*db) == T_Lrgint) {   /* iand(integer,bignum) */
  1144.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  1145.       alen = LEN(LrgInt(&td));
  1146.       blen = LEN(LrgInt(db));
  1147.       xlen = alen > blen ? alen : blen;
  1148.       a = LrgInt(&td);
  1149.       b = LrgInt(db);
  1150.       Protect(x = alcbignum(blen), return Error);
  1151.  
  1152.       if (alen == xlen && !a->sign)
  1153.          ad = DIG(a,0);
  1154.       else {
  1155.          Protect(tad = alcbignum(xlen), return Error);
  1156.          ad = DIG(tad,0);
  1157.          bdzero(ad, xlen - alen);
  1158.          bdcopy(DIG(a,0),
  1159.                 &ad[xlen-alen], alen);
  1160.          if (a->sign)
  1161.         compl1(ad, ad, xlen);
  1162.          }
  1163.  
  1164.       if (blen == xlen && !b->sign)
  1165.          bd = DIG(b,0);
  1166.       else {
  1167.          Protect(tbd = alcbignum(xlen), return Error);
  1168.          bd = DIG(tbd,0);
  1169.          bdzero(bd, xlen - blen);
  1170.          bdcopy(DIG(b,0),
  1171.                 &bd[xlen-blen], blen);
  1172.          if (b->sign)
  1173.         compl1(bd, bd, xlen);
  1174.          }
  1175.  
  1176.       for (i = 0; i < xlen; i++)
  1177.          *DIG(x,i) =
  1178.             ad[i] & bd[i];
  1179.  
  1180.       if (a->sign & b->sign) {
  1181.          x->sign = 1;
  1182.          compl1(DIG(x,0),
  1183.                 DIG(x,0),
  1184.                 xlen);
  1185.          }
  1186.       }
  1187.    /* not called for iand(integer,integer) */
  1188.  
  1189.    return mkdesc(x, dx);
  1190. }
  1191.  
  1192. /*
  1193.  *  ior(da, db) -> dx
  1194.  */
  1195.  
  1196. int bigor(da, db, dx)
  1197. dptr da, db, dx;
  1198. {
  1199.    tended struct b_bignum *a, *b, *x, *tad, *tbd;
  1200.    word alen, blen, xlen;
  1201.    word i;
  1202.    DIGIT *ad, *bd;
  1203.    struct descrip td;
  1204.    char tdigits[INTBIGBLK];
  1205.  
  1206.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  1207.       alen = LEN(LrgInt(da));
  1208.       blen = LEN(LrgInt(db));
  1209.       xlen = alen > blen ? alen : blen;
  1210.       a = LrgInt(da);
  1211.       b = LrgInt(db);
  1212.       Protect(x = alcbignum(xlen), return Error);
  1213.  
  1214.       if (alen == xlen && !a->sign)
  1215.          ad = DIG(a,0);
  1216.       else {
  1217.          Protect(tad = alcbignum(xlen), return Error);
  1218.          ad = DIG(tad,0);
  1219.          bdzero(ad, xlen - alen);
  1220.          bdcopy(DIG(a,0),
  1221.                 &ad[xlen-alen], alen);
  1222.          if (a->sign)
  1223.         compl1(ad, ad, xlen);
  1224.          }
  1225.  
  1226.       if (blen == xlen && !b->sign)
  1227.          bd = DIG(b,0);
  1228.       else {
  1229.          Protect(tbd = alcbignum(xlen), return Error);
  1230.          bd = DIG(tbd,0);
  1231.          bdzero(bd, xlen - blen);
  1232.          bdcopy(DIG(b,0),
  1233.                 &bd[xlen-blen], blen);
  1234.          if (b->sign)
  1235.         compl1(bd, bd, xlen);
  1236.          }
  1237.  
  1238.       for (i = 0; i < xlen; i++)
  1239.          *DIG(x,i) =
  1240.             ad[i] | bd[i];
  1241.  
  1242.       if (a->sign | b->sign) {
  1243.          x->sign = 1;
  1244.          compl1(DIG(x,0),
  1245.                 DIG(x,0),
  1246.                 xlen);
  1247.          }
  1248.       }
  1249.    else if (Type(*da) == T_Lrgint) {   /* ior(bignum,integer) */
  1250.       itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
  1251.       alen = LEN(LrgInt(da));
  1252.       blen = LEN(LrgInt(&td));
  1253.       xlen = alen > blen ? alen : blen;
  1254.       a = LrgInt(da);
  1255.       b = LrgInt(&td);
  1256.       Protect(x = alcbignum(alen), return Error);
  1257.  
  1258.       if (alen == xlen && !a->sign)
  1259.          ad = DIG(a,0);
  1260.       else {
  1261.          Protect(tad = alcbignum(xlen), return Error);
  1262.          ad = DIG(tad,0);
  1263.          bdzero(ad, xlen - alen);
  1264.          bdcopy(DIG(a,0),
  1265.                 &ad[xlen-alen], alen);
  1266.          if (a->sign)
  1267.         compl1(ad, ad, xlen);
  1268.          }
  1269.  
  1270.       if (blen == xlen && !b->sign)
  1271.          bd = DIG(b,0);
  1272.       else {
  1273.          Protect(tbd = alcbignum(xlen), return Error);
  1274.          bd = DIG(tbd,0);
  1275.          bdzero(bd, xlen - blen);
  1276.          bdcopy(DIG(b,0),
  1277.                 &bd[xlen-blen], blen);
  1278.          if (b->sign)
  1279.         compl1(bd, bd, xlen);
  1280.          }
  1281.  
  1282.       for (i = 0; i < xlen; i++)
  1283.          *DIG(x,i) =
  1284.             ad[i] | bd[i];
  1285.  
  1286.       if (a->sign | b->sign) {
  1287.          x->sign = 1;
  1288.          compl1(DIG(x,0),
  1289.                 DIG(x,0),
  1290.                 xlen);
  1291.          }
  1292.       }
  1293.    else if (Type(*db) == T_Lrgint) {   /* ior(integer,bignym) */
  1294.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  1295.       alen = LEN(LrgInt(&td));
  1296.       blen = LEN(LrgInt(db));
  1297.       xlen = alen > blen ? alen : blen;
  1298.       a = LrgInt(&td);
  1299.       b = LrgInt(db);
  1300.       Protect(x = alcbignum(blen), return Error);
  1301.  
  1302.       if (alen == xlen && !a->sign)
  1303.          ad = DIG(a,0);
  1304.       else {
  1305.          Protect(tad = alcbignum(xlen), return Error);
  1306.          ad = DIG(tad,0);
  1307.          bdzero(ad, xlen - alen);
  1308.          bdcopy(DIG(a,0),
  1309.                 &ad[xlen-alen], alen);
  1310.          if (a->sign)
  1311.         compl1(ad, ad, xlen);
  1312.          }
  1313.  
  1314.       if (blen == xlen && !b->sign)
  1315.          bd = DIG(b,0);
  1316.       else {
  1317.          Protect(tbd = alcbignum(xlen), return Error);
  1318.          bd = DIG(tbd,0);
  1319.          bdzero(bd, xlen - blen);
  1320.          bdcopy(DIG(b,0),
  1321.                 &bd[xlen-blen], blen);
  1322.          if (b->sign)
  1323.         compl1(bd, bd, xlen);
  1324.          }
  1325.  
  1326.       for (i = 0; i < xlen; i++)
  1327.          *DIG(x,i) =
  1328.             ad[i] | bd[i];
  1329.  
  1330.       if (a->sign | b->sign) {
  1331.          x->sign = 1;
  1332.          compl1(DIG(x,0),
  1333.                 DIG(x,0),
  1334.                 xlen);
  1335.          }
  1336.       }
  1337.    /* not called for ior(integer,integer) */
  1338.  
  1339.    return mkdesc(x, dx);
  1340. }
  1341.  
  1342. /*
  1343.  *  xor(da, db) -> dx
  1344.  */
  1345.  
  1346. int bigxor(da, db, dx)
  1347. dptr da, db, dx;
  1348. {
  1349.    tended struct b_bignum *a, *b, *x, *tad, *tbd;
  1350.    word alen, blen, xlen;
  1351.    word i;
  1352.    DIGIT *ad, *bd;
  1353.    struct descrip td;
  1354.    char tdigits[INTBIGBLK];
  1355.  
  1356.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  1357.       alen = LEN(LrgInt(da));
  1358.       blen = LEN(LrgInt(db));
  1359.       xlen = alen > blen ? alen : blen;
  1360.       a = LrgInt(da);
  1361.       b = LrgInt(db);
  1362.       Protect(x = alcbignum(xlen), return Error);
  1363.  
  1364.       if (alen == xlen && !a->sign)
  1365.          ad = DIG(a,0);
  1366.       else {
  1367.          Protect(tad = alcbignum(xlen), return Error);
  1368.          ad = DIG(tad,0);
  1369.          bdzero(ad, xlen - alen);
  1370.          bdcopy(DIG(a,0),
  1371.                 &ad[xlen-alen], alen);
  1372.          if (a->sign)
  1373.         compl1(ad, ad, xlen);
  1374.          }
  1375.  
  1376.       if (blen == xlen && !b->sign)
  1377.          bd = DIG(b,0);
  1378.       else {
  1379.          Protect(tbd = alcbignum(xlen), return Error);
  1380.          bd = DIG(tbd,0);
  1381.          bdzero(bd, xlen - blen);
  1382.          bdcopy(DIG(b,0),
  1383.                 &bd[xlen-blen], blen);
  1384.          if (b->sign)
  1385.         compl1(bd, bd, xlen);
  1386.          }
  1387.  
  1388.       for (i = 0; i < xlen; i++)
  1389.          *DIG(x,i) =
  1390.             ad[i] ^ bd[i];
  1391.  
  1392.       if (a->sign ^ b->sign) {
  1393.          x->sign = 1;
  1394.          compl1(DIG(x,0),
  1395.                 DIG(x,0),
  1396.                 xlen);
  1397.          }
  1398.       }
  1399.    else if (Type(*da) == T_Lrgint) {   /* ixor(bignum,integer) */
  1400.       itobig(IntVal(*db), (struct b_bignum *)tdigits, &td);
  1401.       alen = LEN(LrgInt(da));
  1402.       blen = LEN(LrgInt(&td));
  1403.       xlen = alen > blen ? alen : blen;
  1404.       a = LrgInt(da);
  1405.       b = LrgInt(&td);
  1406.       Protect(x = alcbignum(alen), return Error);
  1407.  
  1408.       if (alen == xlen && !a->sign)
  1409.          ad = DIG(a,0);
  1410.       else {
  1411.          Protect(tad = alcbignum(xlen), return Error);
  1412.          ad = DIG(tad,0);
  1413.          bdzero(ad, xlen - alen);
  1414.          bdcopy(DIG(a,0),
  1415.                 &ad[xlen-alen], alen);
  1416.          if (a->sign)
  1417.         compl1(ad, ad, xlen);
  1418.          }
  1419.  
  1420.       if (blen == xlen && !b->sign)
  1421.          bd = DIG(b,0);
  1422.       else {
  1423.          Protect(tbd = alcbignum(xlen), return Error);
  1424.          bd = DIG(tbd,0);
  1425.          bdzero(bd, xlen - blen);
  1426.          bdcopy(DIG(b,0),
  1427.                 &bd[xlen-blen], blen);
  1428.          if (b->sign)
  1429.         compl1(bd, bd, xlen);
  1430.          }
  1431.  
  1432.       for (i = 0; i < xlen; i++)
  1433.          *DIG(x,i) =
  1434.             ad[i] ^ bd[i];
  1435.  
  1436.       if (a->sign ^ b->sign) {
  1437.          x->sign = 1;
  1438.          compl1(DIG(x,0),
  1439.                 DIG(x,0),
  1440.                 xlen);
  1441.          }
  1442.       }
  1443.    else if (Type(*db) == T_Lrgint) {   /* ixor(integer,bignum) */
  1444.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  1445.       alen = LEN(LrgInt(&td));
  1446.       blen = LEN(LrgInt(db));
  1447.       xlen = alen > blen ? alen : blen;
  1448.       a = LrgInt(&td);
  1449.       b = LrgInt(db);
  1450.       Protect(x = alcbignum(blen), return Error);
  1451.  
  1452.       if (alen == xlen && !a->sign)
  1453.          ad = DIG(a,0);
  1454.       else {
  1455.          Protect(tad = alcbignum(xlen), return Error);
  1456.          ad = DIG(tad,0);
  1457.          bdzero(ad, xlen - alen);
  1458.          bdcopy(DIG(a,0),
  1459.                 &ad[xlen-alen], alen);
  1460.          if (a->sign)
  1461.         compl1(ad, ad, xlen);
  1462.          }
  1463.  
  1464.       if (blen == xlen && !b->sign)
  1465.          bd = DIG(b,0);
  1466.       else {
  1467.          Protect(tbd = alcbignum(xlen), return Error);
  1468.          bd = DIG(tbd,0);
  1469.          bdzero(bd, xlen - blen);
  1470.          bdcopy(DIG(b,0),
  1471.                 &bd[xlen-blen], blen);
  1472.          if (b->sign)
  1473.         compl1(bd, bd, xlen);
  1474.          }
  1475.  
  1476.       for (i = 0; i < xlen; i++)
  1477.          *DIG(x,i) =
  1478.             ad[i] ^ bd[i];
  1479.  
  1480.       if (a->sign ^ b->sign) {
  1481.          x->sign = 1;
  1482.          compl1(DIG(x,0),
  1483.                 DIG(x,0),
  1484.                 xlen);
  1485.          }
  1486.       }
  1487.    /* not called for ixor(integer,integer) */
  1488.  
  1489.    return mkdesc(x, dx);
  1490. }
  1491.  
  1492. /*
  1493.  *  bigshift(da, db) -> dx
  1494.  */
  1495.  
  1496. int bigshift(da, db, dx)
  1497. dptr da, db, dx;
  1498. {
  1499.    tended struct b_bignum *a, *x, *tad;
  1500.    word alen;
  1501.    word r = IntVal(*db) % NB;
  1502.    word q = (r >= 0 ? IntVal(*db) : (IntVal(*db) - (r += NB))) / NB;
  1503.    word xlen;
  1504.    DIGIT *ad;
  1505.    struct descrip td;
  1506.    char tdigits[INTBIGBLK];
  1507.  
  1508.    if (Type(*da) == T_Integer) {
  1509.       itobig(IntVal(*da), (struct b_bignum *)tdigits, &td);
  1510.       da = &td;
  1511.       }
  1512.  
  1513.    alen = LEN(LrgInt(da));
  1514.    xlen = alen + q + 1;
  1515.    if (xlen <= 0) {
  1516.       MakeInt(-LrgInt(da)->sign, dx);
  1517.       return Succeeded;
  1518.       }
  1519.    else {
  1520.       a = LrgInt(da);
  1521.       Protect(x = alcbignum(xlen), return Error);
  1522.  
  1523.       if (a->sign) {
  1524.          Protect(tad = alcbignum(alen), return Error);
  1525.          ad = DIG(tad,0);
  1526.          bdcopy(DIG(a,0),
  1527.                 ad, alen);
  1528.          compl1(ad, ad, alen);
  1529.          }
  1530.       else
  1531.          ad = DIG(a,0);
  1532.  
  1533.       if (q >= 0) {
  1534.          *DIG(x,0) =
  1535.             shifti1(ad, r, (DIGIT)0,
  1536.                     DIG(x,1),
  1537.                     alen);
  1538.          bdzero(DIG(x,alen+1),
  1539.                 q);
  1540.          }
  1541.       else
  1542.          *DIG(x,0) =
  1543.             shifti1(ad, r, ad[alen+q] >> (NB-r),
  1544.                     DIG(x,1), alen+q);
  1545.  
  1546.       if (a->sign) {
  1547.          x->sign = 1;
  1548.          *DIG(x,0) |=
  1549.             B - (1 << r);
  1550.          compl1(DIG(x,0),
  1551.                 DIG(x,0),
  1552.                 xlen);
  1553.          }
  1554.       return mkdesc(x, dx);
  1555.       }
  1556.    }
  1557.  
  1558. /*
  1559.  *  negative if da < db
  1560.  *  zero if da == db
  1561.  *  positive if da > db
  1562.  */
  1563.  
  1564. word bigcmp(da, db)
  1565. dptr da, db;
  1566. {
  1567.    struct b_bignum *a = LrgInt(da);
  1568.    struct b_bignum *b = LrgInt(db);
  1569.    word alen, blen;
  1570.  
  1571.    if (Type(*da) == T_Lrgint && Type(*db) == T_Lrgint) {
  1572.       if (a->sign != b->sign)
  1573.          return (b->sign - a->sign);
  1574.       alen = LEN(a);
  1575.       blen = LEN(b);
  1576.       if (alen != blen)
  1577.          return (a->sign ? blen - alen : alen - blen);
  1578.  
  1579.       if (a->sign)
  1580.          return cmp1(DIG(b,0),
  1581.                      DIG(a,0),
  1582.                      alen);
  1583.       else
  1584.          return cmp1(DIG(a,0),
  1585.                      DIG(b,0),
  1586.                      alen);
  1587.       }
  1588.    else if (Type(*da) == T_Lrgint)    /* cmp(bignum, integer) */
  1589.       return bigcmpi(da, IntVal(*db));
  1590.    else                               /* cmp(integer, bignum) */
  1591.       return -bigcmpi(db, IntVal(*da));
  1592. }
  1593.  
  1594. /*
  1595.  *  ?da -> dx
  1596.  */
  1597.  
  1598. int bigrand(da, dx)
  1599. dptr da, dx;
  1600. {
  1601.    tended struct b_bignum *x, *a, *td, *tu, *tv;
  1602.    word alen = LEN(LrgInt(da));
  1603.    DIGIT *d;
  1604.    word i;
  1605.    double rval;
  1606.  
  1607.    Protect(x = alcbignum(alen), return Error);
  1608.    Protect(td = alcbignum(alen + 1), return Error);
  1609.    d = DIG(td,0);
  1610.    a = LrgInt(da);
  1611.  
  1612.    for (i = alen; i >= 0; i--) {
  1613.       rval = RandVal;
  1614.       d[i] = rval * B;
  1615.       }
  1616.  
  1617.    Protect(tu = alcbignum(alen + 2), return Error);
  1618.    Protect(tv = alcbignum(alen), return Error);
  1619.    if (div1(d, DIG(a,0),
  1620.             NULL,
  1621.             DIG(x,0),
  1622.             (word)1, alen, tu, tv) == Error)
  1623.       return Error;
  1624.    addi1(DIG(x,0),
  1625.          (word)1,
  1626.          DIG(x,0),
  1627.          alen);
  1628.    return mkdesc(x, dx);
  1629. }
  1630.  
  1631. /*
  1632.  *  da + i -> dx
  1633.  */
  1634.  
  1635. static int bigaddi(da, i, dx)
  1636. dptr da, dx;
  1637. word i;
  1638. {
  1639.    tended struct b_bignum *a;
  1640.    struct b_bignum *x;
  1641.    word alen;
  1642.  
  1643.    if (i < 0 && i > MinLong)
  1644.       return bigsubi(da, -i, dx);
  1645.    else if (i < 0 || i >= B ) {
  1646.       struct descrip td;
  1647.       char tdigits[INTBIGBLK];
  1648.  
  1649.       itobig(i, (struct b_bignum *)tdigits, &td);
  1650.       return bigadd(da, &td, dx);
  1651.       }
  1652.    else {
  1653.       alen = LEN(LrgInt(da));
  1654.       a = LrgInt(da);
  1655.       if (a->sign) {
  1656.      Protect(x = alcbignum(alen), return Error);
  1657.          subi1(DIG(a,0),
  1658.                i,
  1659.                DIG(x,0),
  1660.                alen);
  1661.          }
  1662.       else {
  1663.          Protect(x = alcbignum(alen + 1), return Error);
  1664.          *DIG(x,0) =
  1665.             addi1(DIG(a,0),
  1666.                   i,
  1667.                   DIG(x,1),
  1668.                   alen);
  1669.          }
  1670.       x->sign = a->sign;
  1671.       return mkdesc(x, dx);
  1672.       }
  1673. }
  1674.  
  1675. /*
  1676.  *  da - i -> dx
  1677.  */
  1678.  
  1679. static int bigsubi(da, i, dx)
  1680. dptr da, dx;
  1681. word i;
  1682. {
  1683.    tended struct b_bignum *a;
  1684.    struct b_bignum *x;
  1685.    word alen;
  1686.  
  1687.    if (i < 0 && i > MinLong)
  1688.       return bigaddi(da, -i, dx);
  1689.    else if (i < 0 || i >= B) {
  1690.       struct descrip td;
  1691.       char tdigits[INTBIGBLK];
  1692.  
  1693.       itobig(i, (struct b_bignum *)tdigits, &td);
  1694.       return bigsub(da, &td, dx);
  1695.       }
  1696.    else {
  1697.       alen = LEN(LrgInt(da));
  1698.       a = LrgInt(da);
  1699.       if (a->sign) {
  1700.          Protect(x = alcbignum(alen + 1), return Error);
  1701.          *DIG(x,0) =
  1702.             addi1(DIG(a,0),
  1703.                   i,
  1704.                   DIG(x,1),
  1705.                   alen);
  1706.          }
  1707.       else {
  1708.          Protect(x = alcbignum(alen), return Error);
  1709.          subi1(DIG(a,0),
  1710.                i,
  1711.                DIG(x,0),
  1712.                alen);
  1713.          }
  1714.       x->sign = a->sign;
  1715.       return mkdesc(x, dx);
  1716.       }
  1717. }
  1718.  
  1719. /*
  1720.  *  da * i -> dx
  1721.  */
  1722.  
  1723. static int bigmuli(da, i, dx)
  1724. dptr da, dx;
  1725. word i;
  1726. {
  1727.    tended struct b_bignum *a;
  1728.    struct b_bignum *x;
  1729.    word alen;
  1730.  
  1731.    if (i <= -B || i >= B) {
  1732.       struct descrip td;
  1733.       char tdigits[INTBIGBLK];
  1734.  
  1735.       itobig(i, (struct b_bignum *)tdigits, &td);
  1736.       return bigmul(da, &td, dx);
  1737.       }
  1738.    else {
  1739.       alen = LEN(LrgInt(da));
  1740.       a = LrgInt(da);
  1741.       Protect(x = alcbignum(alen + 1), return Error);
  1742.       if (i >= 0)
  1743.          x->sign = a->sign;
  1744.       else {
  1745.          x->sign = 1 ^ a->sign;
  1746.          i = -i;
  1747.          }
  1748.       *DIG(x,0) =
  1749.          muli1(DIG(a,0),
  1750.                i, 0,
  1751.                DIG(x,1),
  1752.                alen);
  1753.       return mkdesc(x, dx);
  1754.       }
  1755. }
  1756.  
  1757. /*
  1758.  *  da / i -> dx
  1759.  */
  1760.  
  1761. static int bigdivi(da, i, dx)
  1762. dptr da, dx;
  1763. word i;
  1764. {
  1765.    tended struct b_bignum *a;
  1766.    struct b_bignum *x;
  1767.    word alen;
  1768.  
  1769.    if (i <= -B || i >= B) {
  1770.       struct descrip td;
  1771.       char tdigits[INTBIGBLK];
  1772.  
  1773.       itobig(i, (struct b_bignum *)tdigits, &td);
  1774.       return bigdiv(da, &td, dx);
  1775.       }
  1776.    else {
  1777.       alen = LEN(LrgInt(da));
  1778.       a = LrgInt(da);
  1779.       Protect(x = alcbignum(alen), return Error);
  1780.       if (i >= 0)
  1781.          x->sign = a->sign;
  1782.       else {
  1783.          x->sign = 1 ^ a->sign;
  1784.          i = -i;
  1785.          }
  1786.       divi1(DIG(a,0),
  1787.             i,
  1788.             DIG(x,0),
  1789.             alen);
  1790.       return mkdesc(x, dx);
  1791.       }
  1792. }
  1793.  
  1794. /*
  1795.  *  da % i -> dx
  1796.  */
  1797.  
  1798. static int bigmodi(da, i, dx)
  1799. dptr da, dx;
  1800. word i;
  1801. {
  1802.    tended struct b_bignum *a, *temp;
  1803.    word alen;
  1804.    word x;
  1805.  
  1806.    if (i <= -B || i >= B) {
  1807.       struct descrip td;
  1808.       char tdigits[INTBIGBLK];
  1809.  
  1810.       itobig(i, (struct b_bignum *)tdigits, &td);
  1811.       return bigmod(da, &td, dx);
  1812.       }
  1813.    else {
  1814.       alen = LEN(LrgInt(da));
  1815.       a = LrgInt(da);
  1816.       temp = a;            /* avoid trash pointer */
  1817.       Protect(temp = alcbignum(alen), return Error);
  1818.       x = divi1(DIG(a,0),
  1819.                 Abs(i),
  1820.                 DIG(temp,0),
  1821.                 alen);
  1822.       if (a->sign)
  1823.      x = -x;
  1824.       MakeInt(x, dx);
  1825.       return Succeeded;
  1826.       }
  1827. }
  1828.  
  1829. /*
  1830.  *  da ^ i -> dx
  1831.  */
  1832.  
  1833. static int bigpowi(da, i, dx)
  1834. dptr da, dx;
  1835. word i;
  1836. {
  1837.    int n = WordBits;
  1838.  
  1839.    if (i > 0) {
  1840.       /* scan bits left to right.  skip leading 1. */
  1841.       while (--n >= 0)
  1842.          if (i & ((word)1 << n))
  1843.         break;
  1844.       /* then, for each zero, square the partial result;
  1845.          for each one, square it and multiply it by a */
  1846.       *dx = *da;
  1847.       while (--n >= 0) {
  1848.          if (bigmul(dx, dx, dx) == Error)
  1849.         return Error;
  1850.          if (i & ((word)1 << n))
  1851.             if (bigmul(dx, da, dx) == Error)
  1852.            return Error;
  1853.          }
  1854.       }
  1855.    else if (i == 0) {
  1856.       MakeInt(1, dx);
  1857.       }
  1858.    else {
  1859.       MakeInt(0, dx);
  1860.       }
  1861.    return Succeeded;
  1862. }
  1863.  
  1864. /*
  1865.  *  a ^ i -> dx
  1866.  */
  1867.  
  1868. static int bigpowii(a, i, dx)
  1869. word a, i;
  1870. dptr dx;
  1871. {
  1872.    word x, y;
  1873.    int n = WordBits;
  1874.    int isbig = 0;
  1875.  
  1876.    if (a == 0 || i <= 0) {              /* special cases */
  1877.       if (a == 0 && i <= 0)             /* 0 ^ negative -> error */
  1878.          ReturnErrNum(204,Error);
  1879.       if (i == 0) {
  1880.          MakeInt(1, dx);
  1881.          return Succeeded;
  1882.          }
  1883.       if (a == -1) {                    /* -1 ^ [odd,even] -> [-1,+1] */
  1884.          if (!(i & 1))
  1885.         a = 1;
  1886.          }
  1887.       else if (a != 1) {                /* 1 ^ any -> 1 */
  1888.          a = 0;
  1889.          }                   /* others ^ negative -> 0 */
  1890.       MakeInt(a, dx);
  1891.       }
  1892.    else {
  1893.       struct descrip td;
  1894.       char tdigits[INTBIGBLK];
  1895.  
  1896.       /* scan bits left to right.  skip leading 1. */
  1897.       while (--n >= 0)
  1898.          if (i & ((word)1 << n))
  1899.         break;
  1900.       /* then, for each zero, square the partial result;
  1901.          for each one, square it and multiply it by a */
  1902.       x = a;
  1903.       while (--n >= 0) {
  1904.          if (isbig) {
  1905.             if (bigmul(dx, dx, dx) == Error)
  1906.            return Error;
  1907.         }
  1908.          else {
  1909.             y = mul(x, x);
  1910.             if (!over_flow)
  1911.                x = y;
  1912.             else {
  1913.                itobig(x, (struct b_bignum *)tdigits, &td);
  1914.                if (bigmul(&td, &td, dx) == Error)
  1915.               return Error;
  1916.                isbig = (Type(*dx) == T_Lrgint);
  1917.                }
  1918.             }
  1919.          if (i & ((word)1 << n)) {
  1920.             if (isbig) {
  1921.                if (bigmuli(dx, a, dx) == Error)
  1922.           return Error;
  1923.            }
  1924.             else {
  1925.                y = mul(x, a);
  1926.                if (!over_flow)
  1927.                   x = y;
  1928.                else {
  1929.                   itobig(x, (struct b_bignum *)tdigits, &td);
  1930.                   if (bigmuli(&td, a, dx) == Error)
  1931.           return Error;
  1932.                   isbig = (Type(*dx) == T_Lrgint);
  1933.                   }
  1934.                }
  1935.             }
  1936.          }
  1937.       if (!isbig) {
  1938.      MakeInt(x, dx);
  1939.      }
  1940.       }
  1941.    return Succeeded;
  1942. }
  1943.  
  1944. /*
  1945.  *  negative if da < i
  1946.  *  zero if da == i
  1947.  *  positive if da > i
  1948.  */
  1949.  
  1950. static word bigcmpi(da, i)
  1951. dptr da;
  1952. word i;
  1953. {
  1954.    struct b_bignum *a = LrgInt(da);
  1955.    word alen = LEN(a);
  1956.  
  1957.    if (i > -B && i < B) {
  1958.       if (i >= 0)
  1959.          if (a->sign)
  1960.         return -1;
  1961.          else
  1962.         return cmpi1(DIG(a,0),
  1963.                      i, alen);
  1964.       else
  1965.          if (a->sign)
  1966.         return -cmpi1(DIG(a,0),
  1967.                       -i, alen);
  1968.          else
  1969.         return 1;
  1970.       }
  1971.    else {
  1972.       struct descrip td;
  1973.       char tdigits[INTBIGBLK];
  1974.  
  1975.       itobig(i, (struct b_bignum *)tdigits, &td);
  1976.       return bigcmp(da, &td);
  1977.       }
  1978. }
  1979.  
  1980.  
  1981. /* These are all straight out of Knuth vol. 2, Sec. 4.3.1. */
  1982.  
  1983. /*
  1984.  *  (u,n) + (v,n) -> (w,n)
  1985.  *
  1986.  *  returns carry, 0 or 1
  1987.  */
  1988.  
  1989. static DIGIT add1(u, v, w, n)
  1990. DIGIT *u, *v, *w;
  1991. word n;
  1992. {
  1993.    uword dig, carry;
  1994.    word i;
  1995.  
  1996.    carry = 0;
  1997.    for (i = n; --i >= 0; ) {
  1998.       dig = (uword)u[i] + v[i] + carry;
  1999.       w[i] = lo(dig);
  2000.       carry = hi(dig);
  2001.       }
  2002.    return carry;
  2003. }
  2004.  
  2005. /*
  2006.  *  (u,n) - (v,n) -> (w,n)
  2007.  *
  2008.  *  returns carry, 0 or -1
  2009.  */
  2010.  
  2011. static word sub1(u, v, w, n)
  2012. DIGIT *u, *v, *w;
  2013. word n;
  2014. {
  2015.    uword dig, carry;
  2016.    word i;
  2017.  
  2018.    carry = 0;
  2019.    for (i = n; --i >= 0; ) {
  2020.       dig = (uword)u[i] - v[i] + carry;
  2021.       w[i] = lo(dig);
  2022.       carry = signed_hi(dig);
  2023.       }
  2024.    return carry;
  2025. }
  2026.  
  2027. /*
  2028.  *  (u,n) * (v,m) -> (w,m+n)
  2029.  */
  2030.  
  2031. static void mul1(u, v, w, n, m)
  2032. DIGIT *u, *v, *w;
  2033. word n, m;
  2034. {
  2035.    word i, j;
  2036.    uword dig, carry;
  2037.  
  2038.    bdzero(&w[m], n);
  2039.  
  2040.    for (j = m; --j >= 0; ) {
  2041.       carry = 0;
  2042.       for (i = n; --i >= 0; ) {
  2043.          dig = (uword)u[i] * v[j] + w[i+j+1] + carry;
  2044.          w[i+j+1] = lo(dig);
  2045.          carry = hi(dig);
  2046.          }
  2047.       w[j] = carry;
  2048.       }
  2049. }
  2050.  
  2051. /*
  2052.  *  (a,m+n) / (b,n) -> (q,m+1) (r,n)
  2053.  *
  2054.  *  if q or r is NULL, the quotient or remainder is discarded
  2055.  */
  2056.  
  2057. static int div1(a, b, q, r, m, n, tu, tv)
  2058. DIGIT *a, *b, *q, *r;
  2059. word m, n;
  2060. struct b_bignum *tu, *tv;
  2061. {
  2062.    uword qhat, rhat;
  2063.    uword dig, carry;
  2064.    DIGIT *u, *v;
  2065.    word d;
  2066.    word i, j;
  2067.  
  2068.    u = DIG(tu,0);
  2069.    v = DIG(tv,0);
  2070.  
  2071.    /* D1 */
  2072.    for (d = 0; d < NB; d++)
  2073.       if (b[0] & (1 << (NB - 1 - d)))
  2074.          break;
  2075.  
  2076.    u[0] = shifti1(a, d, (DIGIT)0, &u[1], m+n);
  2077.    shifti1(b, d, (DIGIT)0, v, n);
  2078.  
  2079.    /* D2, D7 */
  2080.    for (j = 0; j <= m; j++) {
  2081.       /* D3 */
  2082.       if (u[j] == v[0]) {
  2083.          qhat = B - 1;
  2084.          rhat = (uword)v[0] + u[j+1];
  2085.          }
  2086.       else {
  2087.          uword numerator = dbl(u[j], u[j+1]);
  2088.          qhat = numerator / (uword)v[0];
  2089.          rhat = numerator % (uword)v[0];
  2090.          }
  2091.  
  2092.       while (rhat < (uword)B && qhat * (uword)v[1] > (uword)dbl(rhat, u[j+2])) {
  2093.          qhat -= 1;
  2094.          rhat += v[0];
  2095.          }
  2096.  
  2097.       /* D4 */
  2098.       carry = 0;
  2099.       for (i = n; i > 0; i--) {
  2100.          dig = u[i+j] - v[i-1] * qhat + carry;       /* -BSQ+B .. B-1 */
  2101.          u[i+j] = lo(dig);
  2102.          if ((uword)dig < (uword)B)
  2103.             carry = hi(dig);
  2104.          else carry = hi(dig) | -B;
  2105.          }
  2106.       carry = (word)(carry + u[j]) < 0;
  2107.  
  2108.       /* D5 */
  2109.       if (q)
  2110.      q[j] = qhat;
  2111.  
  2112.       /* D6 */
  2113.       if (carry) {
  2114.          if (q)
  2115.         q[j] -= 1;
  2116.          carry = 0;
  2117.          for (i = n; i > 0; i--) {
  2118.             dig = (uword)u[i+j] + v[i-1] + carry;
  2119.             u[i+j] = lo(dig);
  2120.             carry = hi(dig);
  2121.             }
  2122.          }
  2123.       }
  2124.  
  2125.    if (r) {
  2126.       if (d == 0)
  2127.          shifti1(&u[m+1], (word)d, (DIGIT)0, r, n);
  2128.       else
  2129.          r[0] = shifti1(&u[m+1], (word)(NB - d), u[m+n]>>d, &r[1], n - 1);
  2130.       }
  2131.    return Succeeded;
  2132. }
  2133.  
  2134. /*
  2135.  *  - (u,n) -> (w,n)
  2136.  *
  2137.  */
  2138.  
  2139. static void compl1(u, w, n)
  2140. DIGIT *u, *w;
  2141. word n;
  2142. {
  2143.    uword dig, carry = 0;
  2144.    word i;
  2145.  
  2146.    for (i = n; --i >= 0; ) {
  2147.       dig = carry - u[i];
  2148.       w[i] = lo(dig);
  2149.       carry = signed_hi(dig);
  2150.       }
  2151. }
  2152.  
  2153. /*
  2154.  *  (u,n) : (v,n)
  2155.  */
  2156.  
  2157. static word cmp1(u, v, n)
  2158. DIGIT *u, *v;
  2159. word n;
  2160. {
  2161.    word i;
  2162.  
  2163.    for (i = 0; i < n; i++)
  2164.       if (u[i] != v[i])
  2165.          return u[i] > v[i] ? 1 : -1;
  2166.    return 0;
  2167. }
  2168.  
  2169. /*
  2170.  *  (u,n) + k -> (w,n)
  2171.  *
  2172.  *  k in 0 .. B-1
  2173.  *  returns carry, 0 or 1
  2174.  */
  2175.  
  2176. static DIGIT addi1(u, k, w, n)
  2177. DIGIT *u, *w;
  2178. word k;
  2179. word n;
  2180. {
  2181.    uword dig, carry;
  2182.    word i;
  2183.  
  2184.    carry = k;
  2185.    for (i = n; --i >= 0; ) {
  2186.       dig = (uword)u[i] + carry;
  2187.       w[i] = lo(dig);
  2188.       carry = hi(dig);
  2189.       }
  2190.    return carry;
  2191. }
  2192.  
  2193. /*
  2194.  *  (u,n) - k -> (w,n)
  2195.  *
  2196.  *  k in 0 .. B-1
  2197.  *  u must be greater than k
  2198.  */
  2199.  
  2200. static void subi1(u, k, w, n)
  2201. DIGIT *u, *w;
  2202. word k;
  2203. word n;
  2204. {
  2205.    uword dig, carry;
  2206.    word i;
  2207.  
  2208.    carry = -k;
  2209.    for (i = n; --i >= 0; ) {
  2210.       dig = (uword)u[i] + carry;
  2211.       w[i] = lo(dig);
  2212.       carry = signed_hi(dig);
  2213.       }
  2214. }
  2215.  
  2216. /*
  2217.  *  (u,n) * k + c -> (w,n)
  2218.  *
  2219.  *  k in 0 .. B-1
  2220.  *  returns carry, 0 .. B-1
  2221.  */
  2222.  
  2223. static DIGIT muli1(u, k, c, w, n)
  2224. DIGIT *u, *w;
  2225. word k;
  2226. int c;
  2227. word n;
  2228. {
  2229.    uword dig, carry;
  2230.    word i;
  2231.  
  2232.    carry = c;
  2233.    for (i = n; --i >= 0; ) {
  2234.       dig = (uword)k * u[i] + carry;
  2235.       w[i] = lo(dig);
  2236.       carry = hi(dig);
  2237.       }
  2238.    return carry;
  2239. }
  2240.  
  2241. /*
  2242.  *  (u,n) / k -> (w,n)
  2243.  *
  2244.  *  k in 0 .. B-1
  2245.  *  returns remainder, 0 .. B-1
  2246.  */
  2247.  
  2248. static DIGIT divi1(u, k, w, n)
  2249. DIGIT *u, *w;
  2250. word k;
  2251. word n;
  2252. {
  2253.    uword dig, remain;
  2254.    word i;
  2255.  
  2256.    remain = 0;
  2257.    for (i = 0; i < n; i++) {
  2258.       dig = dbl(remain, u[i]);
  2259.       w[i] = dig / k;
  2260.       remain = dig % k;
  2261.       }
  2262.    return remain;
  2263. }
  2264.  
  2265. /*
  2266.  *  ((u,n) << k) + c -> (w,n)
  2267.  *
  2268.  *  k in 0 .. NB-1
  2269.  *  c in 0 .. B-1
  2270.  *  returns carry, 0 .. B-1
  2271.  */
  2272.  
  2273. static DIGIT shifti1(u, k, c, w, n)
  2274. DIGIT *u, c, *w;
  2275. word k;
  2276. word n;
  2277. {
  2278.    uword dig;
  2279.    word i;
  2280.  
  2281.    if (k == 0) {
  2282.       bdcopy(u, w, n);
  2283.       return 0;
  2284.       }
  2285.  
  2286.    for (i = n; --i >= 0; ) {
  2287.       dig = ((uword)u[i] << k) + c;
  2288.       w[i] = lo(dig);
  2289.       c = hi(dig);
  2290.       }
  2291.    return c;
  2292. }
  2293.  
  2294. /*
  2295.  *  (u,n) : k
  2296.  *
  2297.  *  k in 0 .. B-1
  2298.  */
  2299.  
  2300. static word cmpi1(u, k, n)
  2301. DIGIT *u;
  2302. word k;
  2303. word n;
  2304. {
  2305.    word i;
  2306.  
  2307.    for (i = 0; i < n-1; i++)
  2308.       if (u[i])
  2309.      return 1;
  2310.    if (u[n - 1] == (DIGIT)k)
  2311.       return 0;
  2312.    return u[n - 1] > (DIGIT)k ? 1 : -1;
  2313. }
  2314.  
  2315. #else                    /* LargeInts */
  2316. static char x;            /* prevent empty module */
  2317. #endif                    /* LargeInts */
  2318.