2.11BSD#431 - /usr/src/new/PORT/B/src/bsmall/B1num.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */
   2: /* $Header: B1num.c,v 1.1 84/06/28 00:48:56 timo Exp $ */
   3: 
   4: /* B numbers, small version */
   5: 
   6:     /*
   7: 	 * THIS VERSION SHOULD ONLY BE USED IF
   8: 	 * THE SYSTEM IS TOO LARGE OTHERWISE.
   9: 	 * IT USES FLOATING POINT ARITHMETIC FOR EXACT NUMBERS
  10: 	 * INSTEAD OF ARBITRARY LENGTH RATIONAL ARITHMETIC.
  11: 	 */
  12: 
  13: #include "b.h"
  14: #include "b0con.h"
  15: #include "b1obj.h"
  16: #include "b2syn.h" /* for def of Keymark() only */
  17: #include "B1num.h"
  18: 
  19: value numerator(v) value v; {
  20:     Checknum(v);
  21:     if (!Exact(v)) error("*/ on approximate number");
  22:     return mk_int(Numerator(v));
  23: }
  24: 
  25: value denominator(v) value v; {
  26:     Checknum(v);
  27:     if (!Exact(v)) error("/* on approximate number"); /* */
  28:     return mk_int(Denominator(v));
  29: }
  30: 
  31: double numval(v) value v; {
  32:     Checknum(v);
  33:     return Numval(v);
  34: }
  35: 
  36: checkint(v) value v; {
  37:     Checknum(v);
  38:     if (Denominator(v) != One) error("number not an integer");
  39: }
  40: 
  41: bool large(v) value v; {
  42:     checkint(v);
  43:     if (Numerator(v) < -Maxint || Numerator(v) > Maxint) return Yes;
  44:     return No;
  45: }
  46: 
  47: int intval(v) value v; {
  48:     checkint(v);
  49:     return (int)Numerator(v);
  50: }
  51: 
  52: intlet propintlet(i) int i; {
  53:     if (i < -Maxintlet || i > Maxintlet)
  54:         error("exceedingly large integer");
  55:     return i;
  56: }
  57: 
  58: integer gcd(i, j) integer i, j; {
  59:     integer k;
  60:     if (i == Zero && j == Zero) syserr("gcd(0, 0)");
  61:     if (i != floor(i) || j != floor(j))
  62:         syserr("gcd called with non-integer");
  63:     if (i < Zero) i= -i; if (j < Zero) j= -j;
  64:     if (i < j) {
  65:         k= i; i= j; j= k;
  66:     }
  67:     while (j >= One) {
  68:         k= i-j*floor(i/j);
  69:         i= j; j= k;
  70:     }
  71:     if (j != Zero) error(
  72:         "arithmetic overflow while simplifying exact number");
  73:     if (i != floor(i)) syserr("gcd returns non-integer");
  74:     return i;
  75: }
  76: 
  77: value b_zero, b_one, b_minus_one, zero, one;
  78: 
  79: value mk_exact(p, q, len) register integer p, q; intlet len; {
  80:     value v; integer d;
  81:     if (q == One && len ==0) {
  82:         if (p == Zero) return copy(b_zero);
  83:         if (p == One)  return copy(b_one);
  84:         if (p == -One) return copy(b_minus_one);
  85:     }
  86:     v= grab_num(len);
  87:     if (q == One) {
  88:         Numerator(v)= p; Denominator(v)= q;
  89:         return v;
  90:     }
  91:     if (q == Zero) error("attempt to make exact number with denominator 0");
  92:     if (q < Zero) {p= -p; q= -q;}
  93:     d= (q == One ? One : p == One ? One : gcd(p, q));
  94:     Numerator(v)= p/d; Denominator(v)= q/d;
  95:     return v;
  96: }
  97: 
  98: bool integral(v) value v; {
  99:     return Integral(v);
 100: }
 101: 
 102: value mk_integer(p) int p; {
 103:     return mk_exact((integer)p, One, 0);
 104: }
 105: 
 106: value mk_int(p) integer p; {
 107:     return mk_exact(p, One, 0);
 108: }
 109: 
 110: value mk_approx(x) register double x; {
 111:     value v= grab_num(0);
 112:     Approxval(v)= x; Denominator(v)= Zero;
 113:     return v;
 114: }
 115: 
 116: initnum() {
 117:     b_zero= grab_num(0);
 118:         Numerator(b_zero)= Zero; Denominator(b_zero)= One;
 119:     b_one= grab_num(0);
 120:         Numerator(b_one)= One; Denominator(b_one)= One;
 121:     b_minus_one= grab_num(0);
 122:         Numerator(b_minus_one)= -One; Denominator(b_minus_one)= One;
 123:     zero= mk_integer(0);
 124:     one= mk_integer(1);
 125: }
 126: 
 127: value approximate(v) value v; {
 128:     if (!Exact(v)) return copy(v);
 129:     return mk_approx(Numerator(v)/Denominator(v));
 130: }
 131: 
 132: numcomp(v, w) value v, w; {
 133:     double vv= Numval(v), ww= Numval(w);
 134:     if (vv < ww) return -1;
 135:     if (vv > ww) return  1;
 136:     if (Exact(v) && Exact(w)) return 0;
 137:     if (Exact(v)) return -1; /* 1 < 1E0 */
 138:     if (Exact(w)) return  1; /* 1E0 > 1 */
 139:     return 0;
 140: }
 141: 
 142: double numhash(v) value v; {
 143:     number *n= (number *)Ats(v);
 144:     return .123*n->p + .777*n->q;
 145: }
 146: 
 147: #define CONVBUFSIZ 100
 148: char convbuf[CONVBUFSIZ];
 149: 
 150: string convnum(v) value v; {
 151:     double x; string bp; bool prec_loss= No;
 152:     Checknum(v);
 153:     x= Numval(v);
 154:  conv:  if (!prec_loss && Exact(v) && fabs(x) <= LONG &&
 155:         fabs(Numerator(v)) < BIG && fabs(Denominator(v)) < BIG) {
 156:         intlet len= 0 < Length(v) && Length(v) <= MAXNUMDIG ? Length(v) : 0;
 157:         intlet dcnt, sigcnt; bool sig;
 158:         if (Denominator(v) != One) {
 159:             intlet k; double p= 1.0, q;
 160:             prec_loss= Yes;
 161:             for (k= 1; k < MAXNUMDIG; k++) {
 162:                 p*= 10.0;
 163:                 q= p/Denominator(v);
 164:                 if (k >= len && q == floor(q)) {
 165:                     prec_loss= No;
 166:                     break;
 167:                 }
 168:             }
 169:             len= k;
 170:         }
 171:     convex: sprintf(convbuf, "%.*f", len, x);
 172:         dcnt= sigcnt= 0; sig= No;
 173:         for (bp= convbuf; *bp != '\0'; bp++)
 174:             if ('0' <= *bp && *bp <= '9') {
 175:                 dcnt++;
 176:                 if (*bp != '0') sig= Yes;
 177:                 if (sig) sigcnt++;
 178:             }
 179:         if (sigcnt < MINNUMDIG && prec_loss) goto conv;
 180:         if (dcnt > MAXNUMDIG) {
 181:             if (len <= 0) syserr("conversion error 1");
 182:             if (Denominator(v) == One) len= 0;
 183:             else len-= dcnt-MAXNUMDIG;
 184:             if (len < 0) syserr("conversion error 2");
 185:             goto convex;
 186:         }
 187:     } else { /*approx etc*/
 188:         sprintf(convbuf, "%.*e", MAXNUMDIG-5, x);
 189:         for (bp= convbuf; *bp != '\0'; bp++)
 190:         if (*bp == 'e') {
 191:             *bp= 'E';
 192:             break;
 193:         }
 194:     }
 195:     return convbuf;
 196: }
 197: 
 198: value numconst(tx, q) txptr tx, q; {
 199:     bool dig= No; double ex= 0, ap= 1; intlet ndap, len= 0;
 200:     while (tx < q && '0' <= *tx && *tx <= '9') {
 201:         dig= Yes;
 202:         ex= 10*ex+(*tx++ - '0');
 203:     }
 204:     if (tx < q && *tx == '.') {
 205:         tx++; ndap= 0;
 206:         while (tx < q && '0' <= *tx && *tx <= '9') {
 207:             dig= Yes; ndap++;
 208:             len= *tx == '0' ? ndap : 0;
 209:             ex= 10*ex+(*tx++ - '0'); ap*= 10;
 210:         }
 211:         if (!dig) syserr("numconst[1]");
 212:     }
 213:     if (tx < q && *tx == 'E') {
 214:         intlet sign= 1; double expo= 0;
 215:         tx++;
 216:         if (!('0' <= *tx && *tx <= '9') && Keymark(*tx)) {
 217:             tx--;
 218:             goto exact;
 219:         }
 220:         if (!dig) ex= 1;
 221:         if (tx < q && (*tx == '+' || *tx == '-'))
 222:             if (*tx++ == '-') sign= -1;
 223:         dig= No;
 224:         while (tx < q && '0' <= *tx && *tx <= '9') {
 225:             dig= Yes;
 226:             expo= 10*expo+(*tx++ - '0');
 227:         }
 228:         if (!dig) syserr("numconst[2]");
 229:         return mk_approx(ex/ap*exp(sign*expo*log(10.0)));
 230:     }
 231: exact:  return mk_exact(ex, ap, len);
 232: }
 233: 
 234: printnum(f1, v) FILE *f1; value v; {
 235:     FILE *f= f1 ? f1 : stdout;
 236:     if (!Exact(v) || Denominator(v) == One) {
 237:         if (!Exact(v))
 238:             fputc('~', f);
 239:         fputs(convnum(v), f);
 240:     }
 241:     else {
 242:         value w = numerator(v);
 243:         fputs(convnum(w), f);
 244:         release(w);
 245:         fputc('/', f);
 246:         w = denominator(v);
 247:         fputs(convnum(w), f);
 248:         release(w);
 249:     }
 250:     if (!f1) fputc('\n', f); /* Flush buffer for sdb */
 251: }

Defined functions

Defined variables

Defined macros