4.3BSD - /usr/src/new/B/src/bint/b1nuQ.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1nuQ.c,v 1.4 85/08/22 16:51:40 timo Exp $
   5: */
   6: 
   7: #include "b.h"
   8: #include "b1obj.h"
   9: #include "b0con.h"
  10: #include "b1num.h"
  11: 
  12: 
  13: /* Product of integer and single "digit" */
  14: 
  15: Visible integer int1mul(v, n1) integer v; digit n1; {
  16:     integer a;
  17:     digit save, bigcarry, carry = 0;
  18:     double z, zz, n = n1;
  19:     register int i;
  20:     struct integer vv;
  21: 
  22:     FreezeSmallInt(v, vv);
  23: 
  24:     a = (integer) grab_num(Length(v)+2);
  25: 
  26:     for (i = 0; i < Length(v); ++i) {
  27:         z = Digit(v,i) * n;
  28:         bigcarry = zz = floor(z/BASE);
  29:         carry += z - zz*BASE;
  30:         Digit(a,i) = save = Modulo(carry, BASE);
  31:         carry = (carry-save)/BASE + bigcarry;
  32:     }
  33: 
  34:     Digit(a,i) = save = Modulo(carry, BASE);
  35:     Digit(a,i+1) = (carry-save)/BASE;
  36: 
  37:     return int_canon(a);
  38: }
  39: 
  40: 
  41: /* Quotient of positive integer and single "digit" > 0 */
  42: 
  43: Hidden integer int1div(v, n1, prem) integer v; digit n1, *prem; {
  44:     integer q;
  45:     double r_over_n, r = 0, n = n1;
  46:     register int i;
  47:     struct integer vv;
  48: 
  49:     FreezeSmallInt(v, vv);
  50: 
  51:     q = (integer) grab_num(Length(v));
  52:     for (i = Length(v)-1; i >= 0; --i) {
  53:         r = r*BASE + Digit(v,i);
  54:         Digit(q,i) = r_over_n = floor(r/n);
  55:         r -= r_over_n * n;
  56:     }
  57:     if (prem)
  58:         *prem = r;
  59:     return int_canon(q);
  60: }
  61: 
  62: 
  63: /* Long division routine, gives access to division algorithm. */
  64: 
  65: Visible digit int_ldiv(v1, w1, pquot, prem) integer v1, w1, *pquot, *prem; {
  66:     integer a;
  67:     int sign = 1, rel_v = 0, rel_w = 0;
  68:     digit div, rem;
  69:     struct integer vv1, ww1;
  70: 
  71:     if (w1 == int_0) syserr(MESS(1100, "zero division (int_ldiv)"));
  72: 
  73:     /* Make v, w positive */
  74:     if (Msd(v1) < 0) {
  75:         sign = -1;
  76:         ++rel_v;
  77:         v1 = int_neg(v1);
  78:     }
  79: 
  80:     if (Msd(w1) < 0) {
  81:         sign *= -1;
  82:         ++rel_w;
  83:         w1 = int_neg(w1);
  84:     }
  85: 
  86:     FreezeSmallInt(v1, vv1);
  87:     FreezeSmallInt(w1, ww1);
  88: 
  89:     div = sign;
  90: 
  91:     /* Check v << w or single-digit w */
  92:     if (Length(v1) < Length(w1)
  93:         || Length(v1) == Length(w1)
  94:             && Digit(v1, Length(v1)-1) < Digit(w1, Length(w1)-1)) {
  95:         a = int_0;
  96:         if (prem) {
  97:             if (v1 == &vv1) *prem= (integer) MkSmallInt(Digit(v1,0));
  98:             else *prem = (integer) Copy(v1);
  99:         }
 100:     }
 101:     else if (Length(w1) == 1) {
 102:         /* Single-precision division */
 103:         a = int1div(v1, Digit(w1,0), &rem);
 104:         if (prem) *prem = mk_int((double)rem);
 105:     }
 106:     else {
 107:         /* Multi-precision division */
 108:         /* Cf. Knuth II Sec. 4.3.1. Algorithm D */
 109:         /* Note that we count in the reverse direction (not easier!) */
 110: 
 111:         double z, zz;
 112:         digit carry, save, bigcarry;
 113:         double q, d = BASE/(Digit(w1, Length(w1)-1)+1);
 114:         register int i, j, k;
 115:         integer v, w;
 116:         digit vj;
 117: 
 118:         /* Normalize: make Msd(w) >= BASE/2 by multiplying
 119: 		   both v and w by d */
 120: 
 121:         v = int1mul(v1, (digit)d);
 122:             /* v is used as accumulator, must make a copy */
 123:             /* v cannot be int_1 */
 124:             /* (then it would be one of the cases above) */
 125: 
 126:         if (d == 1) w = (integer) Copy(w1);
 127:         else w = int1mul(w1, (digit)d);
 128: 
 129:         a = (integer) grab_num(Length(v1)-Length(w)+1);
 130: 
 131:         /* Division loop */
 132: 
 133:         for (j = Length(v1), k = Length(a)-1; k >= 0; --j, --k) {
 134:             vj = j >= Length(v) ? 0 : Digit(v,j);
 135: 
 136:             /* Find trial digit */
 137: 
 138:             if (vj == Digit(w, Length(w)-1)) q = BASE-1;
 139:             else q = floor( ((double)vj*BASE
 140:                     + Digit(v,j-1)) / Digit(w, Length(w)-1) );
 141: 
 142:             /* Correct trial digit */
 143: 
 144:             while (Digit(w,Length(w)-2) * q >
 145:                 ((double)vj*BASE + Digit(v,j-1)
 146:                     - q*Digit(w, Length(w)-1)) *BASE + Digit(v,j-2))
 147:                 --q;
 148: 
 149:             /* Subtract q*w from v */
 150: 
 151:             carry = 0;
 152:             for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
 153:                 z = Digit(w,i) * q;
 154:                 bigcarry = zz = floor(z/BASE);
 155:                 carry += Digit(v,i+k) - z + zz*BASE;
 156:                 Digit(v,i+k) =
 157:                     save = Modulo(carry, BASE);
 158:                 carry = (carry-save)/BASE - bigcarry;
 159:             }
 160: 
 161:             if (i+k < Length(v))
 162:                 carry += Digit(v, i+k), Digit(v, i+k) = 0;
 163: 
 164:             /* Add back necessary? */
 165: 
 166:                 /* It is very difficult to find test cases
 167: 				   where add back is necessary if BASE is large.
 168: 				   Thanks to Arjen Lenstra, we have v=n*n-1, w=n,
 169: 				   where n = 8109636009903000000 (the last six
 170: 				   digits are not important). */
 171: 
 172:             if (carry == 0)     /* No */
 173:                 Digit(a,k) = q;
 174:             else {      /* Yes, add back */
 175:                 if (carry != -1) syserr(MESS(1101, "int_ldiv internal failure"));
 176:                 Digit(a,k) = q-1;
 177:                 carry = 0;
 178:                 for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
 179:                     carry += Digit(v, i+k) + Digit(w,i);
 180:                     Digit(v,i+k) =
 181:                         save = Modulo(carry, BASE);
 182:                     carry = (carry-save)/BASE;
 183:                 }
 184:             }
 185:         }   /* End for(j) */
 186: 
 187:         if (prem) *prem = int_canon(v); /* Store remainder */
 188:         else release((value) v);
 189:         div = sign*d;   /* Store normalization factor */
 190:         release((value) w);
 191:         a = int_canon(a);
 192:     }
 193: 
 194:     if (rel_v) release((value) v1);
 195:     if (rel_w) release((value) w1);
 196: 
 197:     if (sign < 0) {
 198:         integer temp = a;
 199:         a = int_neg(a);
 200:         release((value) temp);
 201:     }
 202: 
 203:     if (pquot) *pquot = a;
 204:     else release((value) a);
 205:     return div;
 206: }
 207: 
 208: 
 209: Visible integer int_quot(v, w) integer v, w; {
 210:     integer quo;
 211:     VOID int_ldiv(v, w, &quo, (integer*)0);
 212:     return quo;
 213: }
 214: 
 215: Visible integer int_mod(v, w) integer v, w; {
 216:     integer rem;
 217:     digit div;
 218:     bool flag;
 219:     div = int_ldiv(v, w, (integer*)0, &rem); /* Rem. is always positive */
 220:     if (rem == int_0)
 221:         return rem; /* v mod w = 0 */
 222:     flag = (div < 0);
 223:     if (flag || Msd(w) < 0) div = -div;
 224:     if (div != 1) { /* Divide by div to get proper remainder back */
 225:         v = int1div(rem, div, (digit*)0);
 226:         release((value) rem);
 227:         rem = v;
 228:     }
 229:     if (flag) { /* Make same sign as w */
 230:         if (Msd(w) < 0) v = int_sum(w, rem);
 231:         else v = int_diff(w, rem);
 232:         release((value) rem);
 233:         rem = v;
 234:     }
 235:     return rem;
 236: }

Defined functions