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

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.