Return to b1nuG.c CVS log

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

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

1.1     ! root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
        !             2: 
        !             3: /*
        !             4:   $Header: b1nuG.c,v 1.4 85/08/22 16:50:59 timo Exp $
        !             5: */
        !             6: 
        !             7: #include "b.h"
        !             8: #include "b0con.h"
        !             9: #include "b1obj.h"
        !            10: #include "b1num.h"
        !            11: 
        !            12: 
        !            13: /*
        !            14:  * Routines for greatest common divisor calculation
        !            15:  * Cf. Knuth II Sec. 4.5.2. Algorithm B
        !            16:  * "Binary gcd algorithm"
        !            17:  * The labels correspond with those in the book
        !            18:  *
        !            19:  * Assumptions about built-in arithmetic:
        !            20:  * x>>1 == x/2  (if x >= 0)
        !            21:  * 1<<k == 2**k (if it fits in a word)
        !            22:  */
        !            23: 
        !            24: /* Single-precision gcd for integers > 0 */
        !            25: 
        !            26: Hidden digit dig_gcd(u, v) register digit u, v; {
        !            27:        register digit t;
        !            28:        register int k = 0;
        !            29: 
        !            30:        if (u <= 0 || v <= 0) syserr(MESS(900, "dig_gcd of number(s) <= 0"));
        !            31: 
        !            32:  /*B1*/        while (Even(u) && Even(v)) ++k, u >>= 1, v >>= 1;
        !            33: 
        !            34:  /*B2*/        if (Even(u)) t = u;
        !            35:        else {  t = -v;
        !            36:                goto B4;
        !            37:        }
        !            38: 
        !            39:        do {
        !            40:  /*B3*/                do {
        !            41:                        t /= 2;
        !            42:  B4:                   ;
        !            43:                } while (Even(t));
        !            44: 
        !            45:  /*B5*/                if (t > 0) u = t;
        !            46:                else v = -t;
        !            47: 
        !            48:  /*B6*/                t = u-v;
        !            49:        } while (t);
        !            50: 
        !            51:        return u * (1<<k);
        !            52: }
        !            53: 
        !            54: 
        !            55: Visible integer int_half(v) integer v; {
        !            56:        register int i;
        !            57:        register long carry;
        !            58: 
        !            59:        if (IsSmallInt(v))      return (integer) MkSmallInt(SmallIntVal(v) / 2);
        !            60: 
        !            61:        if (Msd(v) < 0) {
        !            62:                i = Length(v)-2;
        !            63:                if (i < 0) {
        !            64:                        release((value) v);
        !            65:                        return int_0;
        !            66:                }
        !            67:                carry = BASE;
        !            68:        }
        !            69:        else {
        !            70:                carry = 0;
        !            71:                i = Length(v)-1;
        !            72:        }
        !            73: 
        !            74:        if (Refcnt(v) > 1) uniql((value *) &v);
        !            75: 
        !            76:        for (; i >= 0; --i) {
        !            77:                carry += Digit(v,i);
        !            78:                Digit(v,i) = carry/2;
        !            79:                carry = carry&1 ? BASE : 0;
        !            80:        }
        !            81: 
        !            82:        return int_canon(v);
        !            83: }
        !            84: 
        !            85: 
        !            86: Hidden integer twice(v) integer v; {
        !            87:        digit carry = 0;
        !            88:        int i;
        !            89: 
        !            90:        if (IsSmallInt(v)) return mk_int(2.0 * SmallIntVal(v));
        !            91: 
        !            92:        if (Refcnt(v) > 1) uniql((value *) &v);
        !            93: 
        !            94:        for (i = 0; i < Length(v); ++i) {
        !            95:                carry += Digit(v,i) * 2;
        !            96:                if (carry >= BASE)
        !            97:                        Digit(v,i) = carry-BASE, carry = 1;
        !            98:                else
        !            99:                        Digit(v,i) = carry, carry = 0;
        !           100:        }
        !           101: 
        !           102:        if (carry) {    /* Enlarge the number */
        !           103:                v = (integer) regrab_num((value) v, Length(v)+1);
        !           104:                Digit(v, Length(v)-1) = carry;
        !           105:        }
        !           106: 
        !           107:        return v;
        !           108: }
        !           109: 
        !           110: 
        !           111: Hidden bool even(u) integer u; {
        !           112:        if (IsSmallInt(u))
        !           113:                return Even(SmallIntVal(u));
        !           114:        return Even(Lsd(u));
        !           115: }
        !           116: 
        !           117: 
        !           118: /* Multi-precision gcd of integers > 0 */
        !           119: 
        !           120: Visible integer int_gcd(u1, v1) integer u1, v1; {
        !           121:        struct integer uu1, vv1;
        !           122: 
        !           123:        if (Msd(u1) <= 0 || Msd(v1) <= 0)
        !           124:                syserr(MESS(901, "gcd of number(s) <= 0"));
        !           125: 
        !           126:        if (u1==int_1 || v1==int_1) return int_1;
        !           127:                /* Speed-up for e.g. 1E-100000 */
        !           128: 
        !           129:        if (IsSmallInt(u1) && IsSmallInt(v1))
        !           130:                return (integer)
        !           131:                        MkSmallInt(dig_gcd(SmallIntVal(u1), SmallIntVal(v1)));
        !           132: 
        !           133:        FreezeSmallInt(u1, uu1);
        !           134:        FreezeSmallInt(v1, vv1);
        !           135: 
        !           136:        /* Multi-precision binary gcd algorithm */
        !           137: 
        !           138:        {       long k = 0;
        !           139:                integer t, u, v;
        !           140: 
        !           141:                u = (integer) Copy((value) u1);
        !           142:                v = (integer) Copy((value) v1);
        !           143: 
        !           144:  /*B1*/                while (even(u) && even(v)) {
        !           145:                        u = int_half(u);
        !           146:                        v = int_half(v);
        !           147:                        if (++k < 0) {
        !           148:                                /*It's a number we can't cope with,
        !           149:                                  with too many common factors 2.
        !           150:                                  Though the user can't help it,
        !           151:                                  the least we can do is to allow
        !           152:                                  continuation of the session.
        !           153:                                */
        !           154:                                error(MESS(902, "exceptionally large rational number"));
        !           155:                                k = 0;
        !           156:                        }
        !           157:                }
        !           158: 
        !           159:  /*B2*/                if (even(u)) t = (integer) Copy(u);
        !           160:                else {
        !           161:                        t = int_neg(v);
        !           162:                        goto B4;
        !           163:                }
        !           164: 
        !           165:                do {
        !           166:  /*B3*/                        do {
        !           167:                                t = int_half(t);
        !           168:  B4:                           ;
        !           169:                        } while (even(t));
        !           170: 
        !           171:  /*B5*/                        if (Msd(t) >= 0) {
        !           172:                                release((value) u);
        !           173:                                u = t;
        !           174:                        }
        !           175:                        else {
        !           176:                                release((value) v);
        !           177:                                v = int_neg(t);
        !           178:                                release((value) t);
        !           179:                        }
        !           180: 
        !           181:  /*B6*/                        t = int_diff(u, v);
        !           182:                        /* t cannot be int_1 since both u and v are odd! */
        !           183:                } while (t != int_0);
        !           184: 
        !           185:                release((value) t);
        !           186:                release((value) v);
        !           187: 
        !           188:                while (--k >= 0) u = twice(u);
        !           189:                
        !           190:                return u;
        !           191:        }
        !           192: }

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