Return to b1fun.c CVS log

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

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

1.1     ! root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
        !             2: 
        !             3: /*
        !             4:   $Header: b1fun.c,v 1.4 85/08/22 16:48:24 timo Exp $
        !             5: */
        !             6: 
        !             7: /* Functions defined on numeric values. */
        !             8: 
        !             9: #include <errno.h> /* For EDOM and ERANGE */
        !            10: 
        !            11: #include "b.h"
        !            12: #include "b0con.h"
        !            13: #include "b1obj.h"
        !            14: #include "b1num.h"
        !            15: 
        !            16: /*
        !            17:  * The visible routines here implement predefined B arithmetic operators,
        !            18:  * taking one or two numeric values as operands, and returning a numeric
        !            19:  * value.
        !            20:  * No type checking of operands is done: this must be done by the caller.
        !            21:  */
        !            22: 
        !            23: typedef value (*valfun)();
        !            24: typedef rational (*ratfun)();
        !            25: typedef real (*appfun)();
        !            26: typedef double (*mathfun)();
        !            27: 
        !            28: /*
        !            29:  * For the arithmetic functions (+, -, *, /) the same action is needed:
        !            30:  * 1) if both operands are Integral, use function from int_* submodule;
        !            31:  * 2) if both are Exact, use function from rat_* submodule (after possibly
        !            32:  *    converting one of them from Integral to Rational);
        !            33:  * 3) otherwise, make both approximate and use function from app_*
        !            34:  *    submodule.
        !            35:  * The functions performing the appropriate action for each of the submodules
        !            36:  * are passed as parameters.
        !            37:  * Division is a slight exception, since i/j can be a rational.
        !            38:  * See `quot' below.
        !            39:  */
        !            40: 
        !            41: Hidden value dyop(u, v, int_fun, rat_fun, app_fun)
        !            42:        value u, v;
        !            43:        valfun int_fun;
        !            44:        ratfun rat_fun;
        !            45:        appfun app_fun;
        !            46: {
        !            47:        if (Integral(u) && Integral(v)) /* Use integral operation */
        !            48:                return (*int_fun)(u, v);
        !            49: 
        !            50:        if (Exact(u) && Exact(v)) {
        !            51:                rational u1, v1, a;
        !            52: 
        !            53:                /* Use rational operation */
        !            54: 
        !            55:                u1 = Integral(u) ? mk_rat((integer)u, int_1, 0) :
        !            56:                                (rational) Copy(u);
        !            57:                v1 = Integral(v) ? mk_rat((integer)v, int_1, 0) :
        !            58:                                (rational) Copy(v);
        !            59:                a = (*rat_fun)(u1, v1);
        !            60:                release((value) u1);
        !            61:                release((value) v1);
        !            62: 
        !            63:                if (Denominator(a) == int_1 && Roundsize(a) == 0) {
        !            64:                        integer b = (integer) Copy(Numerator(a));
        !            65:                        release((value) a);
        !            66:                        return (value) b;
        !            67:                }
        !            68: 
        !            69:                return (value) a;
        !            70:        }
        !            71: 
        !            72:        /* Use approximate operation */
        !            73: 
        !            74:        {
        !            75:                real u1, v1, a;
        !            76:                u1 = Approximate(u) ? (real) Copy(u) : (real) approximate(u);
        !            77:                v1 = Approximate(v) ? (real) Copy(v) : (real) approximate(v);
        !            78:                a = (*app_fun)(u1, v1);
        !            79:                release((value) u1);
        !            80:                release((value) v1);
        !            81: 
        !            82:                return (value) a;
        !            83:        }
        !            84: }
        !            85: 
        !            86: 
        !            87: Hidden integer isum(u, v) integer u, v; {
        !            88:        return int_sum(u, v);
        !            89: }
        !            90: 
        !            91: Visible value sum(u, v) value u, v; {
        !            92:        if (IsSmallInt(u) && IsSmallInt(v))
        !            93:                return mk_integer(SmallIntVal(u) + SmallIntVal(v));
        !            94:        return dyop(u, v, (value (*)())isum, rat_sum, app_sum);
        !            95: }
        !            96: 
        !            97: Hidden integer idiff(u, v) integer u, v; {
        !            98:        return int_diff(u, v);
        !            99: }
        !           100: 
        !           101: Visible value diff(u, v) value u, v; {
        !           102:        if (IsSmallInt(u) && IsSmallInt(v))
        !           103:                return mk_integer(SmallIntVal(u) - SmallIntVal(v));
        !           104:        return dyop(u, v, (value (*)())idiff, rat_diff, app_diff);
        !           105: }
        !           106: 
        !           107: Visible value prod(u, v) value u, v; {
        !           108:        if (IsSmallInt(u) && IsSmallInt(v))
        !           109:                return (value)
        !           110:                        mk_int((double)SmallIntVal(u) * (double)SmallIntVal(v));
        !           111:        return dyop(u, v, (value (*)())int_prod, rat_prod, app_prod);
        !           112: }
        !           113: 
        !           114: 
        !           115: /*
        !           116:  * We cannot use int_quot (which performs integer division with truncation).
        !           117:  * Here is the routine we need.
        !           118:  */
        !           119: 
        !           120: Hidden value xxx_quot(u, v) integer u, v; {
        !           121: 
        !           122:        if (v == int_0) {
        !           123:                error(MESS(400, "in i/j, j is zero"));
        !           124:                return (value) Copy(u);
        !           125:        }
        !           126: 
        !           127:        return mk_exact(u, v, 0);
        !           128: }
        !           129: 
        !           130: Visible value quot(u, v) value u, v; {
        !           131:        return dyop(u, v, xxx_quot, rat_quot, app_quot);
        !           132: }
        !           133: 
        !           134: 
        !           135: /*
        !           136:  * Unary minus and abs follow the same principle but with only one operand.
        !           137:  */
        !           138: 
        !           139: Visible value negated(u) value u; {
        !           140:        if (IsSmallInt(u)) return mk_integer(-SmallIntVal(u));
        !           141:        if (Integral(u))
        !           142:                return (value) int_neg((integer)u);
        !           143:        if (Rational(u))
        !           144:                return (value) rat_neg((rational)u);
        !           145:        return (value) app_neg((real)u);
        !           146: }
        !           147: 
        !           148: 
        !           149: Visible value absval(u) value u; {
        !           150:        if (Integral(u)) {
        !           151:                if (Msd((integer)u) < 0)
        !           152:                        return (value) int_neg((integer)u);
        !           153:        } else if (Rational(u)) {
        !           154:                if (Msd(Numerator((rational)u)) < 0)
        !           155:                        return (value) rat_neg((rational)u);
        !           156:        } else if (Approximate(u) && Frac((real)u) < 0)
        !           157:                return (value) app_neg((real)u);
        !           158: 
        !           159:        return Copy(u);
        !           160: }
        !           161: 
        !           162: 
        !           163: /*
        !           164:  * The remaining operators follow less similar paths and some of
        !           165:  * them contain quite subtle code.
        !           166:  */
        !           167: 
        !           168: Visible value mod(u, v) value u, v; {
        !           169:        value p, q, m, n;
        !           170: 
        !           171:        if (v == (value)int_0 ||
        !           172:                Rational(v) && Numerator((rational)v) == int_0 ||
        !           173:                Approximate(v) && Frac((real)v) == 0) {
        !           174:                error(MESS(401, "in x mod y, y is zero"));
        !           175:                return Copy(u);
        !           176:        }
        !           177: 
        !           178:        if (Integral(u) && Integral(v))
        !           179:                return (value) int_mod((integer)u, (integer)v);
        !           180: 
        !           181:        /* Compute `u-v*floor(u/v)', as in the formal definition of `u mod v'. */
        !           182: 
        !           183:        q = quot(u, v);
        !           184:        n = floorf(q);
        !           185:        release(q);
        !           186:        p = prod(n, v);
        !           187:        release(n);
        !           188:        m = diff(u, p);
        !           189:        release(p);
        !           190: 
        !           191:        return m;
        !           192: }
        !           193: 
        !           194: 
        !           195: /*
        !           196:  * u**v has the most special cases of all the predefined arithmetic functions.
        !           197:  */
        !           198: 
        !           199: Visible value power(u, v) value u, v; {
        !           200:        real ru, rv, rw;
        !           201:        if (Exact(u) && (Integral(v) ||
        !           202:                        /* Next check catches for integers disguised as rationals: */
        !           203:                        Rational(v) && Denominator((rational)v) == int_1)) {
        !           204:                rational a;
        !           205:                integer b = Integral(v) ? (integer)v : Numerator((rational)v);
        !           206:                        /* Now b is really an integer. */
        !           207: 
        !           208:                u = Integral(u) ? (value) mk_rat((integer)u, int_1, 0) :
        !           209:                                Copy(u);
        !           210:                        
        !           211:                a = rat_power((rational)u, b);
        !           212:                release(u);
        !           213:                if (Denominator(a) == int_1) { /* Make integral result */
        !           214:                        b = (integer) Copy(Numerator(a));
        !           215:                        release((value) a);
        !           216:                        return (value)b;
        !           217:                }
        !           218:                return (value)a;
        !           219:        }
        !           220: 
        !           221:        if (Exact(v)) {
        !           222:                integer vn, vd;
        !           223:                int s;
        !           224:                ru = (real) approximate(u);
        !           225:                s = (Frac(ru) > 0) - (Frac(ru) < 0);
        !           226: 
        !           227:                if (s < 0) rv = app_neg(ru), release((value)ru), ru = rv;
        !           228:                if (Integral(v)) {
        !           229:                        vn = (integer)v;
        !           230:                        vd = int_1;
        !           231:                } else {
        !           232:                        vd = Denominator((rational)v);
        !           233:                        if (s < 0 && Even(Lsd(vd)))
        !           234:                                error(MESS(402, "in x**(p/q), x is negative and q is even"));
        !           235:                        vn = Numerator((rational)v);
        !           236:                }
        !           237:                if (vn == int_0) {
        !           238:                        release((value)ru);
        !           239:                        return one;
        !           240:                }
        !           241:                if (s == 0 && Msd(vn) < 0) {
        !           242:                        error(MESS(403, "in 0**y, y is negative"));
        !           243:                        return (value) ru;
        !           244:                }
        !           245:                if (s < 0 && Even(Lsd(vn)))
        !           246:                        s = 1;
        !           247:                rv = (real) approximate(v);
        !           248:                rw = app_power(ru, rv);
        !           249:                release((value)ru), release((value)rv);
        !           250:                if (s < 0) ru = app_neg(rw), release((value)rw), rw = ru;
        !           251:                return (value) rw;
        !           252:        }
        !           253: 
        !           254:        /* Everything else: we now know u or v is approximate */
        !           255: 
        !           256:        ru = (real) approximate(u);
        !           257:        if (Frac(ru) < 0) {
        !           258:                error(MESS(404, "in x**y, x is negative and y is not exact"));
        !           259:                return (value) ru;
        !           260:        }
        !           261:        rv = (real) approximate(v);
        !           262:        if (Frac(ru) == 0 && Frac(rv) < 0) {
        !           263:                error(MESS(405, "in 0**y, y is negative"));
        !           264:                release((value)rv);
        !           265:                return (value) ru;
        !           266:        }
        !           267:        rw = app_power(ru, rv);
        !           268:        release((value)ru), release((value)rv);
        !           269:        return (value) rw;
        !           270: }
        !           271: 
        !           272: 
        !           273: /*
        !           274:  * floor: for approximate numbers app_floor() is used;
        !           275:  * for integers it is a no-op; other exact numbers effectively calculate
        !           276:  * u - (u mod 1).
        !           277:  */
        !           278: 
        !           279: Visible value floorf(u) value u; {
        !           280:        integer quo, rem, v;
        !           281:        digit div;
        !           282: 
        !           283:        if (Integral(u)) return Copy(u);
        !           284:        if (Approximate(u)) return app_floor((real)u);
        !           285: 
        !           286:        /* It is a rational number */
        !           287: 
        !           288:        div = int_ldiv(Numerator((rational)u), Denominator((rational)u),
        !           289:                &quo, &rem);
        !           290:        if (div < 0 && rem != int_0) { /* Correction for negative noninteger */
        !           291:                v = int_diff(quo, int_1);
        !           292:                release((value) quo);
        !           293:                quo = v;
        !           294:        }
        !           295:        release((value) rem);
        !           296:        return (value) quo;
        !           297: }
        !           298: 
        !           299: 
        !           300: /*
        !           301:  * ceiling x is defined as -floor(-x);
        !           302:  * and that's how it's implemented, except for integers.
        !           303:  */
        !           304: 
        !           305: Visible value ceilf(u) value u; {
        !           306:        value v;
        !           307:        if (Integral(u)) return Copy(u);
        !           308:        u = negated(u);
        !           309:        v = floorf(u);
        !           310:        release(u);
        !           311:        u = negated(v);
        !           312:        release(v);
        !           313:        return u;
        !           314: }
        !           315: 
        !           316: 
        !           317: /*
        !           318:  * round u is defined as floor(u+0.5), which is what is done here,
        !           319:  * except for integers which are left unchanged.
        !           320:  */
        !           321: 
        !           322: Visible value round1(u) value u; {
        !           323:        value v, w;
        !           324: 
        !           325:        if (Integral(u)) return Copy(u);
        !           326: 
        !           327:        v = sum(u, (value)rat_half);
        !           328:        w = floorf(v);
        !           329:        release(v);
        !           330: 
        !           331:        return w;
        !           332: }
        !           333: 
        !           334: 
        !           335: /*
        !           336:  * u round v is defined as 10**-u * round(v*10**u).
        !           337:  * A complication is that u round v is always printed with exactly u digits
        !           338:  * after the decimal point, even if this involves trailing zeros,
        !           339:  * or if v is an integer.
        !           340:  * Consequently, the result is always kept as a rational, even if it can be
        !           341:  * simplified to an integer, and the size field of the rational number
        !           342:  * (which is made negative to distinguish it from integers, and < -1 to
        !           343:  * distinguish it from approximate numbers) is used to store the number of
        !           344:  * significant digits.
        !           345:  * Thus a size of -2 means a normal rational number, and a size < -2
        !           346:  * means a rounded number to be printed with (-2 - length) digits
        !           347:  * after the decimal point.  This last expression can be retrieved using
        !           348:  * the macro Roundsize(v) which should only be applied to Rational
        !           349:  * numbers.
        !           350:  */
        !           351: 
        !           352: Visible value round2(u, v) value u, v; {
        !           353:        value w, tenp;
        !           354:        int i;
        !           355: 
        !           356:        if (!Integral(u)) {
        !           357:                error(MESS(406, "in n round x, n is not an integer"));
        !           358:                i = 0;
        !           359:        } else
        !           360:                i = propintlet(intval(u));
        !           361: 
        !           362:        tenp = tento(i);
        !           363: 
        !           364:        v = prod(v, tenp);
        !           365:        w = round1(v);
        !           366: 
        !           367:        release(v);
        !           368:        v = quot(w, tenp);
        !           369:        release(w);
        !           370:        release(tenp);
        !           371: 
        !           372:        if (i > 0) {    /* Set number of digits to be printed */
        !           373:                if (propintlet(-2 - i) < -2) {
        !           374:                        if (Rational(v))
        !           375:                                Length(v) = -2 - i;
        !           376:                        else if (Integral(v)) {
        !           377:                                w = v;
        !           378:                                v = mk_exact((integer) w, int_1, i);
        !           379:                                release(w);
        !           380:                        }
        !           381:                }
        !           382:        }
        !           383: 
        !           384:        return v;
        !           385: }
        !           386: 
        !           387: 
        !           388: /*
        !           389:  * sign u inspects the sign of either u, u's numerator or u's fractional part.
        !           390:  */
        !           391: 
        !           392: Visible value signum(u) value u; {
        !           393:        int s;
        !           394: 
        !           395:        if (Exact(u)) {
        !           396:                if (Rational(u))
        !           397:                        u = (value) Numerator((rational)u);
        !           398:                s = u==(value)int_0 ? 0 : Msd((integer)u) < 0 ? -1 : 1;
        !           399:        } else
        !           400:                s = Frac((real)u) > 0 ? 1 : Frac((real)u) < 0 ? -1 : 0;
        !           401: 
        !           402:        return MkSmallInt(s);
        !           403: }
        !           404: 
        !           405: 
        !           406: /*
        !           407:  * ~u makes an approximate number of any numerical value.
        !           408:  */
        !           409: 
        !           410: Visible value approximate(u) value u; {
        !           411:        double x, e, expo;
        !           412:        register int i;
        !           413:        if (Approximate(u)) return Copy(u);
        !           414:        if (IsSmallInt(u))
        !           415:                x = SmallIntVal(u), expo = 0;
        !           416:        else if (Rational(u)) {
        !           417:                rational r = (rational) u;
        !           418:                integer p = Numerator(r), q;
        !           419:                double xp, xq;
        !           420:                int lp, lq;
        !           421:                if (p == int_0)
        !           422:                        return Copy((value)app_0);
        !           423:                q = Denominator(r);
        !           424:                lp = IsSmallInt(p) ? 1 : Length(p);
        !           425:                lq = IsSmallInt(q) ? 1 : Length(q);
        !           426:                xp = 0, xq = 0;
        !           427:                expo = lp - lq;
        !           428:                if (IsSmallInt(p)) { xp = SmallIntVal(p); --expo; }
        !           429:                else {
        !           430:                        for (i = Length(p)-1; i >= 0; --i) {
        !           431:                                xp = xp*BASE + Digit(p, i);
        !           432:                                --expo;
        !           433:                                if (xp < -Maxreal/BASE || xp > Maxreal/BASE)
        !           434:                                        break;
        !           435:                        }
        !           436:                }
        !           437:                if (IsSmallInt(q)) { xq = SmallIntVal(q); ++expo; }
        !           438:                else {
        !           439:                        for (i = Length(q)-1; i >= 0; --i) {
        !           440:                                xq = xq*BASE + Digit(q, i);
        !           441:                                ++expo;
        !           442:                                if (xq > Maxreal/BASE)
        !           443:                                        break;
        !           444:                        }
        !           445:                }
        !           446:                x = xp/xq;
        !           447:        }
        !           448:        else if (Integral(u)) {
        !           449:                integer p = (integer) u;
        !           450:                if (Length(p) == 0)
        !           451:                        return Copy((value)app_0);
        !           452:                x = 0;
        !           453:                expo = Length(p);
        !           454:                for (i = Length(p)-1; i >= 0; --i) {
        !           455:                        x = x*BASE + Digit(p, i);
        !           456:                        --expo;
        !           457:                        if (x < -Maxreal/BASE || x > Maxreal/BASE)
        !           458:                                break;
        !           459:                }
        !           460:        }
        !           461:        while (expo < 0 && fabs(x) > BASE/Maxreal) ++expo, x /= BASE;
        !           462:        while (expo > 0 && fabs(x) < Maxreal/BASE) --expo, x *= BASE;
        !           463:        if (expo != 0) {
        !           464:                expo = expo*twologBASE + 1;
        !           465:                e = floor(expo);
        !           466:                x *= .5 * exp((expo-e) * logtwo);
        !           467:        }
        !           468:        else
        !           469:                e = 0;
        !           470:        return (value) mk_approx(x, e);
        !           471: }
        !           472: 
        !           473: 
        !           474: /*
        !           475:  * numerator v returns the numerator of v, whenever v is an exact number.
        !           476:  * For integers, that is v itself.
        !           477:  */
        !           478: 
        !           479: Visible value numerator(v) value v; {
        !           480:        if (!Exact(v)) {
        !           481:                error(MESS(407, "*/ on approximate number"));
        !           482:                return zero;
        !           483:        }
        !           484: 
        !           485:        if (Integral(v)) return Copy(v);
        !           486: 
        !           487:        return (value) Copy(Numerator((rational)v));
        !           488: }
        !           489: 
        !           490: 
        !           491: /*
        !           492:  * /*v returns the denominator of v, whenever v is an exact number.
        !           493:  * For integers, that is 1.
        !           494:  */
        !           495: 
        !           496: Visible value denominator(v) value v; {
        !           497:        if (!Exact(v)) {
        !           498:                error(MESS(408, "/* on approximate number"));
        !           499:                return zero;
        !           500:        }
        !           501: 
        !           502:        if (Integral(v)) return one;
        !           503: 
        !           504:        return (value) Copy(Denominator((rational)v));
        !           505: }
        !           506: 
        !           507: 
        !           508: /*
        !           509:  * u root v is defined as v**(1/u), where u is usually but need not be
        !           510:  * an integer.
        !           511:  */
        !           512: 
        !           513: Visible value root2(u, v) value u, v; {
        !           514:        if (u == (value)int_0 ||
        !           515:                Rational(u) && Numerator((rational)u) == int_0 ||
        !           516:                Approximate(u) && Frac((real)u) == 0) {
        !           517:                error(MESS(409, "in n root x, n is zero"));
        !           518:                v = Copy(v);
        !           519:        } else {
        !           520:                u = quot((value)int_1, u);
        !           521:                v = power(v, u);
        !           522:                release(u);
        !           523:        }
        !           524: 
        !           525:        return v;
        !           526: }
        !           527: 
        !           528: /* root x is computed more exactly than n root x, by doing
        !           529:    one iteration step extra.  This ~guarantees root(n**2) = n. */
        !           530: 
        !           531: Visible value root1(v) value v; {
        !           532:        value r = root2((value)int_2, v);
        !           533:        value v_over_r, theirsum, result;
        !           534:        if (Approximate(r) && Frac((real)r) == 0.0) return (value)r;
        !           535:        v_over_r = quot(v, r);
        !           536:        theirsum = sum(r, v_over_r), release(r), release(v_over_r);
        !           537:        result = quot(theirsum, (value)int_2), release(theirsum);
        !           538:        return result;
        !           539: }
        !           540: 
        !           541: /* The rest of the mathematical functions */
        !           542: 
        !           543: Visible value pi() { return (value) mk_approx(3.141592653589793238463, 0.0); }
        !           544: Visible value e() { return (value) mk_approx(2.718281828459045235360, 0.0); }
        !           545: 
        !           546: Hidden value trig(v, fun, zeroflag) value v; mathfun fun; bool zeroflag; {
        !           547:        real w = (real) approximate(v);
        !           548:        double expo = Expo(w), frac = Frac(w), x, result;
        !           549:        extern int errno;
        !           550:        if (expo <= Minexpo/2) {
        !           551:                if (zeroflag) return (value) w; /* sin small x = x, etc. */
        !           552:                frac = 0, expo = 0;
        !           553:        }
        !           554:        release((value)w);
        !           555:        if (expo > Maxexpo) errno = EDOM;
        !           556:        else {
        !           557:                x = ldexp(frac, (int)expo);
        !           558:                if (x >= Maxtrig || x <= -Maxtrig) errno = EDOM;
        !           559:                else {
        !           560:                        errno = 0;
        !           561:                        result = (*fun)(x);
        !           562:                }
        !           563:        }
        !           564:        if (errno != 0) {
        !           565:                if (errno == ERANGE)
        !           566:                        error(MESS(410, "the result is too large to be representable"));
        !           567:                else if (errno == EDOM)
        !           568:                        error(MESS(411, "x is too large to compute a meaningful result"));
        !           569:                else error(MESS(412, "math library error"));
        !           570:                return Copy((value)app_0);
        !           571:        }
        !           572:        return (value) mk_approx(result, 0.0);
        !           573: }
        !           574: 
        !           575: Visible value sin1(v) value v; { return trig(v, sin, Yes); }
        !           576: Visible value cos1(v) value v; { return trig(v, cos, No); }
        !           577: Visible value tan1(v) value v; { return trig(v, tan, Yes); }
        !           578: 
        !           579: Visible value atn1(v) value v; {
        !           580:        real w = (real) approximate(v);
        !           581:        double expo = Expo(w), frac = Frac(w);
        !           582:        if (expo <= Minexpo + 2) return (value) w; /* atan of small x = x */
        !           583:        release((value)w);
        !           584:        if (expo > Maxexpo) expo = Maxexpo;
        !           585:        return (value) mk_approx(atan(ldexp(frac, (int)expo)), 0.0);
        !           586: }
        !           587: 
        !           588: Visible value exp1(v) value v; {
        !           589:        real w = (real) approximate(v);
        !           590:        real x = app_exp(w);
        !           591:        release((value)w);
        !           592:        return (value) x;
        !           593: }
        !           594: 
        !           595: Visible value log1(v) value v; {
        !           596:        real w = (real) approximate(v);
        !           597:        real x = app_log(w);
        !           598:        release((value)w);
        !           599:        return (value) x;
        !           600: }
        !           601: 
        !           602: Visible value log2(u, v) value u, v;{
        !           603:        value w;
        !           604:        u = log1(u);
        !           605:        v = log1(v);
        !           606:        w = quot(v, u);
        !           607:        release(u), release(v);
        !           608:        return w;
        !           609: }
        !           610: 
        !           611: Visible value atn2(u, v) value u, v; {
        !           612:        real ru = (real) approximate(u), rv = (real) approximate(v);
        !           613:        double uexpo = Expo(ru), ufrac = Frac(ru);
        !           614:        double vexpo = Expo(rv), vfrac = Frac(rv);
        !           615:        release((value)ru), release((value)rv);
        !           616:        if (uexpo > Maxexpo) uexpo = Maxexpo;
        !           617:        if (vexpo > Maxexpo) vexpo = Maxexpo;
        !           618:        return (value) mk_approx(
        !           619:                atan2(
        !           620:                        vexpo < Minexpo ? 0.0 : ldexp(vfrac, (int)vexpo),
        !           621:                        uexpo < Minexpo ? 0.0 : ldexp(ufrac, (int)uexpo)),
        !           622:                0.0);
        !           623: }

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