Return to mdiv.c CVS log

Up to [Research Unix] / researchv10no / cmd / lfactor

Annotation of researchv10no/cmd/lfactor/mdiv.c, revision 1.1

1.1     ! root        1: #include "mp.h"
        !             2: 
        !             3: mint
        !             4: mdiv(a,b,r)
        !             5: mint a, b, *r;
        !             6: {
        !             7:        mint c;
        !             8:        mint num, den;
        !             9:        mint quot, rem;
        !            10:        mint temp, trial;
        !            11:        mint zero, one;
        !            12:        double tquot;
        !            13:        int nsign, dsign;
        !            14: 
        !            15:        zero = itom(0);
        !            16:        one = itom(1);
        !            17:        nsign = 1;
        !            18:        dsign = 1;
        !            19:        if((b.low==0) && (b.high==0)){
        !            20:                printf("mdiv: Zero divide.\n");
        !            21:                abort();
        !            22:        }
        !            23:        num = a;
        !            24:        den = b;
        !            25:        if((num.low<0) || (num.high<0)){
        !            26:                num.high = -num.high;
        !            27:                num.low = -num.low;
        !            28:                nsign = -1;
        !            29:        }
        !            30:        if((den.low<0) || (den.high<0)){
        !            31:                den.high = -den.high;
        !            32:                den.low = -den.low;
        !            33:                dsign = -1;
        !            34:        }
        !            35:        if(mcmp(num,den) < 0){
        !            36:                *r = a;
        !            37:                return(zero);
        !            38:        }
        !            39:        tquot = (num.high*e16+num.low) / (den.high*e16+den.low);
        !            40:        modf(tquot , &tquot);
        !            41:        if(tquot < e16){
        !            42:                c.high = 0.;
        !            43:                c.low = tquot;
        !            44:                temp = mult(c, den);
        !            45:                *r = msub(num, temp);
        !            46:                if(mcmp(*r,den) >= 0){
        !            47:                        c = madd(c, one);
        !            48:                        *r = msub(*r, den);
        !            49:                }
        !            50:                if(mcmp(*r,zero) < 0){
        !            51:                        c = msub(c,one);
        !            52:                        *r = madd(*r, den);
        !            53:                }
        !            54:                if(nsign == -1){
        !            55:                        c.low = -c.low;
        !            56:                        r->high = -r->high;
        !            57:                        r->low = -r->low;
        !            58:                }
        !            59:                if(dsign == -1){
        !            60:                        c.low = -c.low;
        !            61:                }
        !            62:        return(c);
        !            63:        }
        !            64:        modf( tquot/e16 , &trial.high);
        !            65:        trial.low = tquot - trial.high*e16;
        !            66:        temp = mult(trial, den);
        !            67:        temp = msub(num, temp);
        !            68:        quot = mdiv(temp, den, &rem);
        !            69:        c = madd(trial, quot);
        !            70:        if(mcmp(rem,den) >= 0){
        !            71:                rem = msub(rem, den);
        !            72:                c = madd(c, one);
        !            73:        }
        !            74:        if(mcmp(rem,zero) < 0){
        !            75:                rem = madd(rem, den);
        !            76:                c = msub(c, one);
        !            77:        }
        !            78:        *r = rem;
        !            79:        if(nsign == -1){
        !            80:                c.high = -c.high;
        !            81:                c.low = -c.low;
        !            82:                r->high = -r->high;
        !            83:                r->low = -r->low;
        !            84:        }
        !            85:        if(dsign == -1){
        !            86:                c.high = -c.high;
        !            87:                c.low = -c.low;
        !            88:        }
        !            89:        return(c);
        !            90: }
        !            91: 
        !            92: mint
        !            93: sdiv(a,b,r)
        !            94: mint a;
        !            95: int b, *r;
        !            96: {
        !            97:        mint c;
        !            98:        mint bb, rr;
        !            99:        double temp1;
        !           100: 
        !           101:        if(a.high==0){
        !           102:                c.high = 0.;
        !           103:                modf(a.low/b , &temp1);
        !           104:                *r = a.low - temp1*b;
        !           105:                c.low = temp1;
        !           106:                return(c);
        !           107:        }
        !           108:        bb = itom(b);
        !           109:        c = mdiv(a,bb,&rr);
        !           110:        *r = rr.low;
        !           111:        return(c);
        !           112: }
        !           113: 
        !           114: mint
        !           115: msqrt(a,r)
        !           116: mint a, *r;
        !           117: {
        !           118:        mint b;
        !           119:        double temp, sqrt();
        !           120:        mint dtemp;
        !           121: 
        !           122:        if((a.high<0) || (a.low<0)){
        !           123:                printf("msqrt: Negative argument.\n");
        !           124:                abort();
        !           125:        }
        !           126:        temp = sqrt(a.high*e16+a.low);
        !           127:        modf(temp , &temp);
        !           128:        b.high = 0.;
        !           129:        b.low = temp;
        !           130:        dtemp.high = 0.;
        !           131:        dtemp.low = temp;
        !           132:        dtemp = mult(dtemp, dtemp);
        !           133:        *r = msub(a, dtemp);
        !           134:        return(b);
        !           135: }
        !           136: 
        !           137: mint
        !           138: mcbrt(a,r)
        !           139: mint a, *r;
        !           140: {
        !           141:        mint b;
        !           142:        double temp;
        !           143:        double log(), exp();
        !           144:        mint dtemp;
        !           145:        int sign = 1;
        !           146:        mint zero;
        !           147: 
        !           148:        zero.low = 0.0;
        !           149:        zero.high = 0.0;
        !           150:        if((a.high<0) || (a.low<0)){
        !           151:                sign = -1;
        !           152:                a = mchs(a);
        !           153:        }
        !           154: 
        !           155:        temp = exp(log(a.high*e16 + a.low)/3.0);
        !           156:        modf(temp, &temp);
        !           157: retry:
        !           158:        b.high = 0;
        !           159:        b.low = temp;
        !           160: 
        !           161:        dtemp.high = 0;
        !           162:        dtemp.low = temp;
        !           163:        dtemp = mult(dtemp,mult(dtemp,dtemp));
        !           164:        *r = msub(a, dtemp);
        !           165: 
        !           166:        if(mcmp(*r, zero) < 0){
        !           167:                temp = temp - 1;
        !           168:                goto retry;
        !           169:        }
        !           170:        temp = temp + 1;
        !           171:        dtemp.high = 0;
        !           172:        dtemp.low = temp;
        !           173:        dtemp = mult(dtemp,mult(dtemp,dtemp));
        !           174:        if(mcmp(a,dtemp) >= 0){
        !           175:                goto retry;
        !           176:        }
        !           177:        return(b);
        !           178: }

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