/* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */ /* \$Header: b1num.c,v 1.4 85/08/22 16:51:59 timo Exp \$ */ /* B numbers, basic external interface */ #include "b.h" #include "b0con.h" #include "b1obj.h" #include "b1num.h" /* * This file contains operations on numbers that are not predefined * B functions (the latter are defined in `bfun.c'). * This includes conversion of numeric B values to C `int's and * `double's, (numval() and intval()), * but also utilities for comparing numeric values and hashing a numeric * value to something usable as initialization for the random generator * without chance of overflow (so numval(v) is unusable). * It is also possible to build numbers of all types using mk_integer, * mk_exact (or mk_rat) and mk_approx. Note the rather irregular * (historical!) argument structure for these: mk_approx has a * `double' argument, where mk_exact and mk_rat have two values * which must be `integer' (not `int's!) as arguments. * The random number generator used by the DRAW and CHOOSE statements * is also in this file. */ /* * ival is used internally by intval() and large(). * It converts an integer to double precision, yielding * a huge value when overflow occurs (but giving no error). */ Hidden double ival(u) integer u; { double x = 0; register int i; if (IsSmallInt(u)) return SmallIntVal(u); for (i = Length(u)-1; i >= 0; --i) { if (x >= Maxreal/BASE) return Msd(u) < 0 ? -Maxreal : Maxreal; x = x*BASE + Digit(u, i); } return x; } /* Determine if a value may be converted to an int */ Visible bool large(v) value v; { double r; if (!Is_number(v) || !integral(v)) { error(MESS(1300, "number not an integer")); return No; } if (Rational(v)) v= (value) Numerator((rational)v); r= ival((integer)v); if (r > Maxint) return Yes; return No; } /* return the C `int' value of a B numeric value, if it exists. */ Visible int intval(v) value v; { /* v must be an Integral number or a Rational with Denominator==1 which may result from n round x [via mk_exact]!. */ double i; if (IsSmallInt(v)) i= SmallIntVal(v); else { if (!Is_number(v)) syserr(MESS(1301, "intval on non-number")); if (!integral(v)) { error(MESS(1302, "number not an integer")); return 0; } if (Rational(v)) v= (value) Numerator((rational)v); i= ival((integer)v); } if (i > Maxint || i < -Maxint) { error(MESS(1303, "exceedingly large integer")); return 0; } return (int) i; } /* convert an int to an intlet */ Visible intlet propintlet(i) int i; { if (i > Maxintlet || i < -Maxintlet) { error(MESS(1304, "exceedingly large integer")); return 0; } return i; } /* * determine if a number is integer */ Visible bool integral(v) value v; { if (Integral(v) || (Rational(v) && Denominator((rational)v) == int_1)) return Yes; else return No; } /* * mk_exact makes an exact number out of two integers. * The third argument is the desired number of digits after the decimal * point when the number is printed (see comments in `bfun.c' for * `round' and code in `bnuC.c'). * This printing size (if positive) is hidden in the otherwise nearly * unused * 'Size' field of the value struct * (cf. definition of Rational(v) etc.). */ Visible value mk_exact(p, q, len) integer p, q; int len; { rational r = mk_rat(p, q, len); if (Denominator(r) == int_1 && len <= 0) { integer i = (integer) Copy(Numerator(r)); release((value) r); return (value) i; } return (value) r; } #define uSMALL 1 #define uINT 2 #define uRAT 4 #define uAPP 8 #define vSMALL 16 #define vINT 32 #define vRAT 64 #define vAPP 128 Visible relation numcomp(u, v) value u, v; { int tu, tv; relation s; if (IsSmallInt(u)) tu = uSMALL; else if (Length(u) >= 0) tu = uINT; else if (Length(u) <= -2) tu = uRAT; else tu = uAPP; if (IsSmallInt(v)) tv = vSMALL; else if (Length(v) >= 0) tv = vINT; else if (Length(v) <= -2) tv = vRAT; else tv = vAPP; switch (tu|tv) { case uSMALL|vSMALL: return SmallIntVal(u) - SmallIntVal(v); case uSMALL|vINT: case uINT|vSMALL: case uINT|vINT: return int_comp((integer)u, (integer)v); case uSMALL|vRAT: case uINT|vRAT: u = (value) mk_rat((integer)u, int_1, 0); s = rat_comp((rational)u, (rational)v); release(u); return s; case uRAT|vRAT: return rat_comp((rational)u, (rational)v); case uRAT|vSMALL: case uRAT|vINT: v = (value) mk_rat((integer)v, int_1, 0); s = rat_comp((rational)u, (rational)v); release(v); return s; case uSMALL|vAPP: case uINT|vAPP: case uRAT|vAPP: u = approximate(u); s = app_comp((real)u, (real)v); release(u); return s == 0 ? -1 : s; /* u < ~u = v */ case uAPP|vAPP: return app_comp((real)u, (real)v); case uAPP|vSMALL: case uAPP|vINT: case uAPP|vRAT: v = approximate(v); s = app_comp((real)u, (real)v); release(v); return s == 0 ? 1 : s; /* u = ~v > v */ default: syserr(MESS(1305, "num_comp")); /* NOTREACHED */ } } /* * Deliver 10**n, where n is a (maybe negative!) C integer. * The result is a value (integer or rational, actually). */ Visible value tento(n) int n; { if (n < 0) { integer i= int_tento(-n); value v= (value) mk_exact(int_1, i, 0); release((value) i); return v; } return (value) int_tento(n); } /* * numval returns the numerical value of any numeric B value * as a C `double'. */ Visible double numval(u) value u; { double expo, frac; if (!Is_number(u)) { error(MESS(1306, "value not a number")); return 0.0; } u = approximate(u); expo = Expo((real)u), frac = Frac((real)u); release(u); if (expo > Maxexpo) { error(MESS(1307, "approximate number too large to be handled")); return 0.0; } if(expo < Minexpo) return 0.0; return ldexp(frac, (int)expo); } /* * Random numbers */ /* * numhash produces a `double' number that depends on the value of * v, useful for initializing the random generator. * Needs rewriting, so that for instance numhash(n) never equals n. * The magic numbers here are chosen at random. */ Visible double numhash(v) value v; { if (Integral(v)) { double d = 0; register int i; if (IsSmallInt(v)) return SmallIntVal(v); for (i = Length(v) - 1; i >= 0; --i) { d *= 2; d += Digit((integer)v, i); } return d; } if (Rational(v)) return .777 * numhash((value) Numerator((rational)v)) + .123 * numhash((value) Denominator((rational)v)); return numval(v); /* Which fails for HUGE reals. Never mind. */ } /* Initialize the random generator */ double lastran; Hidden Procedure setran (seed) double seed; { double x; x = seed >= 0 ? seed : -seed; while (x >= 1) x /= 10; lastran = floor(67108864.0*x); } Visible Procedure set_random(v) value v; { setran((double) hash(v)); } /* Return a random number in [0, 1). */ Visible value random() { double p; p = 26353589.0 * lastran + 1; lastran = p - 67108864.0*floor(p/67108864.0); return (value) mk_approx(lastran / 67108864.0, 0.0); } Visible Procedure initnum() { rat_init(); setran((double) SEED); } Visible Procedure endnum() { endrat(); }