[BACK]Return to fctr.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/fctr.c, Revision 1.1

1.1     ! noro        1: /* $OpenXM: OpenXM/src/asir99/builtin/fctr.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
        !             2: #include "ca.h"
        !             3: #include "parse.h"
        !             4:
        !             5: void Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint();
        !             6: void Pptozp(), Pcont();
        !             7: void Pafctr(), Pagcd();
        !             8: void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab();
        !             9: void Pirred_check(), Pnfctr_mod();
        !            10:
        !            11: struct ftab fctr_tab[] = {
        !            12:        {"fctr",Pfctr,1},
        !            13:        {"gcd",Pgcd,-3},
        !            14:        {"gcdz",Pgcdz,2},
        !            15:        {"lcm",Plcm,2},
        !            16:        {"sqfr",Psqfr,1},
        !            17:        {"ufctrhint",Pufctrhint,2},
        !            18:        {"ptozp",Pptozp,1},
        !            19:        {"cont",Pcont,-2},
        !            20:        {"afctr",Pafctr,2},
        !            21:        {"agcd",Pagcd,3},
        !            22:        {"modsqfr",Pmodsqfr,2},
        !            23:        {"modfctr",Pmodfctr,2},
        !            24: #if 0
        !            25:        {"ddd",Pddd,2},
        !            26:        {"newddd",Pnewddd,2},
        !            27: #endif
        !            28:        {"ddd_tab",Pddd_tab,2},
        !            29:        {"irred_check",Pirred_check,2},
        !            30:        {"nfctr_mod",Pnfctr_mod,2},
        !            31:        {0,0,0},
        !            32: };
        !            33:
        !            34: void Pfctr(arg,rp)
        !            35: NODE arg;
        !            36: LIST *rp;
        !            37: {
        !            38:        DCP dc;
        !            39:
        !            40:        asir_assert(ARG0(arg),O_P,"fctr");
        !            41:        fctrp(CO,(P)ARG0(arg),&dc);
        !            42:        dcptolist(dc,rp);
        !            43: }
        !            44:
        !            45: void Pgcd(arg,rp)
        !            46: NODE arg;
        !            47: P *rp;
        !            48: {
        !            49:        P p1,p2,g1,g2,g;
        !            50:        Num m;
        !            51:        int mod;
        !            52:
        !            53:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
        !            54:        asir_assert(p1,O_P,"gcd");
        !            55:        asir_assert(p2,O_P,"gcd");
        !            56:        if ( !p1 )
        !            57:                *rp = p2;
        !            58:        else if ( !p2 )
        !            59:                *rp = p1;
        !            60:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
        !            61:                error("gcd : invalid argument");
        !            62:        else if ( argc(arg) == 2 )
        !            63:                ezgcdp(CO,p1,p2,rp);
        !            64:        else {
        !            65:                m = (Num)ARG2(arg);
        !            66:                asir_assert(m,O_P,"gcd");
        !            67:                mod = QTOS((Q)m);
        !            68:                ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
        !            69:                gcdprsmp(CO,mod,g1,g2,&g);
        !            70:                mptop(g,rp);
        !            71:        }
        !            72: }
        !            73:
        !            74: void Pgcdz(arg,rp)
        !            75: NODE arg;
        !            76: P *rp;
        !            77: {
        !            78:        P p1,p2,t;
        !            79:        Q c1,c2;
        !            80:        N n;
        !            81:
        !            82:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
        !            83:        asir_assert(p1,O_P,"gcdz");
        !            84:        asir_assert(p2,O_P,"gcdz");
        !            85:        if ( !p1 )
        !            86:                *rp = p2;
        !            87:        else if ( !p2 )
        !            88:                *rp = p1;
        !            89:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
        !            90:                error("gcdz : invalid argument");
        !            91:        else if ( NUM(p1) || NUM(p2) ) {
        !            92:                if ( NUM(p1) )
        !            93:                        c1 = (Q)p1;
        !            94:                else
        !            95:                        ptozp(p1,1,&c1,&t);
        !            96:                if ( NUM(p2) )
        !            97:                        c2 = (Q)p2;
        !            98:                else
        !            99:                        ptozp(p2,1,&c2,&t);
        !           100:                gcdn(NM(c1),NM(c2),&n); NTOQ(n,1,c1); *rp = (P)c1;
        !           101:        } else {
        !           102: #if 0
        !           103:                w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
        !           104: #endif
        !           105:                ezgcdpz(CO,p1,p2,rp);
        !           106:        }
        !           107: }
        !           108:
        !           109: void Plcm(arg,rp)
        !           110: NODE arg;
        !           111: P *rp;
        !           112: {
        !           113:        P t1,t2,p1,p2,g,q;
        !           114:        Q c;
        !           115:
        !           116:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
        !           117:        asir_assert(p1,O_P,"lcm");
        !           118:        asir_assert(p2,O_P,"lcm");
        !           119:        if ( !p1 || !p2 )
        !           120:                *rp = 0;
        !           121:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
        !           122:                error("lcm : invalid argument");
        !           123:        else {
        !           124:                ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
        !           125:                ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
        !           126:        }
        !           127: }
        !           128:
        !           129: void Psqfr(arg,rp)
        !           130: NODE arg;
        !           131: LIST *rp;
        !           132: {
        !           133:        DCP dc;
        !           134:
        !           135:        asir_assert(ARG0(arg),O_P,"sqfr");
        !           136:        sqfrp(CO,(P)ARG0(arg),&dc);
        !           137:        dcptolist(dc,rp);
        !           138: }
        !           139:
        !           140: void Pufctrhint(arg,rp)
        !           141: NODE arg;
        !           142: LIST *rp;
        !           143: {
        !           144:        DCP dc;
        !           145:
        !           146:        asir_assert(ARG0(arg),O_P,"ufctrhint");
        !           147:        asir_assert(ARG1(arg),O_N,"ufctrhint");
        !           148:        ufctr((P)ARG0(arg),QTOS((Q)ARG1(arg)),&dc);
        !           149:        dcptolist(dc,rp);
        !           150: }
        !           151:
        !           152: #if 0
        !           153: Pmgcd(arg,rp)
        !           154: NODE arg;
        !           155: Obj *rp;
        !           156: {
        !           157:        NODE node,tn;
        !           158:        int i,m;
        !           159:        P *l;
        !           160:
        !           161:        node = BDY((LIST)ARG0(arg));
        !           162:        for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
        !           163:        m = i; l = (P *)ALLOCA(m*sizeof(P));
        !           164:        for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
        !           165:                l[i] = (P)BDY(tn);
        !           166:        nezgcdnpz(CO,l,m,rp);
        !           167: }
        !           168: #endif
        !           169:
        !           170: void Pcont(arg,rp)
        !           171: NODE arg;
        !           172: P *rp;
        !           173: {
        !           174:        DCP dc;
        !           175:        int m;
        !           176:        P p,p1;
        !           177:        P *l;
        !           178:        V v;
        !           179:
        !           180:        asir_assert(ARG0(arg),O_P,"cont");
        !           181:        p = (P)ARG0(arg);
        !           182:        if ( NUM(p) )
        !           183:                *rp = p;
        !           184:        else {
        !           185:                if ( argc(arg) == 2 ) {
        !           186:                        v = VR((P)ARG1(arg));
        !           187:                        change_mvar(CO,p,v,&p1);
        !           188:                        if ( VR(p1) != v ) {
        !           189:                                *rp = p1; return;
        !           190:                        } else
        !           191:                                p = p1;
        !           192:                }
        !           193:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
        !           194:                l = (P *)ALLOCA(m*sizeof(P));
        !           195:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
        !           196:                        l[m] = COEF(dc);
        !           197:                nezgcdnpz(CO,l,m,rp);
        !           198:        }
        !           199: }
        !           200:
        !           201: void Pptozp(arg,rp)
        !           202: NODE arg;
        !           203: P *rp;
        !           204: {
        !           205:        Q t;
        !           206:
        !           207:        asir_assert(ARG0(arg),O_P,"ptozp");
        !           208:        ptozp((P)ARG0(arg),1,&t,rp);
        !           209: }
        !           210:
        !           211: void Pafctr(arg,rp)
        !           212: NODE arg;
        !           213: LIST *rp;
        !           214: {
        !           215:        DCP dc;
        !           216:
        !           217:        asir_assert(ARG0(arg),O_P,"afctr");
        !           218:        asir_assert(ARG1(arg),O_P,"afctr");
        !           219:        afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
        !           220:        dcptolist(dc,rp);
        !           221: }
        !           222:
        !           223: void Pagcd(arg,rp)
        !           224: NODE arg;
        !           225: P *rp;
        !           226: {
        !           227:        asir_assert(ARG0(arg),O_P,"agcd");
        !           228:        asir_assert(ARG1(arg),O_P,"agcd");
        !           229:        asir_assert(ARG2(arg),O_P,"agcd");
        !           230:        gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
        !           231: }
        !           232:
        !           233: #if 1
        !           234: #define Mulum mulum
        !           235: #define Divum divum
        !           236: #define Mulsum mulsum
        !           237: #define Gcdum gcdum
        !           238: #endif
        !           239:
        !           240: void Mulum(), Mulsum(), Gcdum();
        !           241: int Divum();
        !           242:
        !           243: #define FCTR 0 /* berlekamp */
        !           244: #define SQFR 1
        !           245: #define DDD 2  /* Cantor-Zassenhauss */
        !           246: #define NEWDDD 3  /* berlekamp + root-finding by Cantor-Zassenhauss */
        !           247:
        !           248: UM *resberle();
        !           249:
        !           250: void Pmodfctr(arg,rp)
        !           251: NODE arg;
        !           252: LIST *rp;
        !           253: {
        !           254:        DCP dc;
        !           255:        int mod;
        !           256:
        !           257:        mod = QTOS((Q)ARG1(arg));
        !           258:        if ( mod < 0 )
        !           259:                error("modfctr : invalid modulus");
        !           260:        modfctrp(ARG0(arg),mod,NEWDDD,&dc);
        !           261:        if ( !dc ) {
        !           262:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
        !           263:        }
        !           264:        dcptolist(dc,rp);
        !           265: }
        !           266:
        !           267: void Pmodsqfr(arg,rp)
        !           268: NODE arg;
        !           269: LIST *rp;
        !           270: {
        !           271:        DCP dc;
        !           272:
        !           273:        if ( !dc ) {
        !           274:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
        !           275:        }
        !           276:        modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc);
        !           277:        dcptolist(dc,rp);
        !           278: }
        !           279:
        !           280: void Pddd(arg,rp)
        !           281: NODE arg;
        !           282: LIST *rp;
        !           283: {
        !           284:        DCP dc;
        !           285:
        !           286:        if ( !dc ) {
        !           287:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
        !           288:        }
        !           289:        modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc);
        !           290:        dcptolist(dc,rp);
        !           291: }
        !           292:
        !           293: void Pnewddd(arg,rp)
        !           294: NODE arg;
        !           295: LIST *rp;
        !           296: {
        !           297:        DCP dc;
        !           298:
        !           299:        if ( !dc ) {
        !           300:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
        !           301:        }
        !           302:        modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc);
        !           303:        dcptolist(dc,rp);
        !           304: }
        !           305:
        !           306: void Pirred_check(arg,rp)
        !           307: NODE arg;
        !           308: Q *rp;
        !           309: {
        !           310:        P p;
        !           311:        UM mp;
        !           312:        int r,mod;
        !           313:
        !           314:        p = (P)ARG0(arg);
        !           315:        if ( !p ) {
        !           316:                *rp = 0; return;
        !           317:        }
        !           318:        mp = W_UMALLOC(UDEG(p));
        !           319:        mod = QTOS((Q)ARG1(arg));
        !           320:        ptoum(mod,p,mp);
        !           321:        r = irred_check(mp,mod);
        !           322:        if ( r )
        !           323:                *rp = ONE;
        !           324:        else
        !           325:                *rp = 0;
        !           326: }
        !           327:
        !           328: void Pnfctr_mod(arg,rp)
        !           329: NODE arg;
        !           330: Q *rp;
        !           331: {
        !           332:        P p;
        !           333:        UM mp;
        !           334:        int r,mod;
        !           335:
        !           336:        p = (P)ARG0(arg);
        !           337:        if ( !p ) {
        !           338:                *rp = 0; return;
        !           339:        }
        !           340:        mp = W_UMALLOC(UDEG(p));
        !           341:        mod = QTOS((Q)ARG1(arg));
        !           342:        ptoum(mod,p,mp);
        !           343:        r = nfctr_mod(mp,mod);
        !           344:        STOQ(r,*rp);
        !           345: }
        !           346:
        !           347: void Pddd_tab(arg,rp)
        !           348: NODE arg;
        !           349: VECT *rp;
        !           350: {
        !           351:        P p;
        !           352:        UM mp,t,q,r1,w,w1;
        !           353:        UM *r,*s;
        !           354:        int dr,mod,n,i;
        !           355:        VECT result;
        !           356:        V v;
        !           357:
        !           358:        p = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
        !           359:        v = VR(p);
        !           360:        n = UDEG(p); mp = W_UMALLOC(n);
        !           361:        ptoum(mod,p,mp);
        !           362:        r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
        !           363:        r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
        !           364:        t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
        !           365:        DEG(t) = mod; COEF(t)[mod] = 1;
        !           366:        q = W_UMALLOC(mod);
        !           367:        dr = divum(mod,t,mp,q);
        !           368:        DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
        !           369:        s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
        !           370:        w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
        !           371:        w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
        !           372:        for ( i = 1; i < n; i++ ) {
        !           373:                DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
        !           374:                mulum(mod,r1,w,w1);
        !           375:                dr = divum(mod,w1,mp,q); DEG(w1) = dr;
        !           376:                s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
        !           377:        }
        !           378:        for ( i = 2; i < n; i++ ) {
        !           379:                mult_mod_tab(r[i-1],mod,s,w,n);
        !           380:                r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
        !           381:        }
        !           382:        MKVECT(result,n);
        !           383:        for ( i = 0; i < n; i++ )
        !           384:                umtop(v,r[i],(P *)&BDY(result)[i]);
        !           385:        *rp = result;
        !           386: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>