/* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */ /* \$Header: B1fun.c,v 1.1 84/06/28 00:48:54 timo Exp \$ */ /* B functions */ #include "b.h" #include "b1obj.h" #include "b2sem.h" #include "B1num.h" #define Maxlen(x, y) (Length(x) > Length(y) ? Length(x) : Length(y)) #define Sumlen(x, y) (Length(x) + Length(y)) value sum(x, y) register value x, y; { value r, z; Checknum(x); Checknum(y); if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) + Numval(y)); z= mk_exact(Denominator(x), Denominator(y), 0); r= mk_exact(Numerator(x)*Denominator(z)+Numerator(y)*Numerator(z), Denominator(x)*Denominator(z), Maxlen(x, y)); release(z); return r; } value negated(x) register value x; { Checknum(x); if (!Exact(x)) return mk_approx(-Numval(x)); return mk_exact(-Numerator(x), Denominator(x), Length(x)); } value diff(x, y) register value x, y; { value r, my; r= sum(x, my= negated(y)); release(my); return r; } value inverse(x) value x;{ Checknum(x); if (Numval(x) == 0) error("in x/y, y is zero"); if (!Exact(x)) return mk_approx(1.0/Numval(x)); return mk_exact(Denominator(x), Numerator(x), Length(x)); } value prod(x, y) register value x, y; { value a, b, r; Checknum(x); Checknum(y); if (!Exact(x) || !Exact(y)) return mk_approx(Numval(x) * Numval(y)); a= mk_exact(Numerator(x), Denominator(y), 0); b= mk_exact(Numerator(y), Denominator(x), 0); r= mk_exact(Numerator(a)*Numerator(b), Denominator(a)*Denominator(b), Sumlen(x, y)); release(a); release(b); return r; } value quot(x, y) register value x, y; { value r, iy; r= prod(x, iy= inverse(y)); release(iy); return r; } #define Even(x) ((x) == Two*floor((x)/Two)) value power(x, y) register value x, y; { Checknum(x); Checknum(y); if (Exact(y)) { integer py= Numerator(y), qy= Denominator(y); if (Integral(y) && Exact(x)) { integer px, qx, ppx, pqx, Ppx, Pqx; if (py == Zero) return mk_int(One); if (py > Zero) { px= Numerator(x); qx= Denominator(x); } else { py= -py; px= Denominator(x); qx= Numerator(x); } ppx= pqx= One; Ppx= px; Pqx= qx; while (py >= Two) { if (!Even(py)) { ppx*= Ppx; pqx*= Pqx; } Ppx*= Ppx; Pqx*= Pqx; py= floor(py/Two); } ppx*= Ppx; pqx*= Pqx; return mk_exact(ppx, pqx, 0); } /* END Integral(y) && Exact(x) */ else { double vx= Numval(x); short sx= vx < 0 ? -1 : vx == 0 ? 0 : 1; if (sx < 0 && Even(qy)) error("in x**(p/q), x is negative and q is even"); if (sx == 0 && py < Zero) error("0**y with negative y"); if (sx < 0 && Even(py)) sx= 1; return mk_approx(sx * pow(fabs(vx), py/qy)); } } /* END Exact(y) */ else { double vx= Numval(x), vy= Approxval(y); if (vy == 0) return mk_approx(1.0); if (vx < 0) error("in x**y, x is negative and y is not exact"); if (vx == 0 && vy < 0) error("0E0**y with negative y"); return mk_approx(pow(vx, vy)); } } value root2(n, x) register value n, x; { value r, in; Checknum(n); if (Numval(n) == 0) error("in x root y, x is zero"); r= power(x, in= inverse(n)); release(in); return r; } value absval(x) register value x; { Checknum(x); if (!Exact(x)) return mk_approx(fabs(Numval(x))); return mk_exact((integer) fabs((double) Numerator(x)), Denominator(x), Length(x)); } value signum(x) register value x; { double v= numval(x); return mk_int(v < 0 ? -One : v == 0 ? Zero : One); } value floorf(x) register value x; { return mk_int(floor(numval(x))); } value ceilf(x) register value x; { return mk_int(ceil(numval(x))); } value round1(x) register value x; { return mk_int(floor(numval(x) + .5)); } value round2(n, x) register value n, x; { value ten, tenp, xtenp, r0, r; Checknum(n); if (!Integral(n)) error("in n round x, n is not an integer"); ten= mk_integer(10); tenp= power(ten, n); xtenp= prod(x, tenp); r0= round1(xtenp); r= mk_exact(Numerator(r0), Numerator(tenp), propintlet((int) Numerator(n))); release(ten); release(tenp); release(xtenp); release(r0); return r; } value mod(a, n) register value a, n; { value f, p, d; Checknum(a); Checknum(n); f= mk_int(floor(Numval(a) / Numval(n))); p= prod(n, f); d= diff(a, p); release(f); release(p); return d; } double lastran; setran (seed) double seed; {double x; x= seed >= 0 ? seed : -seed; while (x >= 1) x/= 10; lastran= floor(67108864.0 * x); } set_random(v) value v; { setran((double) hash(v)); } value random() /* 0 <= r < 1 */ {double p; p= 26353589.0 * lastran + 1; lastran= p - 67108864.0 * floor (p / 67108864.0); return mk_approx(lastran / 67108864.0); } value root1(v) value v; { value two= mk_integer(2); v= root2(two, v); release(two); return(v); } value pi() { return mk_approx(3.141592653589793238462); } value e() { return mk_approx(exp(1.0)); } value sin1(v) value v; { return mk_approx(sin(numval(v))); } value cos1(v) value v; { return mk_approx(cos(numval(v))); } value tan1(v) value v; { return mk_approx(tan(numval(v))); } value atn1(v) value v; { return mk_approx(atan(numval(v))); } value exp1(v) value v; { return mk_approx(exp(numval(v))); } value log1(v) value v; { return mk_approx(log(numval(v))); } value log2(u, v) value u, v;{ return mk_approx(log(numval(v)) / log(numval(u))); } value atn2(u, v) value u, v; { return mk_approx(atan2(numval(v), numval(u))); }