Return to b1nuI.c CVS log

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

Annotation of 43BSDTahoe/new/B/src/bint/b1nuI.c, revision 1.1

1.1     ! root        1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
        !             2: 
        !             3: /*
        !             4:   $Header: b1nuI.c,v 1.4 85/08/22 16:51:13 timo Exp $
        !             5: */
        !             6: 
        !             7: /* Multi-precision integer arithmetic */
        !             8: 
        !             9: #include "b.h"
        !            10: #include "b1obj.h"
        !            11: #include "b1num.h"
        !            12: #include "b0con.h"
        !            13: #include "b3err.h"
        !            14: 
        !            15: /*
        !            16:  * Number representation:
        !            17:  * ======================
        !            18:  *
        !            19:  * (Think of BASE = 10 for ordinary decimal notation.)
        !            20:  * A number is a sequence of N "digits" b1, b2, ..., bN
        !            21:  * where each bi is in {0..BASE-1}, except for negative numbers,
        !            22:  * where bN = -1.
        !            23:  * The number represented by b1, ..., bN is
        !            24:  *      b1*BASE**(N-1) + b2*BASE**(N-2) + ... + bN .
        !            25:  * The base BASE is chosen so that multiplication of two positive
        !            26:  * integers up to BASE-1 can be multiplied exactly using double
        !            27:  * precision floating point arithmetic.
        !            28:  * Also it must be possible to add two long integers between
        !            29:  * -BASE and +BASE (exclusive), giving a result between -2BASE and
        !            30:  * +2BASE.
        !            31:  * BASE must be even (so we can easily decide whether the whole
        !            32:  * number is even), and positive (to avoid all kinds of other trouble).
        !            33:  * Presently, it is restricted to a power of 10 by the I/O-conversion
        !            34:  * routines (file "b1nuC.c").
        !            35:  *
        !            36:  * Canonical representation:
        !            37:  * bN is never zero (for the number zero itself, N is zero).
        !            38:  * If bN is -1, b[N-1] is never BASE-1 .
        !            39:  * All operands are assumed te be in canonical representation.
        !            40:  * Routine "int_canon" brings a number in canonical representation.
        !            41:  *
        !            42:  * Mapping to C objects:
        !            43:  * A "digit" is an integer of type "digit", probably an "int".
        !            44:  * A number is represented as a "B-integer", i.e. something
        !            45:  * of type "integer" (which is actually a pointer to some struct).
        !            46:  * The number of digits N is extracted through the macro Length(v).
        !            47:  * The i-th digit is extracted through the macro Digit(v,N-i).
        !            48:  * (So in C, we count in a backwards direction from 0 ... n-1 !)
        !            49:  * A number is created through a call to grab_num(N), which sets
        !            50:  * N zero digits (thus not in canonical form!).
        !            51:  */
        !            52: 
        !            53: 
        !            54: /*
        !            55:  * Bring an integer into canonical form.
        !            56:  * Make a SmallInt if at all possible.
        !            57:  * NB: Work done by int_canon is duplicated by mk_integer for optimization;
        !            58:  *     if the strategy here changes, look at mk_integer, too!
        !            59:  */
        !            60: 
        !            61: Visible integer int_canon(v) integer v; {
        !            62:        register int i;
        !            63: 
        !            64:        if (IsSmallInt(v)) return v;
        !            65: 
        !            66:        for (i = Length(v) - 1; i >= 0 && Digit(v,i) == 0; --i)
        !            67:                ;
        !            68: 
        !            69:        if (i < 0) {
        !            70:                release((value) v);
        !            71:                return int_0;
        !            72:        }
        !            73: 
        !            74:        if (i == 0) {
        !            75:                digit dig = Digit(v,0);
        !            76:                release((value) v);
        !            77:                return (integer) MkSmallInt(dig);
        !            78:        }
        !            79: 
        !            80:        if (i > 0 && Digit(v,i) == -1) {
        !            81:                while (i > 0 && Digit(v, i-1) == BASE-1) --i;
        !            82:                if (i == 0) {
        !            83:                        release((value) v);
        !            84:                        return (integer) MkSmallInt(-1);
        !            85:                }
        !            86:                if (i == 1) {
        !            87:                        digit dig = Digit(v,0) - BASE;
        !            88:                        release((value) v);
        !            89:                        return (integer) MkSmallInt(dig);
        !            90:                }
        !            91:                Digit(v,i) = -1;
        !            92:        }
        !            93: 
        !            94:        if (i+1 < Length(v)) return (integer) regrab_num((value) v, i+1);
        !            95: 
        !            96:        return v;
        !            97: }
        !            98: 
        !            99: 
        !           100: /* General add/subtract subroutine */
        !           101: 
        !           102: typedef double twodigit; /* Might be long on 16 bit machines */
        !           103:        /* Should be in b0con.h */
        !           104: 
        !           105: Hidden twodigit fmodulo(x, y) twodigit x, y; {
        !           106:        return x - y * (twodigit) floor((double)x / (double)y);
        !           107: }
        !           108: 
        !           109: Visible Procedure dig_gadd(to, nto, from, nfrom, ffactor)
        !           110:        digit *to, *from; intlet nto, nfrom; digit ffactor; {
        !           111:        twodigit carry= 0;
        !           112:        twodigit factor= ffactor;
        !           113:        digit save;
        !           114: 
        !           115:        nto -= nfrom;
        !           116:        if (nto < 0)
        !           117:                syserr(MESS(1000, "dig_gadd: nto < nfrom"));
        !           118:        for (; nfrom > 0; ++to, ++from, --nfrom) {
        !           119:                carry += *to + *from * factor;
        !           120:                *to= save= fmodulo(carry, (twodigit)BASE);
        !           121:                carry= (carry-save) / BASE;
        !           122:        }
        !           123:        for (; nto > 0; ++to, --nto) {
        !           124:                if (carry == 0)
        !           125:                        return;
        !           126:                carry += *to;
        !           127:                *to= save= fmodulo(carry, (twodigit)BASE);
        !           128:                carry= (carry-save) / BASE;
        !           129:        }
        !           130:        if (carry != 0)
        !           131:                to[-1] += carry*BASE; /* Assume it's -1 */
        !           132: }
        !           133: 
        !           134: 
        !           135: /* Sum or difference of two integers */
        !           136: /* Should have its own version of dig-gadd without double precision */
        !           137: 
        !           138: Visible integer int_gadd(v, w, factor) integer v, w; intlet factor; {
        !           139:        struct integer vv, ww;
        !           140:        integer s;
        !           141:        int len, lenv, i;
        !           142: 
        !           143:        FreezeSmallInt(v, vv);
        !           144:        FreezeSmallInt(w, ww);
        !           145:        lenv= len= Length(v);
        !           146:        if (Length(w) > len)
        !           147:                len= Length(w);
        !           148:        ++len;
        !           149:        s= (integer) grab_num(len);
        !           150:        for (i= 0; i < lenv; ++i)
        !           151:                Digit(s, i)= Digit(v, i);
        !           152:        for (; i < len; ++i)
        !           153:                Digit(s, i)= 0;
        !           154:        dig_gadd(&Digit(s, 0), len, &Digit(w, 0), Length(w), (digit)factor);
        !           155:        return int_canon(s);
        !           156: }
        !           157: 
        !           158: 
        !           159: /* Product of two integers */
        !           160: 
        !           161: Visible integer int_prod(v, w) integer v, w; {
        !           162:        int i;
        !           163:        integer a;
        !           164:        struct integer vv, ww;
        !           165: 
        !           166:        if (v == int_0 || w == int_0) return int_0;
        !           167:        if (v == int_1) return (integer) Copy(w);
        !           168:        if (w == int_1) return (integer) Copy(v);
        !           169: 
        !           170:        FreezeSmallInt(v, vv);
        !           171:        FreezeSmallInt(w, ww);
        !           172: 
        !           173:        a = (integer) grab_num(Length(v) + Length(w));
        !           174: 
        !           175:        for (i= Length(a)-1; i >= 0; --i)
        !           176:                Digit(a, i)= 0;
        !           177:        for (i = 0; i < Length(v) && !interrupted; ++i)
        !           178:                dig_gadd(&Digit(a, i), Length(w)+1, &Digit(w, 0), Length(w), 
        !           179:                        Digit(v, i));
        !           180: 
        !           181:        return int_canon(a);
        !           182: }
        !           183: 
        !           184: 
        !           185: /* Compare two integers */
        !           186: 
        !           187: Visible relation int_comp(v, w) integer v, w; {
        !           188:        int sv, sw;
        !           189:        register int i;
        !           190:        struct integer vv, ww;
        !           191: 
        !           192:        /* 1. Compare pointers and equal SmallInts */
        !           193:        if (v == w) return 0;
        !           194: 
        !           195:        /* 1a. Handle SmallInts */
        !           196:        if (IsSmallInt(v) && IsSmallInt(w))
        !           197:                return SmallIntVal(v) - SmallIntVal(w);
        !           198:        FreezeSmallInt(v, vv);
        !           199:        FreezeSmallInt(w, ww);
        !           200: 
        !           201:        /* 2. Extract signs */
        !           202:        sv = Length(v)==0 ? 0 : Digit(v,Length(v)-1)<0 ? -1 : 1;
        !           203:        sw = Length(w)==0 ? 0 : Digit(w,Length(w)-1)<0 ? -1 : 1;
        !           204: 
        !           205:        /* 3. Compare signs */
        !           206:        if (sv != sw) return (sv>sw) - (sv<sw);
        !           207: 
        !           208:        /* 4. Compare sizes */
        !           209:        if (Length(v) != Length(w))
        !           210:                return sv * ( (Length(v)>Length(w)) - (Length(v)<Length(w)) );
        !           211: 
        !           212:        /* 5. Compare individual digits */
        !           213:        for (i = Length(v)-1; i >= 0 && Digit(v,i) == Digit(w,i); --i)
        !           214:                ;
        !           215: 
        !           216:        /* 6. All digits equal? */
        !           217:        if (i < 0) return 0;  /* Yes */
        !           218: 
        !           219:        /* 7. Compare leftmost different digits */
        !           220:        if (Digit(v,i) < Digit(w,i)) return -1;
        !           221: 
        !           222:        return 1;
        !           223: }
        !           224: 
        !           225: 
        !           226: /* Construct an integer out of a floating point number */
        !           227: 
        !           228: #define GRAN 8 /* Granularity used when requesting more storage */
        !           229:                /* MOVE TO MEM! */
        !           230: Visible integer mk_int(x) double x; {
        !           231:        register integer a;
        !           232:        integer b;
        !           233:        register int i, j;
        !           234:        int negate;
        !           235: 
        !           236:        if (MinSmallInt <= x && x <= MaxSmallInt)
        !           237:                return (integer) MkSmallInt((int)x);
        !           238: 
        !           239:        a = (integer) grab_num(1);
        !           240:        negate = x < 0 ? 1 : 0;
        !           241:        if (negate) x = -x;
        !           242: 
        !           243:        for (i = 0; x != 0; ++i) {
        !           244:                double z = floor(x/BASE);
        !           245:                digit save = Modulo((digit)(x-z*BASE), BASE);
        !           246:                if (i >= Length(a)) {
        !           247:                        a = (integer) regrab_num((value) a, Length(a)+GRAN);
        !           248:                        for (j = Length(a)-1; j > i; --j)
        !           249:                                Digit(a,j) = 0; /* clear higher digits */
        !           250:                }
        !           251:                Digit(a,i) = save;
        !           252:                x = floor((x-save)/BASE);
        !           253:        }
        !           254: 
        !           255:        if (negate) {
        !           256:                b = int_neg(a);
        !           257:                release((value) a);
        !           258:                return b;
        !           259:        }
        !           260: 
        !           261:        return int_canon(a);
        !           262: }
        !           263: 
        !           264: /* Construct an integer out of a C int.  Like mk_int, but optimized. */
        !           265: 
        !           266: Visible value mk_integer(x) int x; {
        !           267:        if (MinSmallInt <= x && x <= MaxSmallInt) return MkSmallInt(x);
        !           268:        return (value) mk_int((double)x);
        !           269: }
        !           270: 
        !           271: 
        !           272: /* Efficiently compute 10**n as a B integer, where n is a C int >= 0 */
        !           273: 
        !           274: Visible integer int_tento(n) int n; {
        !           275:        integer i;
        !           276:        digit msd = 1;
        !           277:        if (n < 0) syserr(MESS(1001, "int_tento(-n)"));
        !           278:        if (n < tenlogBASE) {
        !           279:                while (n != 0) msd *= 10, --n;
        !           280:                return (integer) MkSmallInt(msd);
        !           281:        }
        !           282:        i = (integer) grab_num(1 + (int)(n/tenlogBASE));
        !           283:        n %= tenlogBASE;
        !           284:        while (n != 0) msd *= 10, --n;
        !           285:        Digit(i, Length(i)-1) = msd;
        !           286:        return i;
        !           287: }
        !           288: 
        !           289: #ifdef NOT_USED
        !           290: /* Approximate ceiling(10 log abs(u/v)), as C int.
        !           291:    It only works for v > 0, u, v both integers.
        !           292:    The result may be one too large or too small */
        !           293: 
        !           294: Visible int scale(u, v) integer u, v; {
        !           295:        int s;
        !           296:        double z;
        !           297:        struct integer uu, vv;
        !           298: 
        !           299:        if (Msd(v) <= 0) syserr(MESS(1002, "scale(u,v<=0)"));
        !           300:        if (u == int_0) return 0; /* `Don't care' case */
        !           301:        FreezeSmallInt(u, uu);
        !           302:        FreezeSmallInt(v, vv);
        !           303:        s = (Length(u) - Length(v)) * tenlogBASE;
        !           304:        if (Digit(u, Length(u)-1) >= 0) z = Digit(u, Length(u)-1);
        !           305:        else {
        !           306:                s -= tenlogBASE;
        !           307:                if (Length(u) == 1) z = 1;
        !           308:                else z = BASE - Digit(u, Length(u)-2);
        !           309:        }
        !           310:        z /= Digit(v, Length(v)-1);
        !           311:        while (z >= 10) z /= 10, ++s;
        !           312:        while (z < 1) z *= 10, --s;
        !           313:        return s;
        !           314: }
        !           315: #endif NOT_USED

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