2.11BSD#431 - /usr/src/new/PORT/B/src/bsmall/B1fun.c

   1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */
   2: /* $Header: B1fun.c,v 1.1 84/06/28 00:48:54 timo Exp $ */
   3: 
   4: /* B functions */
   5: #include "b.h"
   6: #include "b1obj.h"
   7: #include "b2sem.h"
   8: #include "B1num.h"
   9: 
  10: #define Maxlen(x, y) (Length(x) > Length(y) ? Length(x) : Length(y))
  11: #define Sumlen(x, y) (Length(x) + Length(y))
  12: 
  13: value sum(x, y) register value x, y; {
  14:     value r, z;
  15:     Checknum(x); Checknum(y);
  16:     if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) + Numval(y));
  17:     z= mk_exact(Denominator(x), Denominator(y), 0);
  18:     r= mk_exact(Numerator(x)*Denominator(z)+Numerator(y)*Numerator(z),
  19:             Denominator(x)*Denominator(z), Maxlen(x, y));
  20:     release(z);
  21:     return r;
  22: }
  23: 
  24: value negated(x) register value x; {
  25:     Checknum(x);
  26:     if (!Exact(x)) return mk_approx(-Numval(x));
  27:     return mk_exact(-Numerator(x), Denominator(x), Length(x));
  28: }
  29: 
  30: value diff(x, y) register value x, y; {
  31:     value r, my;
  32:     r= sum(x, my= negated(y));
  33:     release(my);
  34:     return r;
  35: }
  36: 
  37: value inverse(x) value x;{
  38:     Checknum(x);
  39:     if (Numval(x) == 0) error("in x/y, y is zero");
  40:     if (!Exact(x)) return mk_approx(1.0/Numval(x));
  41:     return mk_exact(Denominator(x), Numerator(x), Length(x));
  42: }
  43: 
  44: value prod(x, y) register value x, y; {
  45:     value a, b, r;
  46:     Checknum(x); Checknum(y);
  47:     if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) * Numval(y));
  48:     a= mk_exact(Numerator(x), Denominator(y), 0);
  49:     b= mk_exact(Numerator(y), Denominator(x), 0);
  50:     r= mk_exact(Numerator(a)*Numerator(b), Denominator(a)*Denominator(b),
  51:                         Sumlen(x, y));
  52:     release(a); release(b);
  53:     return r;
  54: }
  55: 
  56: value quot(x, y) register value x, y; {
  57:     value r, iy;
  58:     r= prod(x, iy= inverse(y));
  59:     release(iy);
  60:     return r;
  61: }
  62: 
  63: #define Even(x) ((x) == Two*floor((x)/Two))
  64: value power(x, y) register value x, y; {
  65:     Checknum(x); Checknum(y);
  66:     if (Exact(y)) {
  67:         integer py= Numerator(y), qy= Denominator(y);
  68:         if (Integral(y) && Exact(x)) {
  69:             integer px, qx, ppx, pqx, Ppx, Pqx;
  70:             if (py == Zero) return mk_int(One);
  71:             if (py > Zero) {
  72:                 px= Numerator(x);
  73:                 qx= Denominator(x);
  74:             } else {
  75:                 py= -py;
  76:                 px= Denominator(x);
  77:                 qx= Numerator(x);
  78:             }
  79:             ppx= pqx= One;
  80:             Ppx= px; Pqx= qx;
  81:             while (py >= Two) {
  82:                 if (!Even(py)) {
  83:                     ppx*= Ppx; pqx*= Pqx;
  84:                 }
  85:                 Ppx*= Ppx; Pqx*= Pqx;
  86:                 py= floor(py/Two);
  87:             }
  88:             ppx*= Ppx; pqx*= Pqx;
  89:             return mk_exact(ppx, pqx, 0);
  90:         } /* END Integral(y) && Exact(x) */
  91:         else {
  92:             double vx= Numval(x);
  93:             short sx= vx < 0 ? -1 : vx == 0 ? 0 : 1;
  94:             if (sx < 0 && Even(qy))
  95:                 error("in x**(p/q), x is negative and q is even");
  96:             if (sx == 0 && py < Zero)
  97:                 error("0**y with negative y");
  98:             if (sx < 0 && Even(py)) sx= 1;
  99:             return mk_approx(sx * pow(fabs(vx), py/qy));
 100:         }
 101:     } /* END Exact(y) */
 102:     else {
 103:         double vx= Numval(x), vy= Approxval(y);
 104:         if (vy == 0) return mk_approx(1.0);
 105:         if (vx < 0)
 106:             error("in x**y, x is negative and y is not exact");
 107:         if (vx == 0 && vy < 0)
 108:             error("0E0**y with negative y");
 109:         return mk_approx(pow(vx, vy));
 110:     }
 111: }
 112: 
 113: value root2(n, x) register value n, x; {
 114:     value r, in;
 115:     Checknum(n);
 116:     if (Numval(n) == 0) error("in x root y, x is zero");
 117:     r= power(x, in= inverse(n));
 118:     release(in);
 119:     return r;
 120: }
 121: 
 122: value absval(x) register value x; {
 123:     Checknum(x);
 124:     if (!Exact(x)) return mk_approx(fabs(Numval(x)));
 125:     return mk_exact((integer) fabs((double) Numerator(x)), Denominator(x), Length(x));
 126: }
 127: 
 128: value signum(x) register value x; {
 129:     double v= numval(x);
 130:     return mk_int(v < 0 ? -One : v == 0 ? Zero : One);
 131: }
 132: 
 133: value floorf(x) register value x; {
 134:     return mk_int(floor(numval(x)));
 135: }
 136: 
 137: value ceilf(x) register value x; {
 138:     return mk_int(ceil(numval(x)));
 139: }
 140: 
 141: value round1(x) register value x; {
 142:     return mk_int(floor(numval(x) + .5));
 143: }
 144: 
 145: value round2(n, x) register value n, x; {
 146:     value ten, tenp, xtenp, r0, r;
 147:     Checknum(n);
 148:     if (!Integral(n)) error("in n round x, n is not an integer");
 149:     ten= mk_integer(10);
 150:     tenp= power(ten, n);
 151:     xtenp= prod(x, tenp);
 152:     r0= round1(xtenp);
 153:     r= mk_exact(Numerator(r0), Numerator(tenp), propintlet((int) Numerator(n)));
 154:     release(ten); release(tenp); release(xtenp); release(r0);
 155:     return r;
 156: }
 157: 
 158: value mod(a, n) register value a, n; {
 159:     value f, p, d;
 160:     Checknum(a); Checknum(n);
 161:     f= mk_int(floor(Numval(a) / Numval(n)));
 162:     p= prod(n, f);
 163:     d= diff(a, p);
 164:     release(f); release(p);
 165:     return d;
 166: }
 167: 
 168: double lastran;
 169: 
 170: setran (seed) double seed;
 171: {double x;
 172:  x= seed >= 0 ? seed : -seed;
 173:  while (x >= 1) x/= 10;
 174:  lastran= floor(67108864.0 * x);
 175: }
 176: 
 177: set_random(v) value v; {
 178:     setran((double) hash(v));
 179: }
 180: 
 181: value random() /* 0 <= r < 1 */
 182: {double p;
 183:  p= 26353589.0 * lastran + 1;
 184:  lastran= p - 67108864.0 * floor (p / 67108864.0);
 185:  return mk_approx(lastran / 67108864.0);
 186: }
 187: 
 188: value root1(v) value v; {
 189:     value two= mk_integer(2);
 190:     v= root2(two, v);
 191:     release(two);
 192:     return(v);
 193: }
 194: 
 195: value pi() { return mk_approx(3.141592653589793238462); }
 196: value e() { return mk_approx(exp(1.0)); }
 197: 
 198: value sin1(v) value v; { return mk_approx(sin(numval(v))); }
 199: value cos1(v) value v; { return mk_approx(cos(numval(v))); }
 200: value tan1(v) value v; { return mk_approx(tan(numval(v))); }
 201: value atn1(v) value v; { return mk_approx(atan(numval(v))); }
 202: value exp1(v) value v; { return mk_approx(exp(numval(v))); }
 203: value log1(v) value v; { return mk_approx(log(numval(v))); }
 204: 
 205: value log2(u, v) value u, v;{
 206:     return mk_approx(log(numval(v)) / log(numval(u)));
 207: }
 208: 
 209: value atn2(u, v) value u, v; {
 210:     return mk_approx(atan2(numval(v), numval(u)));
 211: }

Defined functions

Defined variables

Defined macros