2.11BSD#431 - /usr/src/new/PORT/B/src/bint/b1num.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1num.c,v 1.4 85/08/22 16:51:59 timo Exp $
   5: */
   6: 
   7: /* B numbers, basic external interface */
   8: 
   9: #include "b.h"
  10: #include "b0con.h"
  11: #include "b1obj.h"
  12: #include "b1num.h"
  13: 
  14: /*
  15:  * This file contains operations on numbers that are not predefined
  16:  * B functions (the latter are defined in `bfun.c').
  17:  * This includes conversion of numeric B values to C `int's and
  18:  * `double's, (numval() and intval()),
  19:  * but also utilities for comparing numeric values and hashing a numeric
  20:  * value to something usable as initialization for the random generator
  21:  * without chance of overflow (so numval(v) is unusable).
  22:  * It is also possible to build numbers of all types using mk_integer,
  23:  * mk_exact (or mk_rat) and mk_approx.  Note the rather irregular
  24:  * (historical!) argument structure for these: mk_approx has a
  25:  * `double' argument, where mk_exact and mk_rat have two values
  26:  * which must be `integer' (not `int's!) as arguments.
  27:  * The random number generator used by the DRAW and CHOOSE statements
  28:  * is also in this file.
  29:  */
  30: 
  31: /*
  32:  * ival is used internally by intval() and large().
  33:  * It converts an integer to double precision, yielding
  34:  * a huge value when overflow occurs (but giving no error).
  35:  */
  36: 
  37: Hidden double ival(u) integer u; {
  38:     double x = 0;
  39:     register int i;
  40: 
  41:     if (IsSmallInt(u)) return SmallIntVal(u);
  42:     for (i = Length(u)-1; i >= 0; --i) {
  43:         if (x >= Maxreal/BASE)
  44:             return Msd(u) < 0 ? -Maxreal : Maxreal;
  45:         x = x*BASE + Digit(u, i);
  46:     }
  47: 
  48:     return x;
  49: }
  50: 
  51: 
  52: /* Determine if a value may be converted to an int */
  53: 
  54: Visible bool large(v) value v; {
  55:     double r;
  56:     if (!Is_number(v) || !integral(v)) {
  57:         error(MESS(1300, "number not an integer"));
  58:         return No;
  59:     }
  60:     if (Rational(v)) v= (value) Numerator((rational)v);
  61:     r= ival((integer)v);
  62:     if (r > Maxint) return Yes;
  63:     return No;
  64: }
  65: 
  66: /* return the C `int' value of a B numeric value, if it exists. */
  67: 
  68: Visible int intval(v) value v; {
  69:     /* v must be an Integral number or a Rational with Denominator==1
  70: 	    which may result from n round x [via mk_exact]!. */
  71:     double i;
  72:     if (IsSmallInt(v)) i= SmallIntVal(v);
  73:     else {
  74:         if (!Is_number(v)) syserr(MESS(1301, "intval on non-number"));
  75:         if (!integral(v)) {
  76:             error(MESS(1302, "number not an integer"));
  77:             return 0;
  78:         }
  79:         if (Rational(v)) v= (value) Numerator((rational)v);
  80:         i= ival((integer)v);
  81:     }
  82:     if (i > Maxint || i < -Maxint) {
  83:         error(MESS(1303, "exceedingly large integer"));
  84:         return 0;
  85:     }
  86:     return (int) i;
  87: }
  88: 
  89: 
  90: /* convert an int to an intlet */
  91: 
  92: Visible intlet propintlet(i) int i; {
  93:     if (i > Maxintlet || i < -Maxintlet) {
  94:         error(MESS(1304, "exceedingly large integer"));
  95:         return 0;
  96:     }
  97:     return i;
  98: }
  99: 
 100: 
 101: /*
 102:  * determine if a number is integer
 103:  */
 104: 
 105: Visible bool integral(v) value v; {
 106:     if (Integral(v) || (Rational(v) && Denominator((rational)v) == int_1))
 107:         return Yes;
 108:     else return No;
 109: }
 110: 
 111: /*
 112:  * mk_exact makes an exact number out of two integers.
 113:  * The third argument is the desired number of digits after the decimal
 114:  * point when the number is printed (see comments in `bfun.c' for
 115:  * `round' and code in `bnuC.c').
 116:  * This printing size (if positive) is hidden in the otherwise nearly
 117:  * unused * 'Size' field of the value struct
 118:  * (cf. definition of Rational(v) etc.).
 119:  */
 120: 
 121: Visible value mk_exact(p, q, len) integer p, q; int len; {
 122:     rational r = mk_rat(p, q, len);
 123: 
 124:     if (Denominator(r) == int_1 && len <= 0) {
 125:         integer i = (integer) Copy(Numerator(r));
 126:         release((value) r);
 127:         return (value) i;
 128:     }
 129: 
 130:     return (value) r;
 131: }
 132: 
 133: #define uSMALL 1
 134: #define uINT 2
 135: #define uRAT 4
 136: #define uAPP 8
 137: #define vSMALL 16
 138: #define vINT 32
 139: #define vRAT 64
 140: #define vAPP 128
 141: 
 142: Visible relation numcomp(u, v) value u, v; {
 143:     int tu, tv; relation s;
 144: 
 145:     if (IsSmallInt(u)) tu = uSMALL;
 146:     else if (Length(u) >= 0) tu = uINT;
 147:     else if (Length(u) <= -2) tu = uRAT;
 148:     else tu = uAPP;
 149:     if (IsSmallInt(v)) tv = vSMALL;
 150:     else if (Length(v) >= 0) tv = vINT;
 151:     else if (Length(v) <= -2) tv = vRAT;
 152:     else tv = vAPP;
 153: 
 154:     switch (tu|tv) {
 155: 
 156:     case uSMALL|vSMALL: return SmallIntVal(u) - SmallIntVal(v);
 157: 
 158:     case uSMALL|vINT:
 159:     case uINT|vSMALL:
 160:     case uINT|vINT: return int_comp((integer)u, (integer)v);
 161: 
 162:     case uSMALL|vRAT:
 163:     case uINT|vRAT:
 164:         u = (value) mk_rat((integer)u, int_1, 0);
 165:         s = rat_comp((rational)u, (rational)v);
 166:         release(u);
 167:         return s;
 168: 
 169:     case uRAT|vRAT:
 170:         return rat_comp((rational)u, (rational)v);
 171: 
 172:     case uRAT|vSMALL:
 173:     case uRAT|vINT:
 174:         v = (value) mk_rat((integer)v, int_1, 0);
 175:         s = rat_comp((rational)u, (rational)v);
 176:         release(v);
 177:         return s;
 178: 
 179:     case uSMALL|vAPP:
 180:     case uINT|vAPP:
 181:     case uRAT|vAPP:
 182:         u = approximate(u);
 183:         s = app_comp((real)u, (real)v);
 184:         release(u);
 185:         return s == 0 ? -1 : s; /* u < ~u = v */
 186: 
 187:     case uAPP|vAPP:
 188:         return app_comp((real)u, (real)v);
 189: 
 190:     case uAPP|vSMALL:
 191:     case uAPP|vINT:
 192:     case uAPP|vRAT:
 193:         v = approximate(v);
 194:         s = app_comp((real)u, (real)v);
 195:         release(v);
 196:         return s == 0 ? 1 : s; /* u = ~v > v */
 197: 
 198:     default: syserr(MESS(1305, "num_comp")); /* NOTREACHED */
 199: 
 200:     }
 201: }
 202: 
 203: 
 204: /*
 205:  * Deliver 10**n, where n is a (maybe negative!) C integer.
 206:  * The result is a value (integer or rational, actually).
 207:  */
 208: 
 209: Visible value tento(n) int n; {
 210:     if (n < 0) {
 211:         integer i= int_tento(-n);
 212:         value v= (value) mk_exact(int_1, i, 0);
 213:         release((value) i);
 214:         return v;
 215:     }
 216:     return (value) int_tento(n);
 217: }
 218: 
 219: 
 220: /*
 221:  * numval returns the numerical value of any numeric B value
 222:  * as a C `double'.
 223:  */
 224: 
 225: Visible double numval(u) value u; {
 226:     double expo, frac;
 227: 
 228:     if (!Is_number(u)) {
 229:         error(MESS(1306, "value not a number"));
 230:         return 0.0;
 231:     }
 232:     u = approximate(u);
 233:     expo = Expo((real)u), frac = Frac((real)u);
 234:     release(u);
 235:     if (expo > Maxexpo) {
 236:         error(MESS(1307, "approximate number too large to be handled"));
 237:         return 0.0;
 238:     }
 239:     if(expo < Minexpo) return 0.0;
 240:     return ldexp(frac, (int)expo);
 241: }
 242: 
 243: 
 244: /*
 245:  * Random numbers
 246:  */
 247: 
 248: 
 249: /*
 250:  * numhash produces a `double' number that depends on the value of
 251:  * v, useful for initializing the random generator.
 252:  * Needs rewriting, so that for instance numhash(n) never equals n.
 253:  * The magic numbers here are chosen at random.
 254:  */
 255: 
 256: Visible double numhash(v) value v; {
 257:     if (Integral(v)) {
 258:         double d = 0;
 259:         register int i;
 260: 
 261:         if (IsSmallInt(v)) return SmallIntVal(v);
 262:         for (i = Length(v) - 1; i >= 0; --i) {
 263:             d *= 2;
 264:             d += Digit((integer)v, i);
 265:         }
 266: 
 267:         return d;
 268:     }
 269: 
 270:     if (Rational(v))
 271:         return .777 * numhash((value) Numerator((rational)v)) +
 272:                .123 * numhash((value) Denominator((rational)v));
 273: 
 274:     return numval(v); /* Which fails for HUGE reals.  Never mind. */
 275: }
 276: 
 277: 
 278: /* Initialize the random generator */
 279: 
 280: double lastran;
 281: 
 282: Hidden Procedure setran (seed) double seed; {
 283:     double x;
 284: 
 285:     x = seed >= 0 ? seed : -seed;
 286:     while (x >= 1) x /= 10;
 287:     lastran = floor(67108864.0*x);
 288: }
 289: 
 290: Visible Procedure set_random(v) value v; {
 291:     setran((double) hash(v));
 292: }
 293: 
 294: 
 295: /* Return a random number in [0, 1). */
 296: 
 297: Visible value random() {
 298:     double p;
 299: 
 300:     p = 26353589.0 * lastran + 1;
 301:     lastran = p - 67108864.0*floor(p/67108864.0);
 302: 
 303:     return (value) mk_approx(lastran / 67108864.0, 0.0);
 304: }
 305: 
 306: Visible Procedure initnum() {
 307:     rat_init();
 308:     setran((double) SEED);
 309: }
 310: 
 311: Visible Procedure endnum() {
 312:     endrat();
 313: }

Defined functions

Defined variables

Defined macros