Return to lfactor.c CVS log

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

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

1.1     ! root        1: #include <stdio.h>
        !             2: #include "mp.h"
        !             3: #include "form.h"
        !             4: #define testing 0
        !             5: #define fmax 20
        !             6: #define gridmax 30000
        !             7: form stockf[fmax];
        !             8: struct grid { short a; float b;};
        !             9: union{
        !            10:        struct grid grida[gridmax+1];
        !            11:        form powers[512];
        !            12: }ugly;
        !            13: FILE *ptab;
        !            14: 
        !            15: mint factab[100];
        !            16: int facptr;
        !            17: double flimit;
        !            18: 
        !            19: main(argc, argp)
        !            20: int argc;
        !            21: char *argp[];
        !            22: {
        !            23:        form formt1;
        !            24:        form formtt;
        !            25:        form bigf, bigfi;
        !            26:        form smallf;
        !            27:        form formttt;
        !            28:        form regen;
        !            29:        unsigned giant;
        !            30:        struct grid *j1, *j2;
        !            31:        int ordind;
        !            32:        int prime;
        !            33:        mint nn;
        !            34:        mint nugget;
        !            35:        mint mjunk, mtemp;
        !            36:        mint mtemp1;
        !            37:        mint zero, one, two, four;
        !            38:        double a, b;
        !            39:        mint ma, mb, mc;
        !            40:        double temp;
        !            41:        mint det;
        !            42:        double lastpr;
        !            43:        double m,n;
        !            44:        double prod;
        !            45:        double hbar, h;
        !            46:        double pi = 3.141592653589793;
        !            47:        double ht, hr, it;
        !            48:        double logord;
        !            49:        double maxord;
        !            50:        double order[fmax];
        !            51:        double tt, ttt;
        !            52:        float ftt, fttt;
        !            53:        int baby;
        !            54:        int kbaby;
        !            55:        register j;
        !            56:        int i;
        !            57:        double k;
        !            58:        int goti, gotj, gotit;
        !            59:        int ttn;
        !            60:        int junk;
        !            61:        int uneq;
        !            62:        int jmax;
        !            63:        int fcount;
        !            64:        int parity;
        !            65:        double sqroot, fifth;
        !            66:        int verbose = 0;
        !            67:        int nofact = 0;
        !            68:        int compare();
        !            69:        int compar1();
        !            70:        double getpr();
        !            71:        double modf();
        !            72:        double log(), exp(), sqrt();
        !            73: 
        !            74:        setbuf(stdout,NULL);
        !            75:        verbose = 0;
        !            76:        if(argc > 1){
        !            77:                if(*argp[1] == '+'){
        !            78:                        verbose = 1;
        !            79:                }
        !            80:                if(*argp[1] == '-'){
        !            81:                        nofact = 1;
        !            82:                }
        !            83:        }
        !            84:        zero = itom(0);
        !            85:        one = itom(1);
        !            86:        two = itom(2);
        !            87:        four = itom(4);
        !            88:        facptr = 0;
        !            89: 
        !            90:        if(min(&nn) == EOF) exit(0);
        !            91:        if(mcmp(nn,zero) <= 0) exit(0);
        !            92:        while(1){
        !            93:                mtemp = sdiv(nn,2,&i);
        !            94:                if(i==0){
        !            95:                        printf("Prime factor: 2\n");
        !            96:                        nn = mtemp;
        !            97:                }else break;
        !            98:        }
        !            99: 
        !           100:        ptab = fopen("/usr/lib/ptab","r");
        !           101:        if(ptab == NULL){
        !           102:                printf("ptab?\n");
        !           103:                exit(1);
        !           104:        }
        !           105: retry:
        !           106:        getpr(-1);
        !           107:        lastpr = 2.;
        !           108:        if(mcmp(nn,itom(1))==0)
        !           109:                exit(0);
        !           110:        sqroot = sqrt(nn.high*1e16 + nn.low);
        !           111:        fifth = exp(0.2*log(nn.high*1e16+nn.low));
        !           112:        if(fifth < 1000.) fifth = 1000.;
        !           113:        sdiv(nn,4, &i);
        !           114:        if((i== 1) || (i== 2)){
        !           115:                parity = 0;
        !           116:                prod = 1.;
        !           117:        }
        !           118:        if(i== 3)
        !           119:                parity = 1;
        !           120:        sdiv(nn, 8, &junk);
        !           121:        if(junk == 3) prod = 2./3.;
        !           122:        if(junk == 7) prod = 2.;
        !           123:        fcount = 0;
        !           124:        while(1){
        !           125:                if((m = getpr(0)) < 0)
        !           126:                        break;
        !           127:                lastpr = m;
        !           128:                if(lastpr > sqroot){
        !           129:                        printf("Prime factor: %.0f\n", nn.low);
        !           130:                        exit(0);
        !           131:                }
        !           132:                if(lastpr > 10.*fifth)
        !           133:                        break;
        !           134:                modf(nn.high/m, &temp);
        !           135:                n = nn.high - m*temp;
        !           136:                modf(e16/m, &temp);
        !           137:                temp = e16 - m*temp;
        !           138:                n = n*temp;
        !           139:                modf(nn.low/m, &temp);
        !           140:                temp = nn.low - m*temp;
        !           141:                n = n + temp;
        !           142:                j = jacobi(-n,m);
        !           143:                if((j==0) && (nofact==0)){
        !           144:                        printf("Prime factor: %.0f\n", m);
        !           145:                        mtemp.high = 0;
        !           146:                        mtemp.low = m;
        !           147:                        nn = mdiv(nn, mtemp, &mtemp);
        !           148:                        if(mcmp(mtemp,zero)!=0) abort();
        !           149:                        goto retry;
        !           150:                }
        !           151:                a = m;
        !           152:                if((j==1) && (fcount<fmax)){
        !           153:                        det = mchs(nn);
        !           154:                        sdiv(det, 4, &i);
        !           155:                        if((i== -1) || (i== -2))
        !           156:                                det = smult(det, 4);
        !           157:                        for(b=parity; b<=a; b = b+2){
        !           158:                                ma.high = 0.;
        !           159:                                ma.low = a;
        !           160:                                mb.high = 0.;
        !           161:                                mb.low = b;
        !           162:                                mc.high = 0.;
        !           163:                                mc.low = b*b;
        !           164:                                mc = msub(mc, det);
        !           165:                                mtemp.high = 0.;
        !           166:                                mtemp.low = 4.*a;
        !           167:                                mc = mdiv(mc, mtemp, &mtemp);
        !           168:                                if(mcmp(mtemp,zero) == 0){
        !           169:                                        stockf[fcount].a = ma;
        !           170:                                        stockf[fcount].b = mb;
        !           171:                                        stockf[fcount].c = mc;
        !           172: #if testing
        !           173:                                        if(verbose){
        !           174:                                                formout(stockf[fcount]);
        !           175:                                                printf("\n");
        !           176:                                        }
        !           177: #endif
        !           178:                                        fcount++;
        !           179:                                }
        !           180:                        }
        !           181:                }
        !           182:                prod *= m/(m-j);
        !           183:        }
        !           184:        printf("Trying to factor: n = ");
        !           185:        mout(nn);
        !           186:        printf("Factor limit = %.0f\n", a);
        !           187:        flimit = a;
        !           188:        mtemp1 = mcbrt(nn, &mjunk);
        !           189:        if(mjunk.low == 0){
        !           190:                mout1(nn);
        !           191:                printf(" is a cube.\n");
        !           192:                exit(0);
        !           193:        }
        !           194:        mtemp1 = msqrt(nn, &mjunk);
        !           195:        if(mjunk.low == 0){
        !           196:                mout1(nn);
        !           197:                printf(" is a square.\n");
        !           198:                exit(0);
        !           199:        }
        !           200: 
        !           201:        prime = pseudo(nn, itom(2));
        !           202:        if(prime == 0) prime = pseudo(nn, itom(3));
        !           203:        if(prime == 0) prime = pseudo(nn, itom(5));
        !           204:        if(prime == 0) prime = pseudo(nn, itom(7));
        !           205:        if(prime == 0){
        !           206:                printf("n is pseudoprime (mod 2,3,5,7).\n");
        !           207:        }else{
        !           208:                printf("n is composite.\n");
        !           209:        }
        !           210: 
        !           211:        if(parity==0) prod *= 2.;
        !           212:        modf(prod*mtemp1.low/pi, &hbar);
        !           213: #if testing
        !           214:        printf("h~~ %.0f\n", hbar);
        !           215: #endif
        !           216: 
        !           217:        ugly.grida[0].a = 0;
        !           218:        ugly.grida[0].b = 1.;
        !           219:        if(parity == 1){
        !           220:                baby = 1;
        !           221:                kbaby = 0;
        !           222:        }else{
        !           223:                baby = 2;
        !           224:                kbaby = 0;
        !           225:        }
        !           226: 
        !           227:        sdiv(nn, 8, &i);
        !           228:        if(i == 1){
        !           229:                baby = 4;
        !           230:                kbaby = 0;
        !           231:        }else if((i == 5) && (prime != 0)){
        !           232:                baby = 4;
        !           233:                kbaby = 0;
        !           234:        }
        !           235:        if((parity == 1) && (prime != 0)){
        !           236:                baby = 2;
        !           237:                kbaby = 0;
        !           238:        }
        !           239: 
        !           240:        modf(hbar/baby, &temp);
        !           241:        hbar = baby*temp;
        !           242:        hbar = hbar + kbaby;
        !           243: 
        !           244:        smallf = formpow(stockf[0], (double)baby);
        !           245:        ugly.grida[1].a = 1;
        !           246:        ugly.grida[1].b = smallf.a.low;
        !           247:        formt1 = formpow(stockf[0], hbar);
        !           248: 
        !           249:        if(mcmp(formt1.a,one) == 0){
        !           250:                h = hbar;
        !           251:                goto found;
        !           252:        }
        !           253:        formtt = smallf;
        !           254:        ttn = 1;
        !           255:        temp = sqrt(0.1*hbar/(sqrt(fifth)*baby));
        !           256:        if(temp > gridmax) temp = gridmax;
        !           257:        giant = temp;
        !           258:        printf("Grid = %d by %d\n", giant, baby);
        !           259:        while(ttn <= giant){
        !           260:                if(formtt.a.low == formt1.a.low){
        !           261:                        if(formtt.b.low == formt1.b.low){
        !           262:                                h = hbar - ttn*baby;
        !           263:                                goto found;
        !           264:                        }
        !           265:                        if(formtt.b.low == -formt1.b.low){
        !           266:                                h = hbar + ttn*baby;
        !           267:                                goto found;
        !           268:                        }
        !           269:                }
        !           270:                formtt = compos(formtt, smallf);
        !           271:                ttn += 1;
        !           272:                ugly.grida[ttn].a = ttn;
        !           273:                ugly.grida[ttn].b = formtt.a.low;
        !           274:        }
        !           275:        qsort(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compare);
        !           276:        bigf = formpow(smallf, (double)(2*giant));
        !           277:        bigfi = bigf;
        !           278:        bigfi.b.low = -bigfi.b.low;
        !           279:        formtt = compos(bigf, formt1);
        !           280:        formttt = compos(bigfi, formt1);
        !           281:        for(k=1; k<fifth/2.; k=k+1){
        !           282:                tt = formtt.a.low;
        !           283:                ttt = formttt.a.low;
        !           284:                ftt = formtt.a.low;
        !           285:                fttt = formttt.a.low;
        !           286:                if(search(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compar1, &ftt, &j1, &j2) >=0){
        !           287: b1:
        !           288:                        j = j1->a;
        !           289:                        regen = formpow(smallf, (double)j);
        !           290:                        if(tt == regen.a.low){
        !           291:                                if(formtt.b.low == regen.b.low){
        !           292:                                        h = hbar + (2.*giant*k - j)*baby;
        !           293:                                        goto found;
        !           294:                                }else
        !           295:                                if(formtt.b.low == -regen.b.low){
        !           296:                                        h = hbar + (2.*giant*k + j)*baby;
        !           297:                                        goto found;
        !           298:                                }
        !           299:                        }
        !           300:                        if(j2 > j1+1){
        !           301:                                j1 = j1 + 1;
        !           302:                                goto b1;
        !           303:                        }
        !           304:                }
        !           305:                if(search(&ugly.grida[0].a, giant+1, sizeof(ugly.grida[0]), compar1, &fttt, &j1, &j2) >=0){
        !           306: b2:
        !           307:                        j = j1->a;
        !           308:                        regen = formpow(smallf, (double)j);
        !           309:                        if(ttt == regen.a.low){
        !           310:                                if(formttt.b.low == regen.b.low){
        !           311:                                        h = hbar - (2.*giant*k + j)*baby;
        !           312:                                        goto found;
        !           313:                                }else
        !           314:                                if(formttt.b.low == -regen.b.low){
        !           315:                                        h = hbar - (2.*giant*k - j)*baby;
        !           316:                                        goto found;
        !           317:                                }
        !           318:                        }
        !           319:                        if(j2 > j1 + 1){
        !           320:                                j1 = j1 + 1;
        !           321:                                goto b2;
        !           322:                        }
        !           323:                }
        !           324:                formtt = compos(formtt, bigf);
        !           325:                formttt = compos(formttt, bigfi);
        !           326:        }
        !           327:        printf("Search abandoned at %.0f.\n", 2.*giant*k*baby);
        !           328:        exit(1);
        !           329: found:
        !           330:        printf("h = %.0f\n", h);
        !           331:        temp = hbar - h;
        !           332:        printf("Search succeeded at %.0f (%.3f) ", temp, 1000.*temp/h);
        !           333:        printf("(%.3f)\n",flimit*temp*temp/(h*h));
        !           334: 
        !           335:        ht = 0;
        !           336:        hr =h;
        !           337:        while(1){
        !           338:                modf(hr/2, &temp);
        !           339:                if(hr == 2.*temp){
        !           340:                        hr = temp;
        !           341:                        ht += 1;
        !           342:                }else{
        !           343:                        break;
        !           344:                }
        !           345:        }
        !           346: /*
        !           347:  * h = hr * 2^ht
        !           348:  */
        !           349:        if(ht == 0){
        !           350:                printf("Class number is odd!\n");
        !           351:                mout1(nn);
        !           352:                printf(" is a prime.\n");
        !           353:                exit(0);
        !           354:        }
        !           355: #if testing
        !           356:        printf("Sylow 2-subgroup of order 2^%.0f\n", ht);
        !           357: #endif
        !           358:        for(i=0; i<fmax; i++){
        !           359:                stockf[i] = formpow(stockf[i], hr);
        !           360:        }
        !           361: 
        !           362:        for(i=0; i<fmax; i++){
        !           363:                if(mcmp(stockf[i].b,zero) < 0){
        !           364:                        stockf[i].b = mchs(stockf[i].b);
        !           365:                }
        !           366:        }
        !           367: 
        !           368:        for(i=1,jmax=1; i<fmax; i++){
        !           369:                for(j=0; j<jmax; j++){
        !           370:                        uneq = 0;
        !           371:                        if(mcmp(stockf[i].a,stockf[j].a) != 0) uneq = 1;
        !           372:                        if(mcmp(stockf[i].b,stockf[j].b) != 0) uneq = 1;
        !           373:                        if(mcmp(stockf[i].c,stockf[j].c) != 0) uneq = 1;
        !           374:                        if(uneq == 0) break;
        !           375:                }
        !           376:                if(uneq == 1){
        !           377:                        stockf[jmax++] = stockf[i];
        !           378:                }
        !           379:        }
        !           380: 
        !           381:        for(i=0; i<jmax; i++){
        !           382:                it = 0;
        !           383:                formt1 = stockf[i];
        !           384:                if(mcmp(formt1.a,one) == 0) continue;
        !           385:                while(1){
        !           386:                        formtt = formpow(formt1, (double)2);
        !           387:                        it = it+1;
        !           388:                        if(it>ht){
        !           389:                                printf("main: error no. 1\n");
        !           390:                                abort();
        !           391:                        }
        !           392:                        if(mcmp(formtt.a,one) == 0) break;
        !           393:                        formt1 = formtt;
        !           394:                }
        !           395:                order[i] = it;
        !           396: /*
        !           397:                printf("f%3d: ", i);
        !           398:                formout(stockf[i]);
        !           399:                printf(" order = 2^%.0f: ", it);
        !           400:                formout(formt1);
        !           401:                printf("\n");
        !           402: */
        !           403:        }
        !           404:        logord = 0;
        !           405:        ordind = 0;
        !           406:        for(i=0; i<jmax; i++){
        !           407:                if(order[i] > logord){
        !           408:                        ordind = i;
        !           409:                        logord = order[i];
        !           410:                }
        !           411:        }
        !           412: 
        !           413:        if(logord == ht){
        !           414:                printf("Sylow 2-subgroup is cyclic.\n");
        !           415:                formtt = stockf[ordind];
        !           416:                for(i=0; i<order[ordind]-1; i++){
        !           417:                        formtt = formpow(formtt, (double)2);
        !           418:                }
        !           419:                if(mcmp(formtt.a,two) == 0){
        !           420:                        mout1(nn);
        !           421:                        printf(" is a prime.\n");
        !           422:                        exit(0);
        !           423:                }
        !           424:                if(mcmp(formtt.b,zero) == 0){
        !           425:                        sqtest(formtt.a, flimit);
        !           426:                        sqtest(formtt.c, flimit);
        !           427:                        exit(0);
        !           428:                }
        !           429:                if(mcmp(formtt.a, formtt.b) == 0){
        !           430:                        sqtest(formtt.a, flimit);
        !           431:                        sqtest(msub(mult(itom(4),formtt.c),formtt.a), flimit);
        !           432:                        exit(0);
        !           433:                }
        !           434:                if(mcmp(formtt.a, formtt.c) == 0){
        !           435:                        sqtest(msub(mult(itom(2),formtt.a),formtt.b), flimit);
        !           436:                        sqtest(madd(mult(itom(2),formtt.a),formtt.b), flimit);
        !           437:                        exit(0);
        !           438:                }
        !           439:                printf("main: error no. 4\n");
        !           440:                exit(1);
        !           441:        }
        !           442:        if(logord > 10){
        !           443:                printf("Sylow 2-subgroup too big.\n");
        !           444:                exit(1);
        !           445:        }
        !           446: 
        !           447:        for(i=0,j=1;i<logord;i++){
        !           448:                j = j*2;
        !           449:        }
        !           450:        maxord = j;
        !           451:        printf("\n");
        !           452: 
        !           453:        ugly.powers[0] = stockf[ordind];
        !           454:        for(i=1;i<maxord/2;i++){
        !           455:                ugly.powers[i] = compos(ugly.powers[i-1],stockf[ordind]);
        !           456:        }
        !           457: 
        !           458:        for(i=0; i<jmax; i++){
        !           459:                if(i == ordind) continue;
        !           460:                gotit = 0;
        !           461:                formt1 = stockf[i];
        !           462:                for(j=0; j<=order[i]; j++){
        !           463:                        for(junk=0; junk<maxord/2; junk++){
        !           464:                                if(mcmp(formt1.a,ugly.powers[junk].a) != 0)
        !           465:                                        continue;
        !           466:                                if(mcmp(formt1.c,ugly.powers[junk].c) != 0)
        !           467:                                        continue;
        !           468:                                if(mcmp(formt1.b,ugly.powers[junk].b) != 0){
        !           469:                                        goti = junk + 1;
        !           470:                                }else{
        !           471:                                        goti = maxord - junk - 1;
        !           472:                                }
        !           473:                                for(junk=0,gotj=1; junk<j; junk++){
        !           474:                                        gotj *= 2;
        !           475:                                }
        !           476:                                if(goti%gotj != 0){
        !           477:                                        printf("main: error no. 2\n");
        !           478:                                }
        !           479:                                goti = goti/gotj;
        !           480:                                if(goti > maxord/2){
        !           481:                                        goti = maxord - goti;
        !           482:                                }
        !           483:                                stockf[i] = compos(stockf[i],ugly.powers[goti-1]);
        !           484:                                gotit = 1;
        !           485:                                break;
        !           486:                        }
        !           487:                        if(gotit == 1) break;
        !           488:                        formt1 = compos(formt1,formt1);
        !           489:                }
        !           490:        }
        !           491: 
        !           492:        for(i=0; i<jmax; i++){
        !           493:                if(mcmp(stockf[i].b,zero) < 0){
        !           494:                        stockf[i].b = mchs(stockf[i].b);
        !           495:                }
        !           496:        }
        !           497: 
        !           498:        junk = jmax;
        !           499:        for(i=1,jmax=1; i<junk; i++){
        !           500:                for(j=0; j<jmax; j++){
        !           501:                        uneq = 0;
        !           502:                        if(mcmp(stockf[i].a,stockf[j].a) != 0) uneq = 1;
        !           503:                        if(mcmp(stockf[i].b,stockf[j].b) != 0) uneq = 1;
        !           504:                        if(mcmp(stockf[i].c,stockf[j].c) != 0) uneq = 1;
        !           505:                        if(uneq == 0) break;
        !           506:                }
        !           507:                if(uneq == 1){
        !           508:                        stockf[jmax++] = stockf[i];
        !           509:                }
        !           510:        }
        !           511: 
        !           512:        for(i=0; i<jmax; i++){
        !           513:                it = 0;
        !           514:                formt1 = stockf[i];
        !           515:                if(mcmp(formt1.a,one) == 0) continue;
        !           516:                while(1){
        !           517:                        formtt = formpow(formt1, (double)2);
        !           518:                        it = it+1;
        !           519:                        if(it>ht){
        !           520:                                printf("main: error no. 3\n");
        !           521:                                abort();
        !           522:                        }
        !           523:                        if(mcmp(formtt.a,one) == 0) break;
        !           524:                        formt1 = formtt;
        !           525:                }
        !           526:                order[i] = it;
        !           527:                if(verbose){
        !           528:                        printf("f%3d: ", i);
        !           529:                        formout(stockf[i]);
        !           530:                        printf(" order = 2^%.0f: ", it);
        !           531:                        formout(formt1);
        !           532:                        printf("\n");
        !           533:                }
        !           534:                if(mcmp(formt1.b,zero) == 0){
        !           535:                        putfact(formt1.a);
        !           536:                        putfact(formt1.c);
        !           537:                }else if(mcmp(formt1.a,formt1.b) == 0){
        !           538:                        putfact(formt1.a);
        !           539:                        putfact(msub(mult(four,formt1.c),formt1.a));
        !           540:                }else if(mcmp(formt1.a,formt1.c) == 0){
        !           541:                        putfact(madd(mult(two,formt1.a),formt1.b));
        !           542:                        putfact(msub(mult(two,formt1.a),formt1.b));
        !           543:                }else{
        !           544:                        printf("main: error no. 5\n");
        !           545:                }
        !           546:        }
        !           547:        nugget = nn;
        !           548:        for(i=0;i<facptr;i++){
        !           549:                if(mcmp(factab[i],one) == 0) continue;
        !           550:                if(mcmp(factab[i],nn) >= 0) continue;
        !           551:                if(mcmp(factab[i],nugget) > 0) continue;
        !           552:                prtest(factab[i], flimit);
        !           553:                nugget = mdiv(nugget, factab[i], &mjunk);
        !           554:                if(mcmp(mjunk,zero) != 0){
        !           555:                        printf("main: error no. 6\n");
        !           556:                }
        !           557:        }
        !           558:        if(mcmp(nugget,one) > 0){
        !           559:                prtest(nugget, flimit);
        !           560:        }
        !           561:        exit(0);
        !           562: }
        !           563: 
        !           564: /*
        !           565:  * m must be positive and odd.
        !           566:  */
        !           567: int
        !           568: jacobi(n,m)
        !           569: double n,m;
        !           570: {
        !           571:        int mr2, mrm1, i, j, sign;
        !           572:        double temp;
        !           573:        double modf();
        !           574: 
        !           575:        sign = 1;
        !           576:        if(n==0) return(0);
        !           577:        if(n==1) return(1);
        !           578:        if(m==1) return(1);
        !           579: 
        !           580:        mr2 = 1;
        !           581:        mrm1 = 1;
        !           582:        modf(m/8, &temp);
        !           583:        i = m - 8*temp;
        !           584:        switch(i){
        !           585:        case 3:
        !           586:                mrm1 = -1;
        !           587:        case 5:
        !           588:                mr2 = -1;
        !           589:        case 1:
        !           590:                break;
        !           591:        case 7:
        !           592:                mrm1 = -1;
        !           593:                break;
        !           594:        default:
        !           595:                abort();
        !           596:        }
        !           597:        
        !           598:        modf(n/m, &temp);
        !           599:        n = n - m*temp;
        !           600:        if(n<0) n += m;
        !           601:        if(n==0) return(0);
        !           602:        if(n+n>m){
        !           603:                n = m-n;
        !           604:                sign = sign*mrm1;
        !           605:        }
        !           606:        while(1){
        !           607:                if(modf(.5*n, &temp) != 0) break;
        !           608:                n = temp;
        !           609:                sign = sign * mr2;
        !           610:        }
        !           611:        if(m>=32768.)
        !           612:                return(sign * jacobi(mrm1*m,n));
        !           613:        else{
        !           614:                i = mrm1*m;
        !           615:                j = n;
        !           616:                return(sign * ijac(i,j));
        !           617:        }
        !           618: }
        !           619: /*
        !           620:  * m must be positive and odd.
        !           621:  */
        !           622: int
        !           623: ijac(n,m)
        !           624: int n,m;
        !           625: {
        !           626:        int mr2, mrm1, sign;
        !           627: 
        !           628:        sign = 1;
        !           629:        if(n==0) return(0);
        !           630:        if(n==1) return(1);
        !           631:        if(m==1) return(1);
        !           632: 
        !           633:        mr2 = mrm1 = 1;
        !           634:        switch(m&07){
        !           635:        case 3:
        !           636:                mrm1 = -1;
        !           637:        case 5:
        !           638:                mr2 = -1;
        !           639:        case 1:
        !           640:                break;
        !           641:        case 7:
        !           642:                mrm1 = -1;
        !           643:                break;
        !           644:        default:
        !           645:                abort();
        !           646:        }
        !           647:        
        !           648:        n = n%m;
        !           649:        if(n<0) n += m;
        !           650:        if(n==0) return(0);
        !           651:        if((n&01)==1){
        !           652:                n = m-n;
        !           653:                sign = sign*mrm1;
        !           654:        }
        !           655:        while((n&03) == 0)
        !           656:                n = n>>2;
        !           657:        if((n&01) == 0){
        !           658:                n = n>>1;
        !           659:                sign = sign * mr2;
        !           660:        }
        !           661:        return(sign * ijac(mrm1*m,n));
        !           662: }
        !           663: 
        !           664: int sieve[] = {
        !           665:         1,  7, 11, 13, 17, 19, 23, 29,
        !           666:        31, 37, 41, 43, 47, 49, 53, 59,
        !           667:        61, 67, 71, 73, 77, 79, 83, 89,
        !           668:        91, 97,101,103,107,109,113,119,
        !           669: };
        !           670: double
        !           671: getpr(arg)
        !           672: int arg;
        !           673: {
        !           674:        static int i;
        !           675:        static unsigned word;
        !           676:        static double base;
        !           677: 
        !           678:        if(arg<0){
        !           679:                rewind(ptab);
        !           680:                word = getw(ptab);
        !           681:                if(feof(ptab)) exit(0);
        !           682:                if(ferror(ptab)) exit(0);
        !           683:                base = 0;
        !           684:                i = -2;
        !           685:                return(0.);
        !           686:        }
        !           687:        if((base == 0) && (i < 0)){
        !           688:                if(i == -2){
        !           689:                        i++;
        !           690:                        return(3);
        !           691:                }
        !           692:                if(i == -1){
        !           693:                        i++;
        !           694:                        return(5);
        !           695:                }
        !           696:        }
        !           697:        while(1){
        !           698:                i++;
        !           699:                if(i>(8*sizeof(int))){
        !           700:                        word = getw(ptab);
        !           701:                        if(ferror(ptab)) abort();
        !           702:                        if(feof(ptab)) return((double)-1);
        !           703:                        i -= 8*sizeof(int);
        !           704:                        base += 30*sizeof(int);
        !           705:                }
        !           706:                if((word & (1<<(i-1))) != 0)
        !           707:                        return(base+sieve[i-1]);
        !           708:        }
        !           709: }
        !           710: 
        !           711: int
        !           712: compare(a,b)
        !           713: struct grid *a, *b;
        !           714: {
        !           715:        if(a->b > b->b)
        !           716:                return(1);
        !           717:        if(a->b < b->b)
        !           718:                return(-1);
        !           719:        return(0);
        !           720: }
        !           721: int
        !           722: compar1(a,b)
        !           723: float *a;
        !           724: struct grid *b;
        !           725: {
        !           726:        if(*a > b->b)
        !           727:                return(1);
        !           728:        if(*a < b->b)
        !           729:                return(-1);
        !           730:        return(0);
        !           731: }
        !           732: 
        !           733: void
        !           734: sqtest(nn, flimit)
        !           735: mint nn;
        !           736: double flimit;
        !           737: {
        !           738: 
        !           739:        mint temp1, temp2;
        !           740:        double dtemp;
        !           741: 
        !           742:        temp1.low = flimit;
        !           743:        temp1.high = 0;
        !           744:        if(mcmp(temp1,nn) > 0){
        !           745:                printf("sqtest: error no. 1\n");
        !           746:                exit(1);
        !           747:        }
        !           748: 
        !           749:        dtemp = log(nn.high+nn.low)/log(flimit);
        !           750:        if(dtemp >= 5.0){
        !           751:                printf("sqtest: error no. 2");
        !           752:                exit(1);
        !           753:        }
        !           754: 
        !           755:        temp1 = msqrt(nn, &temp2);
        !           756:        if(temp2.low == 0){
        !           757:                printf("Square factor: ");
        !           758:                mout(temp1);
        !           759:                return;
        !           760:        }
        !           761:        temp1 = mcbrt(nn, &temp2);
        !           762:        if(temp2.low == 0){
        !           763:                printf("Cubic factor: ");
        !           764:                mout(temp1);
        !           765:                return;
        !           766:        }
        !           767:                printf("Prime factor: ");
        !           768:                mout(nn);
        !           769:        return;
        !           770: }
        !           771: putfact(arg)
        !           772: mint arg;
        !           773: {
        !           774:        mint mtemp;
        !           775:        int i;
        !           776: 
        !           777: 
        !           778:        while(1){
        !           779:                mtemp = sdiv(arg,2,&i);
        !           780:                if(i==0){
        !           781:                        arg = mtemp;
        !           782:                }
        !           783:                else break;
        !           784:        }
        !           785:        if(facptr == 0){
        !           786:                factab[facptr++] = arg;
        !           787:                return(0);
        !           788:        }
        !           789:        for(i=0; i<facptr; i++){
        !           790:                if(mcmp(factab[i],arg) < 0) continue;
        !           791:                if(mcmp(factab[i],arg) == 0) return(0);
        !           792:                mtemp = factab[i];
        !           793:                factab[i] = arg;
        !           794:                arg = mtemp;
        !           795:                continue;
        !           796:        }
        !           797:        factab[facptr++] = arg;
        !           798:        return;
        !           799: 
        !           800: }
        !           801: 
        !           802: prtest(arg)
        !           803: mint arg;
        !           804: {
        !           805:        double dtemp;
        !           806: 
        !           807:        dtemp = log(arg.high+arg.low)/log(flimit);
        !           808:        if(dtemp < 2.0){
        !           809:                printf("Prime factor: ");
        !           810:        }else{
        !           811:                printf("Factor: ");
        !           812:        }
        !           813:        mout(arg);
        !           814:        return;
        !           815: }

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