```   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: }
```

#### Defined functions

atn1 defined in line 579; used 1 times
atn2 defined in line 611; used 1 times
cos1 defined in line 576; used 1 times
denominator defined in line 496; used 1 times
dyop defined in line 41; used 4 times
e defined in line 544; used 6 times
exp1 defined in line 588; used 1 times
idiff defined in line 97; used 1 times
isum defined in line 87; used 1 times
• in line 94
log1 defined in line 595; used 3 times
mod defined in line 168; used 1 times
negated defined in line 139; used 3 times
numerator defined in line 479; used 1 times
pi defined in line 543; used 1 times
quot defined in line 130; used 8 times
rational defined in line 24; used 21 times
root1 defined in line 531; used 1 times
root2 defined in line 513; used 2 times
round1 defined in line 322; used 2 times
round2 defined in line 352; used 1 times
signum defined in line 392; used 1 times
sin1 defined in line 575; used 1 times
tan1 defined in line 577; used 1 times
trig defined in line 546; used 3 times
xxx_quot defined in line 120; used 1 times

#### Defined variables

Hidden defined in line 546; never used
Visible defined in line 544; never used
 Last modified: 1985-08-27 Generated: 2016-12-26 Generated by src2html V0.67 page hit count: 3235