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

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
   2: 
   3: /*
   4:   $Header: b1nuA.c,v 1.4 85/08/22 16:50:25 timo Exp $
   5: */
   6: 
   7: /* Approximate arithmetic */
   8: 
   9: #include "b.h"
  10: #include "b0con.h"
  11: #include "b1obj.h"
  12: #include "b1num.h"
  13: #include "b3err.h" /* For still_ok */
  14: 
  15: /*
  16: For various reasons, on some machines (notably the VAX), the range
  17: of the exponent is too small (ca. 1.7E38), and we cope with this by
  18: adding a second word which holds the exponent.
  19: However, on other machines (notably the IBM PC), the range is sufficient
  20: (ca. 1E300), and here we try to save as much code as possible by not
  21: doing our own exponent handling.  (To be fair, we also don't check
  22: certain error conditions, to save more code.)
  23: The difference is made by #defining EXT_RANGE (in b1num.h), meaning we
  24: have to EXTend the RANGE of the exponent.
  25: */
  26: 
  27: #ifdef EXT_RANGE
  28: Hidden struct real app_0_buf = {Num, 1, -1, 0.0, -(double)Maxint};
  29:     /* Exponent must be less than any realistic exponent! */
  30: #else !EXT_RANGE
  31: Hidden struct real app_0_buf = {Num, 1, -1, 0.0};
  32: #endif !EXT_RANGE
  33: 
  34: Visible real app_0 = &app_0_buf;
  35: 
  36: /*
  37:  * Build an approximate number.
  38:  */
  39: 
  40: Visible real mk_approx(frac, expo) double frac, expo; {
  41:     real u;
  42: #ifdef EXT_RANGE
  43:     expint e;
  44:     if (frac != 0) frac = frexp(frac, &e), expo += e;
  45:     if (frac == 0.5) frac = 1, --expo; /* Assert 0.5 < frac <= 1 */
  46:     if (frac == 0 || expo < -BIG) return (real) Copy(app_0);
  47:     if (expo > BIG) {
  48:         error(MESS(700, "approximate number too large"));
  49:         expo = BIG;
  50:     }
  51: #else !EXT_RANGE
  52:     if (frac == 0.0) return (real) Copy(app_0);
  53:     frac= ldexp(frac, (int)expo);
  54: #endif EXT_RANGE
  55:     u = (real) grab_num(-1);
  56:     Frac(u) = frac;
  57: #ifdef EXT_RANGE
  58:     Expo(u) = expo;
  59: #endif EXT_RANGE
  60:     return u;
  61: }
  62: 
  63: /*
  64:  * Approximate arithmetic.
  65:  */
  66: 
  67: Visible real app_sum(u, v) real u, v; {
  68: #ifdef EXT_RANGE
  69:     real w;
  70:     if (Expo(u) < Expo(v)) w = u, u = v, v = w;
  71:     if (Expo(v) - Expo(u) < Minexpo) return (real) Copy(u);
  72:     return mk_approx(Frac(u) + ldexp(Frac(v), (int)(Expo(v) - Expo(u))),
  73:         Expo(u));
  74: #else !EXT_RANGE
  75:     return mk_approx(Frac(u) + Frac(v), 0.0);
  76: #endif !EXT_RANGE
  77: }
  78: 
  79: Visible real app_diff(u, v) real u, v; {
  80: #ifdef EXT_RANGE
  81:     real w;
  82:     int sign = 1;
  83:     if (Expo(u) < Expo(v)) w = u, u = v, v = w, sign = -1;
  84:     if (Expo(v) - Expo(u) < Minexpo)
  85:         return sign < 0 ? app_neg(u) : (real) Copy(u);
  86:     return mk_approx(
  87:         sign * (Frac(u) - ldexp(Frac(v), (int)(Expo(v) - Expo(u)))),
  88:         Expo(u));
  89: #else !EXT_RANGE
  90:     return mk_approx(Frac(u) - Frac(v), 0.0);
  91: #endif !EXT_RANGE
  92: }
  93: 
  94: Visible real app_neg(u) real u; {
  95:     return mk_approx(-Frac(u), Expo(u));
  96: }
  97: 
  98: Visible real app_prod(u, v) real u, v; {
  99:     return mk_approx(Frac(u) * Frac(v), Expo(u) + Expo(v));
 100: }
 101: 
 102: Visible real app_quot(u, v) real u, v; {
 103:     if (Frac(v) == 0.0) {
 104:         error(MESS(701, "in u/v, v is zero"));
 105:         return (real) Copy(u);
 106:     }
 107:     return mk_approx(Frac(u) / Frac(v), Expo(u) - Expo(v));
 108: }
 109: 
 110: /*
 111: 	YIELD log"(frac, expo):
 112: 		CHECK frac > 0
 113: 		RETURN normalize"(expo*logtwo + log(frac), 0)
 114: */
 115: 
 116: Visible real app_log(v) real v; {
 117:     double frac = Frac(v), expo = Expo(v);
 118:     if (frac <= 0) {
 119:         error(MESS(702, "in log x, x is <= 0"));
 120:         return (real) Copy(app_0);
 121:     }
 122:     return mk_approx(expo*logtwo + log(frac), 0.0);
 123: }
 124: 
 125: /*
 126: 	YIELD exp"(frac, expo):
 127: 		IF expo < minexpo: RETURN zero"
 128: 		WHILE expo < 0: PUT frac/2, expo+1 IN frac, expo
 129: 		PUT exp frac IN f
 130: 		PUT normalize"(f, 0) IN f, e
 131: 		WHILE expo > 0:
 132: 			PUT (f, e) prod" (f, e) IN f, e
 133: 			PUT expo-1 IN expo
 134: 		RETURN f, e
 135: */
 136: 
 137: Visible real app_exp(v) real v; {
 138: #ifdef EXT_RANGE
 139:     expint ei;
 140:     double frac = Frac(v), expo = Expo(v), new_expo;
 141:     static double canexp;
 142:     if (!canexp) canexp = floor(log(log(Maxreal)) / logtwo);
 143:     if (expo <= canexp) {
 144:         if (expo < Minexpo) return mk_approx(1.0, 0.0);
 145:         frac = ldexp(frac, (int)expo);
 146:         expo = 0;
 147:     }
 148:     else if (expo >= Maxexpo) {
 149:         /* Definitely too big (the real boundary is much smaller
 150: 		   but here we are in danger of overflowing new_expo
 151: 		   in the loop below) */
 152:         return mk_approx(1.0, Maxreal); /* Force an error! */
 153:     }
 154:     else {
 155:         frac = ldexp(frac, (int)canexp);
 156:         expo -= canexp;
 157:     }
 158:     frac = exp(frac);
 159:     new_expo = 0;
 160:     while (expo > 0 && frac != 0) {
 161:         frac = frexp(frac, &ei);
 162:         new_expo += ei;
 163:         frac *= frac;
 164:         new_expo += new_expo;
 165:         --expo;
 166:     }
 167:     return mk_approx(frac, new_expo);
 168: #else !EXT_RANGE
 169:     return mk_approx(exp(Frac(v)), 0.0);
 170: #endif !EXT_RANGE
 171: }
 172: 
 173: /*
 174: 	YIELD (frac, expo) power" v":
 175: 		\   (f*2**e)**v =
 176: 		\ = f**v * 2**(e*v) =
 177: 		\ = f**v * 2**((e*v) mod 1) * 2**floor(e*v) .
 178: 		PUT exp"(v" prod" normalize"(log frac, 0)) IN temp1" \ = f**v
 179: 		PUT expo*numval(v") IN ev \ = e*v
 180: 		PUT exp(logtwo * (ev - floor ev))) IN temp2 \ = 2**(ev mod 1)
 181: 		PUT temp1" IN f, e
 182: 		RETURN normalize"(f*temp2, e + floor ev)
 183: */
 184: 
 185: Visible real app_power(u, v) real u, v; {
 186:     double frac = Frac(u);
 187: #ifdef EXT_RANGE
 188:     real logfrac, vlogfrac, result;
 189:     double expo = Expo(u), rest;
 190: #endif !EXT_RANGE
 191:     if (frac <= 0) {
 192:         if (frac < 0) error(MESS(703, "in 0**v, v is negative"));
 193:         if (v == app_0) return mk_approx(1.0, 0.0); /* 0**0 = 1 */
 194:         return (real) Copy(app_0); /* 0**x = 0 */
 195:     }
 196: #ifdef EXT_RANGE
 197:     frac *= 2, expo -= 1; /* Renormalize to 1 < frac <= 2, so log frac > 0 */
 198:     logfrac = mk_approx(log(frac), 0.0);
 199:     vlogfrac = app_prod(v, logfrac);
 200:     result = app_exp(vlogfrac);
 201:     /* But what if result overflows but expo is very negative??? */
 202:     if (still_ok) {
 203:         expo *= numval((value)v);
 204:         rest = expo - floor(expo);
 205:         frac = Frac(result) * exp(logtwo*rest);
 206:         expo = Expo(result) + floor(expo);
 207:     }
 208:     release((value)logfrac), release((value)vlogfrac), release((value)result);
 209:     return mk_approx(frac, expo);
 210: #else !EXT_RANGE
 211:     return mk_approx(exp(log(frac) * Frac(v)), 0.0);
 212: #endif !EXT_RANGE
 213: }
 214: 
 215: Visible int app_comp(u, v) real u, v; {
 216:     double xu, xv;
 217: #ifdef EXT_RANGE
 218:     double eu, ev;
 219: #endif EXT_RANGE
 220:     if (u == v) return 0;
 221:     xu = Frac(u), xv = Frac(v);
 222: #ifdef EXT_RANGE
 223:     if (xu*xv > 0) {
 224:         eu = Expo(u), ev = Expo(v);
 225:         if (eu < ev) return xu < 0 ? 1 : -1;
 226:         if (eu > ev) return xu < 0 ? -1 : 1;
 227:     }
 228: #endif EXT_RANGE
 229:     if (xu < xv) return -1;
 230:     if (xu > xv) return 1;
 231:     return 0;
 232: }
 233: 
 234: Visible value app_floor(u) real u; {
 235: #ifdef EXT_RANGE
 236:     integer v, w;
 237:     value twotow, result;
 238:     if (Expo(u) <= Dblbits)
 239:         return (value)
 240:             mk_int(floor(ldexp(Frac(u),
 241:                 (int)(Expo(u) < 0 ? -1 : Expo(u)))));
 242:     v = mk_int(ldexp(Frac(u), Dblbits));
 243:     w = mk_int(Expo(u) - Dblbits);
 244:     twotow = power((value)int_2, (value)w);
 245:     result = prod((value)v, twotow);
 246:     release((value) v), release((value) w), release(twotow);
 247:     if (!Integral(result)) syserr(MESS(704, "app_floor: result not integral"));
 248:     return result;
 249: #else !EXT_RANGE
 250:     return (value) mk_int(floor(Frac(u)));
 251: #endif !EXT_RANGE
 252: }

Defined functions

Defined variables