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