Return to b1nuQ.c CVS log

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

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

1.1     ! root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
        !             2: 
        !             3: /*
        !             4:   $Header: b1nuQ.c,v 1.4 85/08/22 16:51:40 timo Exp $
        !             5: */
        !             6: 
        !             7: #include "b.h"
        !             8: #include "b1obj.h"
        !             9: #include "b0con.h"
        !            10: #include "b1num.h"
        !            11: 
        !            12: 
        !            13: /* Product of integer and single "digit" */
        !            14: 
        !            15: Visible integer int1mul(v, n1) integer v; digit n1; {
        !            16:        integer a;
        !            17:        digit save, bigcarry, carry = 0;
        !            18:        double z, zz, n = n1;
        !            19:        register int i;
        !            20:        struct integer vv;
        !            21: 
        !            22:        FreezeSmallInt(v, vv);
        !            23: 
        !            24:        a = (integer) grab_num(Length(v)+2);
        !            25: 
        !            26:        for (i = 0; i < Length(v); ++i) {
        !            27:                z = Digit(v,i) * n;
        !            28:                bigcarry = zz = floor(z/BASE);
        !            29:                carry += z - zz*BASE;
        !            30:                Digit(a,i) = save = Modulo(carry, BASE);
        !            31:                carry = (carry-save)/BASE + bigcarry;
        !            32:        }
        !            33: 
        !            34:        Digit(a,i) = save = Modulo(carry, BASE);
        !            35:        Digit(a,i+1) = (carry-save)/BASE;
        !            36: 
        !            37:        return int_canon(a);
        !            38: }
        !            39: 
        !            40: 
        !            41: /* Quotient of positive integer and single "digit" > 0 */
        !            42: 
        !            43: Hidden integer int1div(v, n1, prem) integer v; digit n1, *prem; {
        !            44:        integer q;
        !            45:        double r_over_n, r = 0, n = n1;
        !            46:        register int i;
        !            47:        struct integer vv;
        !            48: 
        !            49:        FreezeSmallInt(v, vv);
        !            50: 
        !            51:        q = (integer) grab_num(Length(v));
        !            52:        for (i = Length(v)-1; i >= 0; --i) {
        !            53:                r = r*BASE + Digit(v,i);
        !            54:                Digit(q,i) = r_over_n = floor(r/n);
        !            55:                r -= r_over_n * n;
        !            56:        }
        !            57:        if (prem)
        !            58:                *prem = r;
        !            59:        return int_canon(q);
        !            60: }
        !            61: 
        !            62: 
        !            63: /* Long division routine, gives access to division algorithm. */
        !            64: 
        !            65: Visible digit int_ldiv(v1, w1, pquot, prem) integer v1, w1, *pquot, *prem; {
        !            66:        integer a;
        !            67:        int sign = 1, rel_v = 0, rel_w = 0;
        !            68:        digit div, rem;
        !            69:        struct integer vv1, ww1;
        !            70: 
        !            71:        if (w1 == int_0) syserr(MESS(1100, "zero division (int_ldiv)"));
        !            72: 
        !            73:        /* Make v, w positive */
        !            74:        if (Msd(v1) < 0) {
        !            75:                sign = -1;
        !            76:                ++rel_v;
        !            77:                v1 = int_neg(v1);
        !            78:        }
        !            79: 
        !            80:        if (Msd(w1) < 0) {
        !            81:                sign *= -1;
        !            82:                ++rel_w;
        !            83:                w1 = int_neg(w1);
        !            84:        }
        !            85:        
        !            86:        FreezeSmallInt(v1, vv1);
        !            87:        FreezeSmallInt(w1, ww1);
        !            88: 
        !            89:        div = sign;
        !            90: 
        !            91:        /* Check v << w or single-digit w */
        !            92:        if (Length(v1) < Length(w1)
        !            93:                || Length(v1) == Length(w1)
        !            94:                        && Digit(v1, Length(v1)-1) < Digit(w1, Length(w1)-1)) {
        !            95:                a = int_0;
        !            96:                if (prem) {
        !            97:                        if (v1 == &vv1) *prem= (integer) MkSmallInt(Digit(v1,0));
        !            98:                        else *prem = (integer) Copy(v1);
        !            99:                }
        !           100:        }
        !           101:        else if (Length(w1) == 1) {
        !           102:                /* Single-precision division */
        !           103:                a = int1div(v1, Digit(w1,0), &rem);
        !           104:                if (prem) *prem = mk_int((double)rem);
        !           105:        }
        !           106:        else {
        !           107:                /* Multi-precision division */
        !           108:                /* Cf. Knuth II Sec. 4.3.1. Algorithm D */
        !           109:                /* Note that we count in the reverse direction (not easier!) */
        !           110: 
        !           111:                double z, zz;
        !           112:                digit carry, save, bigcarry;
        !           113:                double q, d = BASE/(Digit(w1, Length(w1)-1)+1);
        !           114:                register int i, j, k;
        !           115:                integer v, w;
        !           116:                digit vj;
        !           117: 
        !           118:                /* Normalize: make Msd(w) >= BASE/2 by multiplying
        !           119:                   both v and w by d */
        !           120: 
        !           121:                v = int1mul(v1, (digit)d);
        !           122:                        /* v is used as accumulator, must make a copy */
        !           123:                        /* v cannot be int_1 */
        !           124:                        /* (then it would be one of the cases above) */
        !           125: 
        !           126:                if (d == 1) w = (integer) Copy(w1);
        !           127:                else w = int1mul(w1, (digit)d);
        !           128: 
        !           129:                a = (integer) grab_num(Length(v1)-Length(w)+1);
        !           130: 
        !           131:                /* Division loop */
        !           132: 
        !           133:                for (j = Length(v1), k = Length(a)-1; k >= 0; --j, --k) {
        !           134:                        vj = j >= Length(v) ? 0 : Digit(v,j);
        !           135: 
        !           136:                        /* Find trial digit */
        !           137: 
        !           138:                        if (vj == Digit(w, Length(w)-1)) q = BASE-1;
        !           139:                        else q = floor( ((double)vj*BASE
        !           140:                                        + Digit(v,j-1)) / Digit(w, Length(w)-1) );
        !           141: 
        !           142:                        /* Correct trial digit */
        !           143: 
        !           144:                        while (Digit(w,Length(w)-2) * q >
        !           145:                                ((double)vj*BASE + Digit(v,j-1)
        !           146:                                        - q*Digit(w, Length(w)-1)) *BASE + Digit(v,j-2))
        !           147:                                --q;
        !           148: 
        !           149:                        /* Subtract q*w from v */
        !           150: 
        !           151:                        carry = 0;
        !           152:                        for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
        !           153:                                z = Digit(w,i) * q;
        !           154:                                bigcarry = zz = floor(z/BASE);
        !           155:                                carry += Digit(v,i+k) - z + zz*BASE;
        !           156:                                Digit(v,i+k) =
        !           157:                                        save = Modulo(carry, BASE);
        !           158:                                carry = (carry-save)/BASE - bigcarry;
        !           159:                        }
        !           160: 
        !           161:                        if (i+k < Length(v))
        !           162:                                carry += Digit(v, i+k), Digit(v, i+k) = 0;
        !           163: 
        !           164:                        /* Add back necessary? */
        !           165: 
        !           166:                                /* It is very difficult to find test cases
        !           167:                                   where add back is necessary if BASE is large.
        !           168:                                   Thanks to Arjen Lenstra, we have v=n*n-1, w=n,
        !           169:                                   where n = 8109636009903000000 (the last six
        !           170:                                   digits are not important). */
        !           171: 
        !           172:                        if (carry == 0)         /* No */
        !           173:                                Digit(a,k) = q;
        !           174:                        else {          /* Yes, add back */
        !           175:                                if (carry != -1) syserr(MESS(1101, "int_ldiv internal failure"));
        !           176:                                Digit(a,k) = q-1;
        !           177:                                carry = 0;
        !           178:                                for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
        !           179:                                        carry += Digit(v, i+k) + Digit(w,i);
        !           180:                                        Digit(v,i+k) =
        !           181:                                                save = Modulo(carry, BASE);
        !           182:                                        carry = (carry-save)/BASE;
        !           183:                                }
        !           184:                        }
        !           185:                }       /* End for(j) */
        !           186: 
        !           187:                if (prem) *prem = int_canon(v); /* Store remainder */
        !           188:                else release((value) v);
        !           189:                div = sign*d;   /* Store normalization factor */
        !           190:                release((value) w);
        !           191:                a = int_canon(a);
        !           192:        }
        !           193: 
        !           194:        if (rel_v) release((value) v1);
        !           195:        if (rel_w) release((value) w1);
        !           196: 
        !           197:        if (sign < 0) {
        !           198:                integer temp = a;
        !           199:                a = int_neg(a);
        !           200:                release((value) temp);
        !           201:        }
        !           202: 
        !           203:        if (pquot) *pquot = a;
        !           204:        else release((value) a);
        !           205:        return div;
        !           206: }
        !           207: 
        !           208: 
        !           209: Visible integer int_quot(v, w) integer v, w; {
        !           210:        integer quo;
        !           211:        VOID int_ldiv(v, w, &quo, (integer*)0);
        !           212:        return quo;
        !           213: }
        !           214: 
        !           215: Visible integer int_mod(v, w) integer v, w; {
        !           216:        integer rem;
        !           217:        digit div;
        !           218:        bool flag;
        !           219:        div = int_ldiv(v, w, (integer*)0, &rem); /* Rem. is always positive */
        !           220:        if (rem == int_0)
        !           221:                return rem; /* v mod w = 0 */
        !           222:        flag = (div < 0);
        !           223:        if (flag || Msd(w) < 0) div = -div;
        !           224:        if (div != 1) { /* Divide by div to get proper remainder back */
        !           225:                v = int1div(rem, div, (digit*)0);
        !           226:                release((value) rem);
        !           227:                rem = v;
        !           228:        }
        !           229:        if (flag) { /* Make same sign as w */
        !           230:                if (Msd(w) < 0) v = int_sum(w, rem);
        !           231:                else v = int_diff(w, rem);
        !           232:                release((value) rem);
        !           233:                rem = v;
        !           234:        }
        !           235:        return rem;
        !           236: }

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