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

##### Functions

    app_comp          X
    app_diff          X
    app_exp           X
    app_floor         X
    app_log           X
    app_neg           X
    app_power         X
    app_prod          X
    app_quot          X
    app_sum           X
    mk_approx         X

##### Variables

    app_0             X
    app_0_buf         X

       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

app_comp defined in line 215; used 3 times

- in /usr/src/new/PORT/B/src/bint/b1num.c line 183-188(2), 194

app_diff defined in line 79; used 1 times

- in /usr/src/new/PORT/B/src/bint/b1num.h line 135

app_exp defined in line 137; used 3 times

- in line 200
- in /usr/src/new/PORT/B/src/bint/b1fun.c line 590
- in /usr/src/new/PORT/B/src/bint/b1num.h line 140

app_floor defined in line 234; used 2 times

- in /usr/src/new/PORT/B/src/bint/b1fun.c line 284
- in /usr/src/new/PORT/B/src/bint/b1num.h line 144

app_log defined in line 116; used 2 times

- in /usr/src/new/PORT/B/src/bint/b1fun.c line 597
- in /usr/src/new/PORT/B/src/bint/b1num.h line 141

app_neg defined in line 94; used 6 times

- in line 85
- in /usr/src/new/PORT/B/src/bint/b1fun.c line 145, 157, 227, 250
- in /usr/src/new/PORT/B/src/bint/b1num.h line 138

app_power defined in line 185; used 3 times

- in /usr/src/new/PORT/B/src/bint/b1fun.c line 248, 267
- in /usr/src/new/PORT/B/src/bint/b1num.h line 142

app_prod defined in line 98; used 2 times

- in line 199
- in /usr/src/new/PORT/B/src/bint/b1num.h line 136

app_quot defined in line 102; used 1 times

- in /usr/src/new/PORT/B/src/bint/b1num.h line 137

app_sum defined in line 67; used 1 times

- in /usr/src/new/PORT/B/src/bint/b1num.h line 134

mk_approx defined in line 40; used 25 times

- in line 72-75(2), 86-99(4), 107, 122, 144, 152, 167-169(2), 193-198(2), 209-211(2)
- in /usr/src/new/PORT/B/src/bint/b1fun.c line 470, 543-544(2), 572, 585, 618
- in /usr/src/new/PORT/B/src/bint/b1nuC.c line 88
- in /usr/src/new/PORT/B/src/bint/b1num.c line 303
- in /usr/src/new/PORT/B/src/bint/b1num.h line 132

#### Defined variables

app_0 defined in line 34; used 5 times

- in line 46, 52, 120, 193-194(2)

app_0_buf defined in line 31; used 1 times

- in line 34

Last modified: 1985-08-27 Generated: 2016-12-26 | Generated by src2html V0.67 page hit count: 912 |