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

   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

Defined variables