```   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
app_diff defined in line 79; used 1 times
app_exp defined in line 137; used 3 times
app_neg defined in line 94; used 6 times
app_prod defined in line 98; used 2 times
app_quot defined in line 102; used 1 times
app_sum defined in line 67; used 1 times

#### Defined variables

app_0 defined in line 34; used 5 times
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: 2682