|
|
BSD 4.3
/* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
/*
$Header: /var/lib/cvsd/repos/CSRG/43BSD/contrib/B/src/bint/b1num.c,v 1.1.1.1 2018/04/24 16:12:54 root 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();
}
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.