[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

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>