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

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1nuI.c,v 1.4 85/08/22 16:51:13 timo Exp $
   5: */
   6: 
   7: /* Multi-precision integer arithmetic */
   8: 
   9: #include "b.h"
  10: #include "b1obj.h"
  11: #include "b1num.h"
  12: #include "b0con.h"
  13: #include "b3err.h"
  14: 
  15: /*
  16:  * Number representation:
  17:  * ======================
  18:  *
  19:  * (Think of BASE = 10 for ordinary decimal notation.)
  20:  * A number is a sequence of N "digits" b1, b2, ..., bN
  21:  * where each bi is in {0..BASE-1}, except for negative numbers,
  22:  * where bN = -1.
  23:  * The number represented by b1, ..., bN is
  24:  *      b1*BASE**(N-1) + b2*BASE**(N-2) + ... + bN .
  25:  * The base BASE is chosen so that multiplication of two positive
  26:  * integers up to BASE-1 can be multiplied exactly using double
  27:  * precision floating point arithmetic.
  28:  * Also it must be possible to add two long integers between
  29:  * -BASE and +BASE (exclusive), giving a result between -2BASE and
  30:  * +2BASE.
  31:  * BASE must be even (so we can easily decide whether the whole
  32:  * number is even), and positive (to avoid all kinds of other trouble).
  33:  * Presently, it is restricted to a power of 10 by the I/O-conversion
  34:  * routines (file "b1nuC.c").
  35:  *
  36:  * Canonical representation:
  37:  * bN is never zero (for the number zero itself, N is zero).
  38:  * If bN is -1, b[N-1] is never BASE-1 .
  39:  * All operands are assumed te be in canonical representation.
  40:  * Routine "int_canon" brings a number in canonical representation.
  41:  *
  42:  * Mapping to C objects:
  43:  * A "digit" is an integer of type "digit", probably an "int".
  44:  * A number is represented as a "B-integer", i.e. something
  45:  * of type "integer" (which is actually a pointer to some struct).
  46:  * The number of digits N is extracted through the macro Length(v).
  47:  * The i-th digit is extracted through the macro Digit(v,N-i).
  48:  * (So in C, we count in a backwards direction from 0 ... n-1 !)
  49:  * A number is created through a call to grab_num(N), which sets
  50:  * N zero digits (thus not in canonical form!).
  51:  */
  52: 
  53: 
  54: /*
  55:  * Bring an integer into canonical form.
  56:  * Make a SmallInt if at all possible.
  57:  * NB: Work done by int_canon is duplicated by mk_integer for optimization;
  58:  *     if the strategy here changes, look at mk_integer, too!
  59:  */
  60: 
  61: Visible integer int_canon(v) integer v; {
  62:     register int i;
  63: 
  64:     if (IsSmallInt(v)) return v;
  65: 
  66:     for (i = Length(v) - 1; i >= 0 && Digit(v,i) == 0; --i)
  67:         ;
  68: 
  69:     if (i < 0) {
  70:         release((value) v);
  71:         return int_0;
  72:     }
  73: 
  74:     if (i == 0) {
  75:         digit dig = Digit(v,0);
  76:         release((value) v);
  77:         return (integer) MkSmallInt(dig);
  78:     }
  79: 
  80:     if (i > 0 && Digit(v,i) == -1) {
  81:         while (i > 0 && Digit(v, i-1) == BASE-1) --i;
  82:         if (i == 0) {
  83:             release((value) v);
  84:             return (integer) MkSmallInt(-1);
  85:         }
  86:         if (i == 1) {
  87:             digit dig = Digit(v,0) - BASE;
  88:             release((value) v);
  89:             return (integer) MkSmallInt(dig);
  90:         }
  91:         Digit(v,i) = -1;
  92:     }
  93: 
  94:     if (i+1 < Length(v)) return (integer) regrab_num((value) v, i+1);
  95: 
  96:     return v;
  97: }
  98: 
  99: 
 100: /* General add/subtract subroutine */
 101: 
 102: typedef double twodigit; /* Might be long on 16 bit machines */
 103:     /* Should be in b0con.h */
 104: 
 105: Hidden twodigit fmodulo(x, y) twodigit x, y; {
 106:     return x - y * (twodigit) floor((double)x / (double)y);
 107: }
 108: 
 109: Visible Procedure dig_gadd(to, nto, from, nfrom, ffactor)
 110:     digit *to, *from; intlet nto, nfrom; digit ffactor; {
 111:     twodigit carry= 0;
 112:     twodigit factor= ffactor;
 113:     digit save;
 114: 
 115:     nto -= nfrom;
 116:     if (nto < 0)
 117:         syserr(MESS(1000, "dig_gadd: nto < nfrom"));
 118:     for (; nfrom > 0; ++to, ++from, --nfrom) {
 119:         carry += *to + *from * factor;
 120:         *to= save= fmodulo(carry, (twodigit)BASE);
 121:         carry= (carry-save) / BASE;
 122:     }
 123:     for (; nto > 0; ++to, --nto) {
 124:         if (carry == 0)
 125:             return;
 126:         carry += *to;
 127:         *to= save= fmodulo(carry, (twodigit)BASE);
 128:         carry= (carry-save) / BASE;
 129:     }
 130:     if (carry != 0)
 131:         to[-1] += carry*BASE; /* Assume it's -1 */
 132: }
 133: 
 134: 
 135: /* Sum or difference of two integers */
 136: /* Should have its own version of dig-gadd without double precision */
 137: 
 138: Visible integer int_gadd(v, w, factor) integer v, w; intlet factor; {
 139:     struct integer vv, ww;
 140:     integer s;
 141:     int len, lenv, i;
 142: 
 143:     FreezeSmallInt(v, vv);
 144:     FreezeSmallInt(w, ww);
 145:     lenv= len= Length(v);
 146:     if (Length(w) > len)
 147:         len= Length(w);
 148:     ++len;
 149:     s= (integer) grab_num(len);
 150:     for (i= 0; i < lenv; ++i)
 151:         Digit(s, i)= Digit(v, i);
 152:     for (; i < len; ++i)
 153:         Digit(s, i)= 0;
 154:     dig_gadd(&Digit(s, 0), len, &Digit(w, 0), Length(w), (digit)factor);
 155:     return int_canon(s);
 156: }
 157: 
 158: 
 159: /* Product of two integers */
 160: 
 161: Visible integer int_prod(v, w) integer v, w; {
 162:     int i;
 163:     integer a;
 164:     struct integer vv, ww;
 165: 
 166:     if (v == int_0 || w == int_0) return int_0;
 167:     if (v == int_1) return (integer) Copy(w);
 168:     if (w == int_1) return (integer) Copy(v);
 169: 
 170:     FreezeSmallInt(v, vv);
 171:     FreezeSmallInt(w, ww);
 172: 
 173:     a = (integer) grab_num(Length(v) + Length(w));
 174: 
 175:     for (i= Length(a)-1; i >= 0; --i)
 176:         Digit(a, i)= 0;
 177:     for (i = 0; i < Length(v) && !interrupted; ++i)
 178:         dig_gadd(&Digit(a, i), Length(w)+1, &Digit(w, 0), Length(w),
 179:             Digit(v, i));
 180: 
 181:     return int_canon(a);
 182: }
 183: 
 184: 
 185: /* Compare two integers */
 186: 
 187: Visible relation int_comp(v, w) integer v, w; {
 188:     int sv, sw;
 189:     register int i;
 190:     struct integer vv, ww;
 191: 
 192:     /* 1. Compare pointers and equal SmallInts */
 193:     if (v == w) return 0;
 194: 
 195:     /* 1a. Handle SmallInts */
 196:     if (IsSmallInt(v) && IsSmallInt(w))
 197:         return SmallIntVal(v) - SmallIntVal(w);
 198:     FreezeSmallInt(v, vv);
 199:     FreezeSmallInt(w, ww);
 200: 
 201:     /* 2. Extract signs */
 202:     sv = Length(v)==0 ? 0 : Digit(v,Length(v)-1)<0 ? -1 : 1;
 203:     sw = Length(w)==0 ? 0 : Digit(w,Length(w)-1)<0 ? -1 : 1;
 204: 
 205:     /* 3. Compare signs */
 206:     if (sv != sw) return (sv>sw) - (sv<sw);
 207: 
 208:     /* 4. Compare sizes */
 209:     if (Length(v) != Length(w))
 210:         return sv * ( (Length(v)>Length(w)) - (Length(v)<Length(w)) );
 211: 
 212:     /* 5. Compare individual digits */
 213:     for (i = Length(v)-1; i >= 0 && Digit(v,i) == Digit(w,i); --i)
 214:         ;
 215: 
 216:     /* 6. All digits equal? */
 217:     if (i < 0) return 0;  /* Yes */
 218: 
 219:     /* 7. Compare leftmost different digits */
 220:     if (Digit(v,i) < Digit(w,i)) return -1;
 221: 
 222:     return 1;
 223: }
 224: 
 225: 
 226: /* Construct an integer out of a floating point number */
 227: 
 228: #define GRAN 8  /* Granularity used when requesting more storage */
 229:         /* MOVE TO MEM! */
 230: Visible integer mk_int(x) double x; {
 231:     register integer a;
 232:     integer b;
 233:     register int i, j;
 234:     int negate;
 235: 
 236:     if (MinSmallInt <= x && x <= MaxSmallInt)
 237:         return (integer) MkSmallInt((int)x);
 238: 
 239:     a = (integer) grab_num(1);
 240:     negate = x < 0 ? 1 : 0;
 241:     if (negate) x = -x;
 242: 
 243:     for (i = 0; x != 0; ++i) {
 244:         double z = floor(x/BASE);
 245:         digit save = Modulo((digit)(x-z*BASE), BASE);
 246:         if (i >= Length(a)) {
 247:             a = (integer) regrab_num((value) a, Length(a)+GRAN);
 248:             for (j = Length(a)-1; j > i; --j)
 249:                 Digit(a,j) = 0; /* clear higher digits */
 250:         }
 251:         Digit(a,i) = save;
 252:         x = floor((x-save)/BASE);
 253:     }
 254: 
 255:     if (negate) {
 256:         b = int_neg(a);
 257:         release((value) a);
 258:         return b;
 259:     }
 260: 
 261:     return int_canon(a);
 262: }
 263: 
 264: /* Construct an integer out of a C int.  Like mk_int, but optimized. */
 265: 
 266: Visible value mk_integer(x) int x; {
 267:     if (MinSmallInt <= x && x <= MaxSmallInt) return MkSmallInt(x);
 268:     return (value) mk_int((double)x);
 269: }
 270: 
 271: 
 272: /* Efficiently compute 10**n as a B integer, where n is a C int >= 0 */
 273: 
 274: Visible integer int_tento(n) int n; {
 275:     integer i;
 276:     digit msd = 1;
 277:     if (n < 0) syserr(MESS(1001, "int_tento(-n)"));
 278:     if (n < tenlogBASE) {
 279:         while (n != 0) msd *= 10, --n;
 280:         return (integer) MkSmallInt(msd);
 281:     }
 282:     i = (integer) grab_num(1 + (int)(n/tenlogBASE));
 283:     n %= tenlogBASE;
 284:     while (n != 0) msd *= 10, --n;
 285:     Digit(i, Length(i)-1) = msd;
 286:     return i;
 287: }
 288: 
 289: #ifdef NOT_USED
 290: /* Approximate ceiling(10 log abs(u/v)), as C int.
 291:    It only works for v > 0, u, v both integers.
 292:    The result may be one too large or too small */
 293: 
 294: Visible int scale(u, v) integer u, v; {
 295:     int s;
 296:     double z;
 297:     struct integer uu, vv;
 298: 
 299:     if (Msd(v) <= 0) syserr(MESS(1002, "scale(u,v<=0)"));
 300:     if (u == int_0) return 0; /* `Don't care' case */
 301:     FreezeSmallInt(u, uu);
 302:     FreezeSmallInt(v, vv);
 303:     s = (Length(u) - Length(v)) * tenlogBASE;
 304:     if (Digit(u, Length(u)-1) >= 0) z = Digit(u, Length(u)-1);
 305:     else {
 306:         s -= tenlogBASE;
 307:         if (Length(u) == 1) z = 1;
 308:         else z = BASE - Digit(u, Length(u)-2);
 309:     }
 310:     z /= Digit(v, Length(v)-1);
 311:     while (z >= 10) z /= 10, ++s;
 312:     while (z < 1) z *= 10, --s;
 313:     return s;
 314: }
 315: #endif NOT_USED

Defined functions

Defined typedef's

Defined macros