Return to mulmod.c CVS log

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

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

1.1     ! root        1: #include "mp.h"
        !             2: mint mulmod();
        !             3: 
        !             4: /*
        !             5:        returns 0 if pseudoprime
        !             6:                1 if composite
        !             7:                2 if composite pseudoprime
        !             8: */
        !             9: int
        !            10: pseudo(nn, arg2)
        !            11: mint nn;
        !            12: mint arg2;
        !            13: {
        !            14:        mint nm1, y, z;
        !            15:        mint zero, one, two;
        !            16:        mint rem;
        !            17:        mint mtemp;
        !            18:        mint lasty;
        !            19:        int i;
        !            20: 
        !            21:        zero.low = zero.high = 0.0;
        !            22:        one = itom(1);
        !            23:        two = itom(2);
        !            24: 
        !            25:                if(mcmp(nn,zero) <= 0) abort(0);
        !            26:                if(mcmp(nn,two) <=0) abort(0);
        !            27:                mdiv(nn, two, &rem);
        !            28:                if(mcmp(rem, zero) == 0.0) abort(0);
        !            29: 
        !            30:                nm1 = msub(nn, one);
        !            31:                y = one;
        !            32:                z = arg2;
        !            33: 
        !            34:                i = 0;
        !            35:                while(1){
        !            36:                        mtemp = mdiv(nm1, two, &rem);
        !            37:                        if(mcmp(rem,zero) == 0){
        !            38:                                nm1 = mtemp;
        !            39:                                i = i + 1;
        !            40:                        }else break;
        !            41:                }
        !            42: 
        !            43:                while(mcmp(nm1,zero) != 0){
        !            44:                        nm1 = mdiv(nm1, two, &rem);
        !            45:                        if(mcmp(rem, one) == 0){
        !            46:                                y = mulmod(y, z, nn);
        !            47:                        }
        !            48:                        z = mulmod(z, z, nn);
        !            49:                }
        !            50:                if(mcmp(y, one) == 0){
        !            51:                        return(0);
        !            52:                }
        !            53:                while(i--){
        !            54:                        lasty = y;
        !            55:                        y = mulmod(y,y,nn);
        !            56:                        if(mcmp(y, one) == 0) break;
        !            57:                }
        !            58:                if(mcmp(y, one) != 0){
        !            59:                        return(1);
        !            60:                }
        !            61:                if(mcmp(lasty, one) == 0){
        !            62:                        printf("pseudo: error.\n");
        !            63:                }
        !            64:                if(mcmp(lasty,msub(nn,one)) != 0){
        !            65:                        return(2);
        !            66:                }
        !            67:                return(0);
        !            68: }
        !            69: mint
        !            70: mulmod(a,b,c)
        !            71: mint a, b, c;
        !            72: {
        !            73: 
        !            74:        mint mjunk;
        !            75:        mint mtemp1, mtemp2;
        !            76:        mint xy0, xy1, xy2, xy3;
        !            77:        mint x1, x2, y1, y2;
        !            78:        int i;
        !            79:        double z0, z1, z2, z3;
        !            80:        double s0, s1, s2;
        !            81:        double tquot, dtemp1, dtemp2;
        !            82:        mint zero;
        !            83: 
        !            84: 
        !            85:        zero.low = 0.0;
        !            86:        zero.high = 0.0;
        !            87:        if(mcmp(c,zero) == 0){
        !            88:                printf("mulmod: divide check\n");
        !            89:                abort(0);
        !            90:        }
        !            91: 
        !            92:        if((mcmp(a,zero)<0) || (mcmp(b,zero)<0)){
        !            93:                printf("mulmod: negative argument.\n");
        !            94:                abort(0);
        !            95:        }
        !            96:        if(mcmp(c,zero) < 0){
        !            97:                printf("mulmod: negative modulus.\n");
        !            98:                abort(0);
        !            99:        }
        !           100:        x1.low = a.low;
        !           101:        x1.high = 0;
        !           102:        x2.low = a.high;
        !           103:        x2.high = 0;
        !           104: 
        !           105:        y1.low = b.low;
        !           106:        y1.high = 0;
        !           107:        y2.low = b.high;
        !           108:        y2.high = 0;
        !           109: 
        !           110:        xy0 = mult(x1,y1);
        !           111:        xy1 = mult(x1,y2);
        !           112:        xy2 = mult(x2,y1);
        !           113:        xy3 = mult(x2,y2);
        !           114: 
        !           115:        z0 = xy0.low;
        !           116: 
        !           117:        mtemp1.low = xy0.high;
        !           118:        mtemp1.high = 0;
        !           119:        mtemp2.low = xy1.low;
        !           120:        mtemp2.high = 0;
        !           121:        mtemp1 = madd(mtemp1, mtemp2);
        !           122:        mtemp2.low = xy2.low;
        !           123:        mtemp1 = madd(mtemp1, mtemp2);
        !           124:        z1 = mtemp1.low;
        !           125: 
        !           126:        mtemp1.low = mtemp1.high;
        !           127:        mtemp1.high = 0;
        !           128:        mtemp2.low = xy1.high;
        !           129:        mtemp2.high = 0;
        !           130:        mtemp1 = madd(mtemp1, mtemp2);
        !           131:        mtemp2.low = xy2.high;
        !           132:        mtemp1 = madd(mtemp1, mtemp2);
        !           133:        mtemp2.low = xy3.low;
        !           134:        mtemp1 = madd(mtemp1, mtemp2);
        !           135:        z2 = mtemp1.low;
        !           136: 
        !           137:        mtemp1.low = mtemp1.high;
        !           138:        mtemp1.high = 0;
        !           139:        mtemp2.low = xy3.high;
        !           140:        mtemp2.high = 0;
        !           141:        mtemp1 = madd(mtemp1, mtemp2);
        !           142:        z3 = mtemp1.low;
        !           143: 
        !           144: 
        !           145:        if(mtemp1.high != 0.0){
        !           146:                printf("mulmod: error\n");
        !           147:        }
        !           148: 
        !           149:        if(c.high == 0.0){
        !           150:                mtemp1.high = 0.0;
        !           151:                mtemp1.low = z3;
        !           152:                mjunk = mdiv(mtemp1, c, &mtemp2);
        !           153:                z3 = mtemp2.low;
        !           154:                if(mtemp2.high != 0.0) trouble(12);
        !           155:                mtemp1.high = z3;
        !           156:                mtemp1.low = z2;
        !           157:                mjunk = mdiv(mtemp1, c, &mtemp2);
        !           158:                mtemp1.high = mtemp2.low;
        !           159:                mtemp1.low = z1;
        !           160:                mjunk = mdiv(mtemp1, c, &mtemp2);
        !           161:                mtemp1.high = mtemp2.low;
        !           162:                mtemp1.low = z0;
        !           163:                mjunk = mdiv(mtemp1, c, &mtemp2);
        !           164:                if(mcmp(mtemp2, c) >= 0) trouble(10);
        !           165:                if(mcmp(mtemp2, zero) < 0) trouble(11);
        !           166:                return(mtemp2);
        !           167:        }
        !           168: 
        !           169:        mtemp1.high = z3;
        !           170:        mtemp1.low = z2;
        !           171:        if(mcmp(mtemp1, c) >= 0){
        !           172:                mjunk = mdiv(mtemp1, c, &mtemp2);
        !           173:                z3 = mtemp2.high;
        !           174:                z2 = mtemp2.low;
        !           175:        }
        !           176:        dtemp1 = (z3*e16 + z2) + z1/e16;
        !           177:        tquot = dtemp1/(c.high+c.low/e16);
        !           178:        modf(tquot, &tquot);
        !           179:        y1.low = tquot;
        !           180:        y1.high = 0.0;
        !           181:        x1.low = c.low;
        !           182:        x1.high = 0.0;
        !           183:        x2.low = c.high;
        !           184:        x2.high = 0;
        !           185:        xy0 = mult(x1, y1);
        !           186:        xy1 = mult(x2, y1);
        !           187: 
        !           188:        s0 = xy0.low;
        !           189: 
        !           190:        mtemp1.low = xy0.high;
        !           191:        mtemp1.high = 0.0;
        !           192:        mtemp1 = madd(mtemp1, xy1);
        !           193:        s1 = mtemp1.low;
        !           194:        s2 = mtemp1.high;
        !           195:        s0 = z1 - s0;
        !           196:        s1 = z2 - s1;
        !           197:        s2 = z3 - s2;
        !           198: 
        !           199:        if(s2 > 0.0) {s2 = s2 - 1; s1 = s1 + e16;}
        !           200:        if(s2 < 0.0) {s2 = s2 + 1; s1 = s1 - e16;}
        !           201:        if((s1<0.0) && (s0>0.0)){
        !           202:                s1 = s1 + 1;
        !           203:                s0 = s0 - e16;
        !           204:        }else
        !           205:        if((s1>0.0) && (s0<0.0)){
        !           206:                s1 = s1 - 1;
        !           207:                s0 = s0 + e16;
        !           208:        }
        !           209:        if((s1*e16+s0) < 0.0){
        !           210:                mtemp1.low = s0;
        !           211:                mtemp1.high = s1;
        !           212:                mtemp1 = madd(mtemp1, c);
        !           213:                s1 = mtemp1.high;
        !           214:                s0 = mtemp1.low;
        !           215:        }
        !           216:        if(s2 != 0.0) trouble(1);
        !           217:        mtemp1.low = s0;
        !           218:        mtemp1.high = s1;
        !           219:        if(mcmp(mtemp1,zero) < 0) trouble(2);
        !           220:        if(mcmp(mtemp1,c) >= 0){
        !           221:                mtemp1.high = s1;
        !           222:                mtemp1.low = s0;
        !           223:                mtemp1 = msub(mtemp1, c);
        !           224:                s1 = mtemp1.high;
        !           225:                s0 = mtemp1.low;
        !           226:        }
        !           227:        if(mcmp(mtemp1,c) >0) trouble(3);
        !           228: 
        !           229:        z3 = s2;
        !           230:        z2 = s1;
        !           231:        z1 = s0;
        !           232: 
        !           233:        dtemp1 = (z2*e16 + z1) + z0/e16;
        !           234:        tquot = dtemp1/(c.high + c.low/e16);
        !           235:        modf(tquot, &tquot);
        !           236:        y1.low = tquot;
        !           237:        y1.high = 0.0;
        !           238:        x1.low = c.low;
        !           239:        x1.high = 0.0;
        !           240:        x2.low = c.high;
        !           241:        x2.high = 0.0;
        !           242:        xy0 = mult(x1, y1);
        !           243:        xy1 = mult(x2, y1);
        !           244: 
        !           245:        s0 = xy0.low;
        !           246: 
        !           247:        mtemp1.low = xy0.high;
        !           248:        mtemp1.high = 0.0;
        !           249:        mtemp1 = madd(mtemp1, xy1);
        !           250:        s1 = mtemp1.low;
        !           251:        s2 = mtemp1.high;
        !           252: 
        !           253:        s0 = z0 - s0;
        !           254:        s1 = z1 - s1;
        !           255:        s2 = z2 - s2;
        !           256:        if(s2 > 0.0) {s2 = s2 - 1; s1 = s1 + e16;}
        !           257:        if(s2 < 0.0) {s2 = s2 + 1; s1 = s1 - e16;}
        !           258:        if((s1<0.0) && (s0>0.0)){
        !           259:                s1 = s1 + 1;
        !           260:                s0 = s0 - e16;
        !           261:        }else
        !           262:        if((s1>0.0) && (s0<0.0)){
        !           263:                s1 = s1 - 1;
        !           264:                s0 = s0 + e16;
        !           265:        }
        !           266:        if((s1*e16+s0) < 0.0){
        !           267:                mtemp1.low = s0;
        !           268:                mtemp1.high = s1;
        !           269:                mtemp1 = madd(mtemp1, c);
        !           270:                s1 = mtemp1.high;
        !           271:                s0 = mtemp1.low;
        !           272:        }
        !           273:        if(s2 != 0.0) trouble(4);
        !           274:        mtemp1.low = s0;
        !           275:        mtemp1.high = s1;
        !           276:        if(mcmp(mtemp1,zero) < 0) trouble(5);
        !           277:        if(mcmp(mtemp1,c) >= 0){
        !           278:                mtemp1.high = s1;
        !           279:                mtemp1.low = s0;
        !           280:                mtemp1 = msub(mtemp1, c);
        !           281:                s1 = mtemp1.high;
        !           282:                s0 = mtemp1.low;
        !           283:        }
        !           284:        if(mcmp(mtemp1,c) > 0) trouble(6);
        !           285: 
        !           286:        z2 = s2;
        !           287:        z1 = s1;
        !           288:        z0 = s0;
        !           289:        mtemp1.high = s1;
        !           290:        mtemp1.low = s0;
        !           291:        return(mtemp1);
        !           292: }
        !           293: 
        !           294: trouble(i)
        !           295: int i;
        !           296: {
        !           297:        printf("mulmod: error %d\n", i);
        !           298:        abort(0);
        !           299: }

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