Annotation of 43BSDTahoe/new/B/src/bint/b1num.c, revision 1.1.1.1

1.1       root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
                      2: 
                      3: /*
                      4:   $Header: b1num.c,v 1.4 85/08/22 16:51:59 timo Exp $
                      5: */
                      6: 
                      7: /* B numbers, basic external interface */
                      8: 
                      9: #include "b.h"
                     10: #include "b0con.h"
                     11: #include "b1obj.h"
                     12: #include "b1num.h"
                     13: 
                     14: /*
                     15:  * This file contains operations on numbers that are not predefined
                     16:  * B functions (the latter are defined in `bfun.c').
                     17:  * This includes conversion of numeric B values to C `int's and
                     18:  * `double's, (numval() and intval()),
                     19:  * but also utilities for comparing numeric values and hashing a numeric
                     20:  * value to something usable as initialization for the random generator
                     21:  * without chance of overflow (so numval(v) is unusable).
                     22:  * It is also possible to build numbers of all types using mk_integer,
                     23:  * mk_exact (or mk_rat) and mk_approx.  Note the rather irregular
                     24:  * (historical!) argument structure for these: mk_approx has a
                     25:  * `double' argument, where mk_exact and mk_rat have two values
                     26:  * which must be `integer' (not `int's!) as arguments.
                     27:  * The random number generator used by the DRAW and CHOOSE statements
                     28:  * is also in this file.
                     29:  */
                     30: 
                     31: /*
                     32:  * ival is used internally by intval() and large().
                     33:  * It converts an integer to double precision, yielding
                     34:  * a huge value when overflow occurs (but giving no error).
                     35:  */
                     36: 
                     37: Hidden double ival(u) integer u; {
                     38:        double x = 0;
                     39:        register int i;
                     40: 
                     41:        if (IsSmallInt(u)) return SmallIntVal(u);
                     42:        for (i = Length(u)-1; i >= 0; --i) {
                     43:                if (x >= Maxreal/BASE)
                     44:                        return Msd(u) < 0 ? -Maxreal : Maxreal;
                     45:                x = x*BASE + Digit(u, i);
                     46:        }
                     47: 
                     48:        return x;
                     49: }
                     50: 
                     51: 
                     52: /* Determine if a value may be converted to an int */
                     53: 
                     54: Visible bool large(v) value v; {
                     55:        double r;
                     56:        if (!Is_number(v) || !integral(v)) {
                     57:                error(MESS(1300, "number not an integer"));
                     58:                return No;
                     59:        }
                     60:        if (Rational(v)) v= (value) Numerator((rational)v);
                     61:        r= ival((integer)v);
                     62:        if (r > Maxint) return Yes;
                     63:        return No;
                     64: }
                     65: 
                     66: /* return the C `int' value of a B numeric value, if it exists. */
                     67: 
                     68: Visible int intval(v) value v; {
                     69:        /* v must be an Integral number or a Rational with Denominator==1
                     70:            which may result from n round x [via mk_exact]!. */
                     71:        double i;
                     72:        if (IsSmallInt(v)) i= SmallIntVal(v);
                     73:        else {
                     74:                if (!Is_number(v)) syserr(MESS(1301, "intval on non-number"));
                     75:                if (!integral(v)) {
                     76:                        error(MESS(1302, "number not an integer"));
                     77:                        return 0;
                     78:                }
                     79:                if (Rational(v)) v= (value) Numerator((rational)v);
                     80:                i= ival((integer)v);
                     81:        }
                     82:        if (i > Maxint || i < -Maxint) {
                     83:                error(MESS(1303, "exceedingly large integer"));
                     84:                return 0;
                     85:        }
                     86:        return (int) i;
                     87: }
                     88: 
                     89: 
                     90: /* convert an int to an intlet */
                     91: 
                     92: Visible intlet propintlet(i) int i; {
                     93:        if (i > Maxintlet || i < -Maxintlet) {
                     94:                error(MESS(1304, "exceedingly large integer"));
                     95:                return 0;
                     96:        }
                     97:        return i;
                     98: }
                     99: 
                    100: 
                    101: /*
                    102:  * determine if a number is integer
                    103:  */
                    104: 
                    105: Visible bool integral(v) value v; {
                    106:        if (Integral(v) || (Rational(v) && Denominator((rational)v) == int_1))
                    107:                return Yes;
                    108:        else return No;
                    109: }
                    110: 
                    111: /*
                    112:  * mk_exact makes an exact number out of two integers.
                    113:  * The third argument is the desired number of digits after the decimal
                    114:  * point when the number is printed (see comments in `bfun.c' for
                    115:  * `round' and code in `bnuC.c').
                    116:  * This printing size (if positive) is hidden in the otherwise nearly
                    117:  * unused * 'Size' field of the value struct
                    118:  * (cf. definition of Rational(v) etc.).
                    119:  */
                    120: 
                    121: Visible value mk_exact(p, q, len) integer p, q; int len; {
                    122:        rational r = mk_rat(p, q, len);
                    123: 
                    124:        if (Denominator(r) == int_1 && len <= 0) {
                    125:                integer i = (integer) Copy(Numerator(r));
                    126:                release((value) r);
                    127:                return (value) i;
                    128:        }
                    129: 
                    130:        return (value) r;
                    131: }
                    132: 
                    133: #define uSMALL 1
                    134: #define uINT 2
                    135: #define uRAT 4
                    136: #define uAPP 8
                    137: #define vSMALL 16
                    138: #define vINT 32
                    139: #define vRAT 64
                    140: #define vAPP 128
                    141: 
                    142: Visible relation numcomp(u, v) value u, v; {
                    143:        int tu, tv; relation s;
                    144: 
                    145:        if (IsSmallInt(u)) tu = uSMALL;
                    146:        else if (Length(u) >= 0) tu = uINT;
                    147:        else if (Length(u) <= -2) tu = uRAT;
                    148:        else tu = uAPP;
                    149:        if (IsSmallInt(v)) tv = vSMALL;
                    150:        else if (Length(v) >= 0) tv = vINT;
                    151:        else if (Length(v) <= -2) tv = vRAT;
                    152:        else tv = vAPP;
                    153: 
                    154:        switch (tu|tv) {
                    155: 
                    156:        case uSMALL|vSMALL: return SmallIntVal(u) - SmallIntVal(v);
                    157: 
                    158:        case uSMALL|vINT:
                    159:        case uINT|vSMALL:
                    160:        case uINT|vINT: return int_comp((integer)u, (integer)v);
                    161: 
                    162:        case uSMALL|vRAT:
                    163:        case uINT|vRAT:
                    164:                u = (value) mk_rat((integer)u, int_1, 0);
                    165:                s = rat_comp((rational)u, (rational)v);
                    166:                release(u);
                    167:                return s;
                    168: 
                    169:        case uRAT|vRAT:
                    170:                return rat_comp((rational)u, (rational)v);
                    171: 
                    172:        case uRAT|vSMALL:
                    173:        case uRAT|vINT:
                    174:                v = (value) mk_rat((integer)v, int_1, 0);
                    175:                s = rat_comp((rational)u, (rational)v);
                    176:                release(v);
                    177:                return s;
                    178: 
                    179:        case uSMALL|vAPP:
                    180:        case uINT|vAPP:
                    181:        case uRAT|vAPP:
                    182:                u = approximate(u);
                    183:                s = app_comp((real)u, (real)v);
                    184:                release(u);
                    185:                return s == 0 ? -1 : s; /* u < ~u = v */
                    186: 
                    187:        case uAPP|vAPP:
                    188:                return app_comp((real)u, (real)v);
                    189: 
                    190:        case uAPP|vSMALL:
                    191:        case uAPP|vINT:
                    192:        case uAPP|vRAT:
                    193:                v = approximate(v);
                    194:                s = app_comp((real)u, (real)v);
                    195:                release(v);
                    196:                return s == 0 ? 1 : s; /* u = ~v > v */
                    197: 
                    198:        default: syserr(MESS(1305, "num_comp")); /* NOTREACHED */
                    199: 
                    200:        }
                    201: }
                    202: 
                    203: 
                    204: /*
                    205:  * Deliver 10**n, where n is a (maybe negative!) C integer.
                    206:  * The result is a value (integer or rational, actually).
                    207:  */
                    208: 
                    209: Visible value tento(n) int n; {
                    210:        if (n < 0) {
                    211:                integer i= int_tento(-n);
                    212:                value v= (value) mk_exact(int_1, i, 0);
                    213:                release((value) i);
                    214:                return v;
                    215:        }
                    216:        return (value) int_tento(n);
                    217: }
                    218: 
                    219: 
                    220: /*
                    221:  * numval returns the numerical value of any numeric B value
                    222:  * as a C `double'.
                    223:  */
                    224: 
                    225: Visible double numval(u) value u; {
                    226:        double expo, frac;
                    227: 
                    228:        if (!Is_number(u)) {
                    229:                error(MESS(1306, "value not a number"));
                    230:                return 0.0;
                    231:        }
                    232:        u = approximate(u);
                    233:        expo = Expo((real)u), frac = Frac((real)u);
                    234:        release(u);
                    235:        if (expo > Maxexpo) {
                    236:                error(MESS(1307, "approximate number too large to be handled"));
                    237:                return 0.0;
                    238:        }
                    239:        if(expo < Minexpo) return 0.0;
                    240:        return ldexp(frac, (int)expo);
                    241: }
                    242: 
                    243: 
                    244: /*
                    245:  * Random numbers
                    246:  */
                    247: 
                    248: 
                    249: /*
                    250:  * numhash produces a `double' number that depends on the value of
                    251:  * v, useful for initializing the random generator.
                    252:  * Needs rewriting, so that for instance numhash(n) never equals n.
                    253:  * The magic numbers here are chosen at random.
                    254:  */
                    255: 
                    256: Visible double numhash(v) value v; {
                    257:        if (Integral(v)) {
                    258:                double d = 0;
                    259:                register int i;
                    260: 
                    261:                if (IsSmallInt(v)) return SmallIntVal(v);
                    262:                for (i = Length(v) - 1; i >= 0; --i) {
                    263:                        d *= 2;
                    264:                        d += Digit((integer)v, i);
                    265:                }
                    266: 
                    267:                return d;
                    268:        }
                    269: 
                    270:        if (Rational(v))
                    271:                return .777 * numhash((value) Numerator((rational)v)) +
                    272:                       .123 * numhash((value) Denominator((rational)v));
                    273: 
                    274:        return numval(v); /* Which fails for HUGE reals.  Never mind. */
                    275: }
                    276: 
                    277: 
                    278: /* Initialize the random generator */
                    279: 
                    280: double lastran;
                    281: 
                    282: Hidden Procedure setran (seed) double seed; {
                    283:        double x;
                    284: 
                    285:        x = seed >= 0 ? seed : -seed;
                    286:        while (x >= 1) x /= 10;
                    287:        lastran = floor(67108864.0*x);
                    288: }
                    289: 
                    290: Visible Procedure set_random(v) value v; {
                    291:        setran((double) hash(v));
                    292: }
                    293: 
                    294: 
                    295: /* Return a random number in [0, 1). */
                    296: 
                    297: Visible value random() {
                    298:        double p;
                    299: 
                    300:        p = 26353589.0 * lastran + 1;
                    301:        lastran = p - 67108864.0*floor(p/67108864.0);
                    302: 
                    303:        return (value) mk_approx(lastran / 67108864.0, 0.0);
                    304: }
                    305: 
                    306: Visible Procedure initnum() {
                    307:        rat_init();
                    308:        setran((double) SEED);
                    309: }
                    310: 
                    311: Visible Procedure endnum() {
                    312:        endrat();
                    313: }

unix.superglobalmegacorp.com

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.