![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to b1nuA.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint|
Return to b1nuA.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint |
1.1 root 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: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.