![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to b1fun.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 b1fun.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: b1fun.c,v 1.4 85/08/22 16:48:24 timo Exp $
5: */
6:
7: /* Functions defined on numeric values. */
8:
9: #include <errno.h> /* For EDOM and ERANGE */
10:
11: #include "b.h"
12: #include "b0con.h"
13: #include "b1obj.h"
14: #include "b1num.h"
15:
16: /*
17: * The visible routines here implement predefined B arithmetic operators,
18: * taking one or two numeric values as operands, and returning a numeric
19: * value.
20: * No type checking of operands is done: this must be done by the caller.
21: */
22:
23: typedef value (*valfun)();
24: typedef rational (*ratfun)();
25: typedef real (*appfun)();
26: typedef double (*mathfun)();
27:
28: /*
29: * For the arithmetic functions (+, -, *, /) the same action is needed:
30: * 1) if both operands are Integral, use function from int_* submodule;
31: * 2) if both are Exact, use function from rat_* submodule (after possibly
32: * converting one of them from Integral to Rational);
33: * 3) otherwise, make both approximate and use function from app_*
34: * submodule.
35: * The functions performing the appropriate action for each of the submodules
36: * are passed as parameters.
37: * Division is a slight exception, since i/j can be a rational.
38: * See `quot' below.
39: */
40:
41: Hidden value dyop(u, v, int_fun, rat_fun, app_fun)
42: value u, v;
43: valfun int_fun;
44: ratfun rat_fun;
45: appfun app_fun;
46: {
47: if (Integral(u) && Integral(v)) /* Use integral operation */
48: return (*int_fun)(u, v);
49:
50: if (Exact(u) && Exact(v)) {
51: rational u1, v1, a;
52:
53: /* Use rational operation */
54:
55: u1 = Integral(u) ? mk_rat((integer)u, int_1, 0) :
56: (rational) Copy(u);
57: v1 = Integral(v) ? mk_rat((integer)v, int_1, 0) :
58: (rational) Copy(v);
59: a = (*rat_fun)(u1, v1);
60: release((value) u1);
61: release((value) v1);
62:
63: if (Denominator(a) == int_1 && Roundsize(a) == 0) {
64: integer b = (integer) Copy(Numerator(a));
65: release((value) a);
66: return (value) b;
67: }
68:
69: return (value) a;
70: }
71:
72: /* Use approximate operation */
73:
74: {
75: real u1, v1, a;
76: u1 = Approximate(u) ? (real) Copy(u) : (real) approximate(u);
77: v1 = Approximate(v) ? (real) Copy(v) : (real) approximate(v);
78: a = (*app_fun)(u1, v1);
79: release((value) u1);
80: release((value) v1);
81:
82: return (value) a;
83: }
84: }
85:
86:
87: Hidden integer isum(u, v) integer u, v; {
88: return int_sum(u, v);
89: }
90:
91: Visible value sum(u, v) value u, v; {
92: if (IsSmallInt(u) && IsSmallInt(v))
93: return mk_integer(SmallIntVal(u) + SmallIntVal(v));
94: return dyop(u, v, (value (*)())isum, rat_sum, app_sum);
95: }
96:
97: Hidden integer idiff(u, v) integer u, v; {
98: return int_diff(u, v);
99: }
100:
101: Visible value diff(u, v) value u, v; {
102: if (IsSmallInt(u) && IsSmallInt(v))
103: return mk_integer(SmallIntVal(u) - SmallIntVal(v));
104: return dyop(u, v, (value (*)())idiff, rat_diff, app_diff);
105: }
106:
107: Visible value prod(u, v) value u, v; {
108: if (IsSmallInt(u) && IsSmallInt(v))
109: return (value)
110: mk_int((double)SmallIntVal(u) * (double)SmallIntVal(v));
111: return dyop(u, v, (value (*)())int_prod, rat_prod, app_prod);
112: }
113:
114:
115: /*
116: * We cannot use int_quot (which performs integer division with truncation).
117: * Here is the routine we need.
118: */
119:
120: Hidden value xxx_quot(u, v) integer u, v; {
121:
122: if (v == int_0) {
123: error(MESS(400, "in i/j, j is zero"));
124: return (value) Copy(u);
125: }
126:
127: return mk_exact(u, v, 0);
128: }
129:
130: Visible value quot(u, v) value u, v; {
131: return dyop(u, v, xxx_quot, rat_quot, app_quot);
132: }
133:
134:
135: /*
136: * Unary minus and abs follow the same principle but with only one operand.
137: */
138:
139: Visible value negated(u) value u; {
140: if (IsSmallInt(u)) return mk_integer(-SmallIntVal(u));
141: if (Integral(u))
142: return (value) int_neg((integer)u);
143: if (Rational(u))
144: return (value) rat_neg((rational)u);
145: return (value) app_neg((real)u);
146: }
147:
148:
149: Visible value absval(u) value u; {
150: if (Integral(u)) {
151: if (Msd((integer)u) < 0)
152: return (value) int_neg((integer)u);
153: } else if (Rational(u)) {
154: if (Msd(Numerator((rational)u)) < 0)
155: return (value) rat_neg((rational)u);
156: } else if (Approximate(u) && Frac((real)u) < 0)
157: return (value) app_neg((real)u);
158:
159: return Copy(u);
160: }
161:
162:
163: /*
164: * The remaining operators follow less similar paths and some of
165: * them contain quite subtle code.
166: */
167:
168: Visible value mod(u, v) value u, v; {
169: value p, q, m, n;
170:
171: if (v == (value)int_0 ||
172: Rational(v) && Numerator((rational)v) == int_0 ||
173: Approximate(v) && Frac((real)v) == 0) {
174: error(MESS(401, "in x mod y, y is zero"));
175: return Copy(u);
176: }
177:
178: if (Integral(u) && Integral(v))
179: return (value) int_mod((integer)u, (integer)v);
180:
181: /* Compute `u-v*floor(u/v)', as in the formal definition of `u mod v'. */
182:
183: q = quot(u, v);
184: n = floorf(q);
185: release(q);
186: p = prod(n, v);
187: release(n);
188: m = diff(u, p);
189: release(p);
190:
191: return m;
192: }
193:
194:
195: /*
196: * u**v has the most special cases of all the predefined arithmetic functions.
197: */
198:
199: Visible value power(u, v) value u, v; {
200: real ru, rv, rw;
201: if (Exact(u) && (Integral(v) ||
202: /* Next check catches for integers disguised as rationals: */
203: Rational(v) && Denominator((rational)v) == int_1)) {
204: rational a;
205: integer b = Integral(v) ? (integer)v : Numerator((rational)v);
206: /* Now b is really an integer. */
207:
208: u = Integral(u) ? (value) mk_rat((integer)u, int_1, 0) :
209: Copy(u);
210:
211: a = rat_power((rational)u, b);
212: release(u);
213: if (Denominator(a) == int_1) { /* Make integral result */
214: b = (integer) Copy(Numerator(a));
215: release((value) a);
216: return (value)b;
217: }
218: return (value)a;
219: }
220:
221: if (Exact(v)) {
222: integer vn, vd;
223: int s;
224: ru = (real) approximate(u);
225: s = (Frac(ru) > 0) - (Frac(ru) < 0);
226:
227: if (s < 0) rv = app_neg(ru), release((value)ru), ru = rv;
228: if (Integral(v)) {
229: vn = (integer)v;
230: vd = int_1;
231: } else {
232: vd = Denominator((rational)v);
233: if (s < 0 && Even(Lsd(vd)))
234: error(MESS(402, "in x**(p/q), x is negative and q is even"));
235: vn = Numerator((rational)v);
236: }
237: if (vn == int_0) {
238: release((value)ru);
239: return one;
240: }
241: if (s == 0 && Msd(vn) < 0) {
242: error(MESS(403, "in 0**y, y is negative"));
243: return (value) ru;
244: }
245: if (s < 0 && Even(Lsd(vn)))
246: s = 1;
247: rv = (real) approximate(v);
248: rw = app_power(ru, rv);
249: release((value)ru), release((value)rv);
250: if (s < 0) ru = app_neg(rw), release((value)rw), rw = ru;
251: return (value) rw;
252: }
253:
254: /* Everything else: we now know u or v is approximate */
255:
256: ru = (real) approximate(u);
257: if (Frac(ru) < 0) {
258: error(MESS(404, "in x**y, x is negative and y is not exact"));
259: return (value) ru;
260: }
261: rv = (real) approximate(v);
262: if (Frac(ru) == 0 && Frac(rv) < 0) {
263: error(MESS(405, "in 0**y, y is negative"));
264: release((value)rv);
265: return (value) ru;
266: }
267: rw = app_power(ru, rv);
268: release((value)ru), release((value)rv);
269: return (value) rw;
270: }
271:
272:
273: /*
274: * floor: for approximate numbers app_floor() is used;
275: * for integers it is a no-op; other exact numbers effectively calculate
276: * u - (u mod 1).
277: */
278:
279: Visible value floorf(u) value u; {
280: integer quo, rem, v;
281: digit div;
282:
283: if (Integral(u)) return Copy(u);
284: if (Approximate(u)) return app_floor((real)u);
285:
286: /* It is a rational number */
287:
288: div = int_ldiv(Numerator((rational)u), Denominator((rational)u),
289: &quo, &rem);
290: if (div < 0 && rem != int_0) { /* Correction for negative noninteger */
291: v = int_diff(quo, int_1);
292: release((value) quo);
293: quo = v;
294: }
295: release((value) rem);
296: return (value) quo;
297: }
298:
299:
300: /*
301: * ceiling x is defined as -floor(-x);
302: * and that's how it's implemented, except for integers.
303: */
304:
305: Visible value ceilf(u) value u; {
306: value v;
307: if (Integral(u)) return Copy(u);
308: u = negated(u);
309: v = floorf(u);
310: release(u);
311: u = negated(v);
312: release(v);
313: return u;
314: }
315:
316:
317: /*
318: * round u is defined as floor(u+0.5), which is what is done here,
319: * except for integers which are left unchanged.
320: */
321:
322: Visible value round1(u) value u; {
323: value v, w;
324:
325: if (Integral(u)) return Copy(u);
326:
327: v = sum(u, (value)rat_half);
328: w = floorf(v);
329: release(v);
330:
331: return w;
332: }
333:
334:
335: /*
336: * u round v is defined as 10**-u * round(v*10**u).
337: * A complication is that u round v is always printed with exactly u digits
338: * after the decimal point, even if this involves trailing zeros,
339: * or if v is an integer.
340: * Consequently, the result is always kept as a rational, even if it can be
341: * simplified to an integer, and the size field of the rational number
342: * (which is made negative to distinguish it from integers, and < -1 to
343: * distinguish it from approximate numbers) is used to store the number of
344: * significant digits.
345: * Thus a size of -2 means a normal rational number, and a size < -2
346: * means a rounded number to be printed with (-2 - length) digits
347: * after the decimal point. This last expression can be retrieved using
348: * the macro Roundsize(v) which should only be applied to Rational
349: * numbers.
350: */
351:
352: Visible value round2(u, v) value u, v; {
353: value w, tenp;
354: int i;
355:
356: if (!Integral(u)) {
357: error(MESS(406, "in n round x, n is not an integer"));
358: i = 0;
359: } else
360: i = propintlet(intval(u));
361:
362: tenp = tento(i);
363:
364: v = prod(v, tenp);
365: w = round1(v);
366:
367: release(v);
368: v = quot(w, tenp);
369: release(w);
370: release(tenp);
371:
372: if (i > 0) { /* Set number of digits to be printed */
373: if (propintlet(-2 - i) < -2) {
374: if (Rational(v))
375: Length(v) = -2 - i;
376: else if (Integral(v)) {
377: w = v;
378: v = mk_exact((integer) w, int_1, i);
379: release(w);
380: }
381: }
382: }
383:
384: return v;
385: }
386:
387:
388: /*
389: * sign u inspects the sign of either u, u's numerator or u's fractional part.
390: */
391:
392: Visible value signum(u) value u; {
393: int s;
394:
395: if (Exact(u)) {
396: if (Rational(u))
397: u = (value) Numerator((rational)u);
398: s = u==(value)int_0 ? 0 : Msd((integer)u) < 0 ? -1 : 1;
399: } else
400: s = Frac((real)u) > 0 ? 1 : Frac((real)u) < 0 ? -1 : 0;
401:
402: return MkSmallInt(s);
403: }
404:
405:
406: /*
407: * ~u makes an approximate number of any numerical value.
408: */
409:
410: Visible value approximate(u) value u; {
411: double x, e, expo;
412: register int i;
413: if (Approximate(u)) return Copy(u);
414: if (IsSmallInt(u))
415: x = SmallIntVal(u), expo = 0;
416: else if (Rational(u)) {
417: rational r = (rational) u;
418: integer p = Numerator(r), q;
419: double xp, xq;
420: int lp, lq;
421: if (p == int_0)
422: return Copy((value)app_0);
423: q = Denominator(r);
424: lp = IsSmallInt(p) ? 1 : Length(p);
425: lq = IsSmallInt(q) ? 1 : Length(q);
426: xp = 0, xq = 0;
427: expo = lp - lq;
428: if (IsSmallInt(p)) { xp = SmallIntVal(p); --expo; }
429: else {
430: for (i = Length(p)-1; i >= 0; --i) {
431: xp = xp*BASE + Digit(p, i);
432: --expo;
433: if (xp < -Maxreal/BASE || xp > Maxreal/BASE)
434: break;
435: }
436: }
437: if (IsSmallInt(q)) { xq = SmallIntVal(q); ++expo; }
438: else {
439: for (i = Length(q)-1; i >= 0; --i) {
440: xq = xq*BASE + Digit(q, i);
441: ++expo;
442: if (xq > Maxreal/BASE)
443: break;
444: }
445: }
446: x = xp/xq;
447: }
448: else if (Integral(u)) {
449: integer p = (integer) u;
450: if (Length(p) == 0)
451: return Copy((value)app_0);
452: x = 0;
453: expo = Length(p);
454: for (i = Length(p)-1; i >= 0; --i) {
455: x = x*BASE + Digit(p, i);
456: --expo;
457: if (x < -Maxreal/BASE || x > Maxreal/BASE)
458: break;
459: }
460: }
461: while (expo < 0 && fabs(x) > BASE/Maxreal) ++expo, x /= BASE;
462: while (expo > 0 && fabs(x) < Maxreal/BASE) --expo, x *= BASE;
463: if (expo != 0) {
464: expo = expo*twologBASE + 1;
465: e = floor(expo);
466: x *= .5 * exp((expo-e) * logtwo);
467: }
468: else
469: e = 0;
470: return (value) mk_approx(x, e);
471: }
472:
473:
474: /*
475: * numerator v returns the numerator of v, whenever v is an exact number.
476: * For integers, that is v itself.
477: */
478:
479: Visible value numerator(v) value v; {
480: if (!Exact(v)) {
481: error(MESS(407, "*/ on approximate number"));
482: return zero;
483: }
484:
485: if (Integral(v)) return Copy(v);
486:
487: return (value) Copy(Numerator((rational)v));
488: }
489:
490:
491: /*
492: * /*v returns the denominator of v, whenever v is an exact number.
493: * For integers, that is 1.
494: */
495:
496: Visible value denominator(v) value v; {
497: if (!Exact(v)) {
498: error(MESS(408, "/* on approximate number"));
499: return zero;
500: }
501:
502: if (Integral(v)) return one;
503:
504: return (value) Copy(Denominator((rational)v));
505: }
506:
507:
508: /*
509: * u root v is defined as v**(1/u), where u is usually but need not be
510: * an integer.
511: */
512:
513: Visible value root2(u, v) value u, v; {
514: if (u == (value)int_0 ||
515: Rational(u) && Numerator((rational)u) == int_0 ||
516: Approximate(u) && Frac((real)u) == 0) {
517: error(MESS(409, "in n root x, n is zero"));
518: v = Copy(v);
519: } else {
520: u = quot((value)int_1, u);
521: v = power(v, u);
522: release(u);
523: }
524:
525: return v;
526: }
527:
528: /* root x is computed more exactly than n root x, by doing
529: one iteration step extra. This ~guarantees root(n**2) = n. */
530:
531: Visible value root1(v) value v; {
532: value r = root2((value)int_2, v);
533: value v_over_r, theirsum, result;
534: if (Approximate(r) && Frac((real)r) == 0.0) return (value)r;
535: v_over_r = quot(v, r);
536: theirsum = sum(r, v_over_r), release(r), release(v_over_r);
537: result = quot(theirsum, (value)int_2), release(theirsum);
538: return result;
539: }
540:
541: /* The rest of the mathematical functions */
542:
543: Visible value pi() { return (value) mk_approx(3.141592653589793238463, 0.0); }
544: Visible value e() { return (value) mk_approx(2.718281828459045235360, 0.0); }
545:
546: Hidden value trig(v, fun, zeroflag) value v; mathfun fun; bool zeroflag; {
547: real w = (real) approximate(v);
548: double expo = Expo(w), frac = Frac(w), x, result;
549: extern int errno;
550: if (expo <= Minexpo/2) {
551: if (zeroflag) return (value) w; /* sin small x = x, etc. */
552: frac = 0, expo = 0;
553: }
554: release((value)w);
555: if (expo > Maxexpo) errno = EDOM;
556: else {
557: x = ldexp(frac, (int)expo);
558: if (x >= Maxtrig || x <= -Maxtrig) errno = EDOM;
559: else {
560: errno = 0;
561: result = (*fun)(x);
562: }
563: }
564: if (errno != 0) {
565: if (errno == ERANGE)
566: error(MESS(410, "the result is too large to be representable"));
567: else if (errno == EDOM)
568: error(MESS(411, "x is too large to compute a meaningful result"));
569: else error(MESS(412, "math library error"));
570: return Copy((value)app_0);
571: }
572: return (value) mk_approx(result, 0.0);
573: }
574:
575: Visible value sin1(v) value v; { return trig(v, sin, Yes); }
576: Visible value cos1(v) value v; { return trig(v, cos, No); }
577: Visible value tan1(v) value v; { return trig(v, tan, Yes); }
578:
579: Visible value atn1(v) value v; {
580: real w = (real) approximate(v);
581: double expo = Expo(w), frac = Frac(w);
582: if (expo <= Minexpo + 2) return (value) w; /* atan of small x = x */
583: release((value)w);
584: if (expo > Maxexpo) expo = Maxexpo;
585: return (value) mk_approx(atan(ldexp(frac, (int)expo)), 0.0);
586: }
587:
588: Visible value exp1(v) value v; {
589: real w = (real) approximate(v);
590: real x = app_exp(w);
591: release((value)w);
592: return (value) x;
593: }
594:
595: Visible value log1(v) value v; {
596: real w = (real) approximate(v);
597: real x = app_log(w);
598: release((value)w);
599: return (value) x;
600: }
601:
602: Visible value log2(u, v) value u, v;{
603: value w;
604: u = log1(u);
605: v = log1(v);
606: w = quot(v, u);
607: release(u), release(v);
608: return w;
609: }
610:
611: Visible value atn2(u, v) value u, v; {
612: real ru = (real) approximate(u), rv = (real) approximate(v);
613: double uexpo = Expo(ru), ufrac = Frac(ru);
614: double vexpo = Expo(rv), vfrac = Frac(rv);
615: release((value)ru), release((value)rv);
616: if (uexpo > Maxexpo) uexpo = Maxexpo;
617: if (vexpo > Maxexpo) vexpo = Maxexpo;
618: return (value) mk_approx(
619: atan2(
620: vexpo < Minexpo ? 0.0 : ldexp(vfrac, (int)vexpo),
621: uexpo < Minexpo ? 0.0 : ldexp(ufrac, (int)uexpo)),
622: 0.0);
623: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.