Return to b1nuA.c CVS log

Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint

Annotation of 43BSD/contrib/B/src/bint/b1nuA.c, revision 1.1

1.1     ! root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
        !             2: 
        !             3: /*
        !             4:   $Header: b1nuA.c,v 1.4 85/08/22 16:50:25 timo Exp $
        !             5: */
        !             6: 
        !             7: /* Approximate arithmetic */
        !             8: 
        !             9: #include "b.h"
        !            10: #include "b0con.h"
        !            11: #include "b1obj.h"
        !            12: #include "b1num.h"
        !            13: #include "b3err.h" /* For still_ok */
        !            14: 
        !            15: /*
        !            16: For various reasons, on some machines (notably the VAX), the range
        !            17: of the exponent is too small (ca. 1.7E38), and we cope with this by
        !            18: adding a second word which holds the exponent.
        !            19: However, on other machines (notably the IBM PC), the range is sufficient
        !            20: (ca. 1E300), and here we try to save as much code as possible by not
        !            21: doing our own exponent handling.  (To be fair, we also don't check
        !            22: certain error conditions, to save more code.)
        !            23: The difference is made by #defining EXT_RANGE (in b1num.h), meaning we
        !            24: have to EXTend the RANGE of the exponent.
        !            25: */
        !            26: 
        !            27: #ifdef EXT_RANGE
        !            28: Hidden struct real app_0_buf = {Num, 1, -1, 0.0, -(double)Maxint};
        !            29:        /* Exponent must be less than any realistic exponent! */
        !            30: #else !EXT_RANGE
        !            31: Hidden struct real app_0_buf = {Num, 1, -1, 0.0};
        !            32: #endif !EXT_RANGE
        !            33: 
        !            34: Visible real app_0 = &app_0_buf;
        !            35: 
        !            36: /*
        !            37:  * Build an approximate number.
        !            38:  */
        !            39: 
        !            40: Visible real mk_approx(frac, expo) double frac, expo; {
        !            41:        real u;
        !            42: #ifdef EXT_RANGE
        !            43:        expint e;
        !            44:        if (frac != 0) frac = frexp(frac, &e), expo += e;
        !            45:        if (frac == 0.5) frac = 1, --expo; /* Assert 0.5 < frac <= 1 */
        !            46:        if (frac == 0 || expo < -BIG) return (real) Copy(app_0);
        !            47:        if (expo > BIG) {
        !            48:                error(MESS(700, "approximate number too large"));
        !            49:                expo = BIG;
        !            50:        }
        !            51: #else !EXT_RANGE
        !            52:        if (frac == 0.0) return (real) Copy(app_0);
        !            53:        frac= ldexp(frac, (int)expo);
        !            54: #endif EXT_RANGE
        !            55:        u = (real) grab_num(-1);
        !            56:        Frac(u) = frac;
        !            57: #ifdef EXT_RANGE
        !            58:        Expo(u) = expo;
        !            59: #endif EXT_RANGE
        !            60:        return u;
        !            61: }
        !            62: 
        !            63: /*
        !            64:  * Approximate arithmetic.
        !            65:  */
        !            66: 
        !            67: Visible real app_sum(u, v) real u, v; {
        !            68: #ifdef EXT_RANGE
        !            69:        real w;
        !            70:        if (Expo(u) < Expo(v)) w = u, u = v, v = w;
        !            71:        if (Expo(v) - Expo(u) < Minexpo) return (real) Copy(u);
        !            72:        return mk_approx(Frac(u) + ldexp(Frac(v), (int)(Expo(v) - Expo(u))),
        !            73:                Expo(u));
        !            74: #else !EXT_RANGE
        !            75:        return mk_approx(Frac(u) + Frac(v), 0.0);
        !            76: #endif !EXT_RANGE
        !            77: }
        !            78: 
        !            79: Visible real app_diff(u, v) real u, v; {
        !            80: #ifdef EXT_RANGE
        !            81:        real w;
        !            82:        int sign = 1;
        !            83:        if (Expo(u) < Expo(v)) w = u, u = v, v = w, sign = -1;
        !            84:        if (Expo(v) - Expo(u) < Minexpo)
        !            85:                return sign < 0 ? app_neg(u) : (real) Copy(u);
        !            86:        return mk_approx(
        !            87:                sign * (Frac(u) - ldexp(Frac(v), (int)(Expo(v) - Expo(u)))),
        !            88:                Expo(u));
        !            89: #else !EXT_RANGE
        !            90:        return mk_approx(Frac(u) - Frac(v), 0.0);
        !            91: #endif !EXT_RANGE
        !            92: }
        !            93: 
        !            94: Visible real app_neg(u) real u; {
        !            95:        return mk_approx(-Frac(u), Expo(u));
        !            96: }
        !            97: 
        !            98: Visible real app_prod(u, v) real u, v; {
        !            99:        return mk_approx(Frac(u) * Frac(v), Expo(u) + Expo(v));
        !           100: }
        !           101: 
        !           102: Visible real app_quot(u, v) real u, v; {
        !           103:        if (Frac(v) == 0.0) {
        !           104:                error(MESS(701, "in u/v, v is zero"));
        !           105:                return (real) Copy(u);
        !           106:        }
        !           107:        return mk_approx(Frac(u) / Frac(v), Expo(u) - Expo(v));
        !           108: }
        !           109: 
        !           110: /*
        !           111:        YIELD log"(frac, expo):
        !           112:                CHECK frac > 0
        !           113:                RETURN normalize"(expo*logtwo + log(frac), 0)
        !           114: */
        !           115: 
        !           116: Visible real app_log(v) real v; {
        !           117:        double frac = Frac(v), expo = Expo(v);
        !           118:        if (frac <= 0) {
        !           119:                error(MESS(702, "in log x, x is <= 0"));
        !           120:                return (real) Copy(app_0);
        !           121:        }
        !           122:        return mk_approx(expo*logtwo + log(frac), 0.0);
        !           123: }
        !           124: 
        !           125: /*
        !           126:        YIELD exp"(frac, expo):
        !           127:                IF expo < minexpo: RETURN zero"
        !           128:                WHILE expo < 0: PUT frac/2, expo+1 IN frac, expo
        !           129:                PUT exp frac IN f
        !           130:                PUT normalize"(f, 0) IN f, e
        !           131:                WHILE expo > 0:
        !           132:                        PUT (f, e) prod" (f, e) IN f, e
        !           133:                        PUT expo-1 IN expo
        !           134:                RETURN f, e
        !           135: */
        !           136: 
        !           137: Visible real app_exp(v) real v; {
        !           138: #ifdef EXT_RANGE
        !           139:        expint ei;
        !           140:        double frac = Frac(v), expo = Expo(v), new_expo;
        !           141:        static double canexp;
        !           142:        if (!canexp) canexp = floor(log(log(Maxreal)) / logtwo);
        !           143:        if (expo <= canexp) {
        !           144:                if (expo < Minexpo) return mk_approx(1.0, 0.0);
        !           145:                frac = ldexp(frac, (int)expo);
        !           146:                expo = 0;
        !           147:        }
        !           148:        else if (expo >= Maxexpo) {
        !           149:                /* Definitely too big (the real boundary is much smaller
        !           150:                   but here we are in danger of overflowing new_expo
        !           151:                   in the loop below) */
        !           152:                return mk_approx(1.0, Maxreal); /* Force an error! */
        !           153:        }
        !           154:        else {
        !           155:                frac = ldexp(frac, (int)canexp);
        !           156:                expo -= canexp;
        !           157:        }
        !           158:        frac = exp(frac);
        !           159:        new_expo = 0;
        !           160:        while (expo > 0 && frac != 0) {
        !           161:                frac = frexp(frac, &ei);
        !           162:                new_expo += ei;
        !           163:                frac *= frac;
        !           164:                new_expo += new_expo;
        !           165:                --expo;
        !           166:        }
        !           167:        return mk_approx(frac, new_expo);
        !           168: #else !EXT_RANGE
        !           169:        return mk_approx(exp(Frac(v)), 0.0);
        !           170: #endif !EXT_RANGE
        !           171: }
        !           172: 
        !           173: /*
        !           174:        YIELD (frac, expo) power" v":
        !           175:                \   (f*2**e)**v =
        !           176:                \ = f**v * 2**(e*v) =
        !           177:                \ = f**v * 2**((e*v) mod 1) * 2**floor(e*v) .
        !           178:                PUT exp"(v" prod" normalize"(log frac, 0)) IN temp1" \ = f**v
        !           179:                PUT expo*numval(v") IN ev \ = e*v
        !           180:                PUT exp(logtwo * (ev - floor ev))) IN temp2 \ = 2**(ev mod 1)
        !           181:                PUT temp1" IN f, e
        !           182:                RETURN normalize"(f*temp2, e + floor ev)
        !           183: */
        !           184: 
        !           185: Visible real app_power(u, v) real u, v; {
        !           186:        double frac = Frac(u);
        !           187: #ifdef EXT_RANGE
        !           188:        real logfrac, vlogfrac, result;
        !           189:        double expo = Expo(u), rest;
        !           190: #endif !EXT_RANGE
        !           191:        if (frac <= 0) {
        !           192:                if (frac < 0) error(MESS(703, "in 0**v, v is negative"));
        !           193:                if (v == app_0) return mk_approx(1.0, 0.0); /* 0**0 = 1 */
        !           194:                return (real) Copy(app_0); /* 0**x = 0 */
        !           195:        }
        !           196: #ifdef EXT_RANGE
        !           197:        frac *= 2, expo -= 1; /* Renormalize to 1 < frac <= 2, so log frac > 0 */
        !           198:        logfrac = mk_approx(log(frac), 0.0);
        !           199:        vlogfrac = app_prod(v, logfrac);
        !           200:        result = app_exp(vlogfrac);
        !           201:        /* But what if result overflows but expo is very negative??? */
        !           202:        if (still_ok) {
        !           203:                expo *= numval((value)v);
        !           204:                rest = expo - floor(expo);
        !           205:                frac = Frac(result) * exp(logtwo*rest);
        !           206:                expo = Expo(result) + floor(expo);
        !           207:        }
        !           208:        release((value)logfrac), release((value)vlogfrac), release((value)result);
        !           209:        return mk_approx(frac, expo);
        !           210: #else !EXT_RANGE
        !           211:        return mk_approx(exp(log(frac) * Frac(v)), 0.0);
        !           212: #endif !EXT_RANGE
        !           213: }
        !           214: 
        !           215: Visible int app_comp(u, v) real u, v; {
        !           216:        double xu, xv;
        !           217: #ifdef EXT_RANGE
        !           218:        double eu, ev;
        !           219: #endif EXT_RANGE
        !           220:        if (u == v) return 0;
        !           221:        xu = Frac(u), xv = Frac(v);
        !           222: #ifdef EXT_RANGE
        !           223:        if (xu*xv > 0) {
        !           224:                eu = Expo(u), ev = Expo(v);
        !           225:                if (eu < ev) return xu < 0 ? 1 : -1;
        !           226:                if (eu > ev) return xu < 0 ? -1 : 1;
        !           227:        }
        !           228: #endif EXT_RANGE
        !           229:        if (xu < xv) return -1;
        !           230:        if (xu > xv) return 1;
        !           231:        return 0;
        !           232: }
        !           233: 
        !           234: Visible value app_floor(u) real u; {
        !           235: #ifdef EXT_RANGE
        !           236:        integer v, w;
        !           237:        value twotow, result;
        !           238:        if (Expo(u) <= Dblbits)
        !           239:                return (value)
        !           240:                        mk_int(floor(ldexp(Frac(u),
        !           241:                                (int)(Expo(u) < 0 ? -1 : Expo(u)))));
        !           242:        v = mk_int(ldexp(Frac(u), Dblbits));
        !           243:        w = mk_int(Expo(u) - Dblbits);
        !           244:        twotow = power((value)int_2, (value)w);
        !           245:        result = prod((value)v, twotow);
        !           246:        release((value) v), release((value) w), release(twotow);
        !           247:        if (!Integral(result)) syserr(MESS(704, "app_floor: result not integral"));
        !           248:        return result;
        !           249: #else !EXT_RANGE
        !           250:        return (value) mk_int(floor(Frac(u)));
        !           251: #endif !EXT_RANGE
        !           252: }

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