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