[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.19

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.19    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/fctr.c,v 1.18 2003/01/06 01:16:38 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "parse.h"
                     52:
                     53: void Pfctr(), Pgcd(), Pgcdz(), Plcm(), Psqfr(), Pufctrhint();
1.17      noro       54: void Pptozp(), Pcont(), Psfcont();
1.1       noro       55: void Pafctr(), Pagcd();
                     56: void Pmodsqfr(),Pmodfctr(),Pddd(),Pnewddd(),Pddd_tab();
1.15      noro       57: void Psfsqfr(),Psffctr(),Psfbfctr(),Psfufctr(),Psfmintdeg(),Psfgcd();
1.1       noro       58: void Pirred_check(), Pnfctr_mod();
1.16      noro       59: void Pbivariate_hensel_special();
1.1       noro       60:
1.11      noro       61: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr);
                     62:
1.1       noro       63: struct ftab fctr_tab[] = {
1.16      noro       64:        {"bivariate_hensel_special",Pbivariate_hensel_special,6},
1.5       noro       65:        {"fctr",Pfctr,-2},
1.1       noro       66:        {"gcd",Pgcd,-3},
                     67:        {"gcdz",Pgcdz,2},
                     68:        {"lcm",Plcm,2},
                     69:        {"sqfr",Psqfr,1},
                     70:        {"ufctrhint",Pufctrhint,2},
                     71:        {"ptozp",Pptozp,1},
                     72:        {"cont",Pcont,-2},
1.17      noro       73:        {"sfcont",Psfcont,-2},
1.1       noro       74:        {"afctr",Pafctr,2},
                     75:        {"agcd",Pagcd,3},
                     76:        {"modsqfr",Pmodsqfr,2},
                     77:        {"modfctr",Pmodfctr,2},
1.10      noro       78:        {"sfsqfr",Psfsqfr,1},
1.15      noro       79:        {"sffctr",Psffctr,1},
1.10      noro       80:        {"sfufctr",Psfufctr,1},
                     81:        {"sfbfctr",Psfbfctr,-4},
1.11      noro       82:        {"sfmintdeg",Psfmintdeg,5},
1.13      noro       83:        {"sfgcd",Psfgcd,2},
1.1       noro       84: #if 0
                     85:        {"ddd",Pddd,2},
                     86:        {"newddd",Pnewddd,2},
                     87: #endif
                     88:        {"ddd_tab",Pddd_tab,2},
                     89:        {"irred_check",Pirred_check,2},
                     90:        {"nfctr_mod",Pnfctr_mod,2},
                     91:        {0,0,0},
                     92: };
1.16      noro       93:
                     94: /* bivariate_hensel_special(f(x,y):monic in x,g0(x),h0(y),x,y,d) */
                     95:
                     96: void Pbivariate_hensel_special(arg,rp)
                     97: NODE arg;
                     98: LIST *rp;
                     99: {
                    100:        DCP dc;
                    101:        struct oVN vn[2];
                    102:        P f,g0,h0,ak,bk,gk,hk;
                    103:        V vx,vy;
                    104:        VL nvl;
                    105:        Q qk,cbd,bb;
                    106:        int d;
                    107:        NODE n;
                    108:
                    109:        f = (P)ARG0(arg);
                    110:        g0 = (P)ARG1(arg);
                    111:        h0 = (P)ARG2(arg);
                    112:        vx = VR((P)ARG3(arg));
                    113:        vy = VR((P)ARG4(arg));
                    114:        d = QTOS((Q)ARG5(arg));
                    115:        NEWVL(nvl); nvl->v = vx;
                    116:        NEWVL(NEXT(nvl)); NEXT(nvl)->v = vy;
                    117:        NEXT(NEXT(nvl)) = 0;
                    118:        vn[0].v = vy; vn[0].n = 0;
                    119:        vn[1].v = 0; vn[1].n = 0;
                    120:        cbound(nvl,f,&cbd);
                    121:        addq(cbd,cbd,&bb);
                    122:        henzq1(g0,h0,bb,&bk,&ak,&qk);
                    123:        henmv(nvl,vn,f,g0,h0,ak,bk,(P)ONE,(P)ONE,(P)ONE,(P)ONE,qk,d,&gk,&hk);
                    124:        n = mknode(2,gk,hk);
                    125:        MKLIST(*rp,n);
                    126: }
1.1       noro      127:
                    128: void Pfctr(arg,rp)
                    129: NODE arg;
                    130: LIST *rp;
                    131: {
                    132:        DCP dc;
                    133:
                    134:        asir_assert(ARG0(arg),O_P,"fctr");
1.5       noro      135:        if ( argc(arg) == 1 )
                    136:                fctrp(CO,(P)ARG0(arg),&dc);
                    137:        else {
                    138:                asir_assert(ARG1(arg),O_P,"fctr");
                    139:                fctr_wrt_v_p(CO,(P)ARG0(arg),VR((P)ARG1(arg)),&dc);
                    140:        }
1.1       noro      141:        dcptolist(dc,rp);
                    142: }
                    143:
                    144: void Pgcd(arg,rp)
                    145: NODE arg;
                    146: P *rp;
                    147: {
                    148:        P p1,p2,g1,g2,g;
                    149:        Num m;
                    150:        int mod;
                    151:
                    152:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
                    153:        asir_assert(p1,O_P,"gcd");
                    154:        asir_assert(p2,O_P,"gcd");
                    155:        if ( !p1 )
                    156:                *rp = p2;
                    157:        else if ( !p2 )
                    158:                *rp = p1;
                    159:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
1.4       noro      160:                gcdprsp(CO,p1,p2,rp);
1.1       noro      161:        else if ( argc(arg) == 2 )
                    162:                ezgcdp(CO,p1,p2,rp);
                    163:        else {
                    164:                m = (Num)ARG2(arg);
                    165:                asir_assert(m,O_P,"gcd");
                    166:                mod = QTOS((Q)m);
                    167:                ptomp(mod,p1,&g1); ptomp(mod,p2,&g2);
                    168:                gcdprsmp(CO,mod,g1,g2,&g);
                    169:                mptop(g,rp);
                    170:        }
                    171: }
                    172:
                    173: void Pgcdz(arg,rp)
                    174: NODE arg;
                    175: P *rp;
                    176: {
                    177:        P p1,p2,t;
                    178:        Q c1,c2;
                    179:        N n;
                    180:
                    181:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
                    182:        asir_assert(p1,O_P,"gcdz");
                    183:        asir_assert(p2,O_P,"gcdz");
                    184:        if ( !p1 )
                    185:                *rp = p2;
                    186:        else if ( !p2 )
                    187:                *rp = p1;
                    188:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
                    189:                error("gcdz : invalid argument");
                    190:        else if ( NUM(p1) || NUM(p2) ) {
                    191:                if ( NUM(p1) )
                    192:                        c1 = (Q)p1;
                    193:                else
                    194:                        ptozp(p1,1,&c1,&t);
                    195:                if ( NUM(p2) )
                    196:                        c2 = (Q)p2;
                    197:                else
                    198:                        ptozp(p2,1,&c2,&t);
                    199:                gcdn(NM(c1),NM(c2),&n); NTOQ(n,1,c1); *rp = (P)c1;
                    200:        } else {
                    201: #if 0
                    202:                w[0] = p1; w[1] = p2; nezgcdnpz(CO,w,2,rp);
                    203: #endif
                    204:                ezgcdpz(CO,p1,p2,rp);
                    205:        }
                    206: }
                    207:
                    208: void Plcm(arg,rp)
                    209: NODE arg;
                    210: P *rp;
                    211: {
                    212:        P t1,t2,p1,p2,g,q;
                    213:        Q c;
                    214:
                    215:        p1 = (P)ARG0(arg); p2 = (P)ARG1(arg);
                    216:        asir_assert(p1,O_P,"lcm");
                    217:        asir_assert(p2,O_P,"lcm");
                    218:        if ( !p1 || !p2 )
                    219:                *rp = 0;
                    220:        else if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) )
                    221:                error("lcm : invalid argument");
                    222:        else {
                    223:                ptozp(p1,1,&c,&t1); ptozp(p2,1,&c,&t2);
                    224:                ezgcdp(CO,t1,t2,&g); divsp(CO,t1,g,&q); mulp(CO,q,t2,rp);
                    225:        }
                    226: }
                    227:
                    228: void Psqfr(arg,rp)
                    229: NODE arg;
                    230: LIST *rp;
                    231: {
                    232:        DCP dc;
                    233:
                    234:        asir_assert(ARG0(arg),O_P,"sqfr");
                    235:        sqfrp(CO,(P)ARG0(arg),&dc);
                    236:        dcptolist(dc,rp);
                    237: }
                    238:
                    239: void Pufctrhint(arg,rp)
                    240: NODE arg;
                    241: LIST *rp;
                    242: {
                    243:        DCP dc;
                    244:
                    245:        asir_assert(ARG0(arg),O_P,"ufctrhint");
                    246:        asir_assert(ARG1(arg),O_N,"ufctrhint");
                    247:        ufctr((P)ARG0(arg),QTOS((Q)ARG1(arg)),&dc);
                    248:        dcptolist(dc,rp);
                    249: }
                    250:
                    251: #if 0
                    252: Pmgcd(arg,rp)
                    253: NODE arg;
                    254: Obj *rp;
                    255: {
                    256:        NODE node,tn;
                    257:        int i,m;
                    258:        P *l;
                    259:
                    260:        node = BDY((LIST)ARG0(arg));
                    261:        for ( i = 0, tn = node; tn; tn = NEXT(tn), i++ );
                    262:        m = i; l = (P *)ALLOCA(m*sizeof(P));
                    263:        for ( i = 0, tn = node; i < m; tn = NEXT(tn), i++ )
                    264:                l[i] = (P)BDY(tn);
                    265:        nezgcdnpz(CO,l,m,rp);
                    266: }
                    267: #endif
                    268:
                    269: void Pcont(arg,rp)
                    270: NODE arg;
                    271: P *rp;
                    272: {
                    273:        DCP dc;
                    274:        int m;
                    275:        P p,p1;
                    276:        P *l;
                    277:        V v;
                    278:
                    279:        asir_assert(ARG0(arg),O_P,"cont");
                    280:        p = (P)ARG0(arg);
                    281:        if ( NUM(p) )
                    282:                *rp = p;
                    283:        else {
                    284:                if ( argc(arg) == 2 ) {
                    285:                        v = VR((P)ARG1(arg));
                    286:                        change_mvar(CO,p,v,&p1);
                    287:                        if ( VR(p1) != v ) {
                    288:                                *rp = p1; return;
                    289:                        } else
                    290:                                p = p1;
                    291:                }
                    292:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
                    293:                l = (P *)ALLOCA(m*sizeof(P));
                    294:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
                    295:                        l[m] = COEF(dc);
                    296:                nezgcdnpz(CO,l,m,rp);
1.17      noro      297:        }
                    298: }
                    299:
                    300: void Psfcont(arg,rp)
                    301: NODE arg;
                    302: P *rp;
                    303: {
                    304:        DCP dc;
1.18      noro      305:        MP mp;
1.17      noro      306:        int m;
1.18      noro      307:        Obj obj;
1.17      noro      308:        P p,p1;
                    309:        P *l;
                    310:        V v;
                    311:
1.18      noro      312:        obj = (Obj)ARG0(arg);
                    313:        if ( !obj || NUM(obj) )
                    314:                *rp = (P)obj;
                    315:        else if ( OID(obj) == O_P ) {
                    316:                p = (P)obj;
1.17      noro      317:                if ( argc(arg) == 2 ) {
                    318:                        v = VR((P)ARG1(arg));
                    319:                        change_mvar(CO,p,v,&p1);
                    320:                        if ( VR(p1) != v ) {
                    321:                                *rp = p1; return;
                    322:                        } else
                    323:                                p = p1;
                    324:                }
                    325:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ );
                    326:                l = (P *)ALLOCA(m*sizeof(P));
                    327:                for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ )
                    328:                        l[m] = COEF(dc);
                    329:                gcdsf(CO,l,m,rp);
1.18      noro      330:        } else if ( OID(obj) == O_DP ) {
                    331:                for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++ );
                    332:                l = (P *)ALLOCA(m*sizeof(P));
                    333:                for ( m = 0, mp = BDY((DP)obj); mp; mp = NEXT(mp), m++)
                    334:                        l[m] = mp->c;
                    335:                gcdsf(CO,l,m,rp);
1.1       noro      336:        }
                    337: }
                    338:
                    339: void Pptozp(arg,rp)
                    340: NODE arg;
                    341: P *rp;
                    342: {
                    343:        Q t;
                    344:
                    345:        asir_assert(ARG0(arg),O_P,"ptozp");
                    346:        ptozp((P)ARG0(arg),1,&t,rp);
                    347: }
                    348:
                    349: void Pafctr(arg,rp)
                    350: NODE arg;
                    351: LIST *rp;
                    352: {
                    353:        DCP dc;
                    354:
                    355:        asir_assert(ARG0(arg),O_P,"afctr");
                    356:        asir_assert(ARG1(arg),O_P,"afctr");
                    357:        afctr(CO,(P)ARG0(arg),(P)ARG1(arg),&dc);
                    358:        dcptolist(dc,rp);
                    359: }
                    360:
                    361: void Pagcd(arg,rp)
                    362: NODE arg;
                    363: P *rp;
                    364: {
                    365:        asir_assert(ARG0(arg),O_P,"agcd");
                    366:        asir_assert(ARG1(arg),O_P,"agcd");
                    367:        asir_assert(ARG2(arg),O_P,"agcd");
                    368:        gcda(CO,(P)ARG0(arg),(P)ARG1(arg),(P)ARG2(arg),rp);
                    369: }
                    370:
                    371: #if 1
                    372: #define Mulum mulum
                    373: #define Divum divum
                    374: #define Mulsum mulsum
                    375: #define Gcdum gcdum
                    376: #endif
                    377:
                    378: void Mulum(), Mulsum(), Gcdum();
                    379: int Divum();
                    380:
                    381: #define FCTR 0 /* berlekamp */
                    382: #define SQFR 1
                    383: #define DDD 2  /* Cantor-Zassenhauss */
                    384: #define NEWDDD 3  /* berlekamp + root-finding by Cantor-Zassenhauss */
                    385:
                    386: UM *resberle();
                    387:
1.18      noro      388: void reduce_sfdc(DCP sfdc, DCP *dc);
                    389:
1.1       noro      390: void Pmodfctr(arg,rp)
                    391: NODE arg;
                    392: LIST *rp;
                    393: {
1.18      noro      394:        DCP dc,dcu;
                    395:        int mod,i,t;
                    396:        P p;
                    397:        Obj u;
                    398:        VL vl;
1.1       noro      399:
                    400:        mod = QTOS((Q)ARG1(arg));
                    401:        if ( mod < 0 )
                    402:                error("modfctr : invalid modulus");
1.18      noro      403:        p = (P)ARG0(arg);
                    404:        clctv(CO,p,&vl);
1.19    ! noro      405:        if ( !vl ) {
        !           406:                NEWDC(dc); COEF(dc) = p; DEG(dc) = ONE; NEXT(dc) = 0;
        !           407:        } else if ( !NEXT(vl) )
1.18      noro      408:                modfctrp(ARG0(arg),mod,NEWDDD,&dc);
                    409:        else {
                    410:                /* XXX 16384 should be replaced by a macro */
                    411:                for ( i = 0, t = 1; t*mod < 16384; t *= mod, i++ );
                    412:                current_ff = FF_GFS;
                    413:                setmod_sf(mod,i);
                    414:                simp_ff((Obj)p,&u);
                    415:                mfctrsf(CO,(P)u,&dcu);
                    416:                reduce_sfdc(dcu,&dc);
                    417:        }
1.6       noro      418:        if ( !dc ) {
                    419:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
                    420:        }
                    421:        dcptolist(dc,rp);
1.13      noro      422: }
                    423:
                    424: void Psfgcd(arg,rp)
                    425: NODE arg;
                    426: LIST *rp;
                    427: {
                    428:        P ps[2];
                    429:
                    430:        ps[0] = (P)ARG0(arg);
                    431:        ps[1] = (P)ARG1(arg);
                    432:        gcdsf(CO,ps,2,rp);
1.6       noro      433: }
                    434:
1.15      noro      435: void Psffctr(arg,rp)
                    436: NODE arg;
                    437: LIST *rp;
                    438: {
                    439:        DCP dc;
                    440:
                    441:        mfctrsf(CO,ARG0(arg),&dc);
                    442:        dcptolist(dc,rp);
                    443: }
                    444:
1.10      noro      445: void Psfsqfr(arg,rp)
                    446: NODE arg;
                    447: LIST *rp;
                    448: {
                    449:        DCP dc;
                    450:
1.14      noro      451:        sqfrsf(CO,ARG0(arg),&dc);
1.10      noro      452:        dcptolist(dc,rp);
                    453: }
                    454:
                    455: void Psfufctr(arg,rp)
1.6       noro      456: NODE arg;
                    457: LIST *rp;
                    458: {
                    459:        DCP dc;
                    460:
1.15      noro      461:        ufctrsf(ARG0(arg),&dc);
1.1       noro      462:        dcptolist(dc,rp);
1.7       noro      463: }
                    464:
                    465: void Psfbfctr(arg,rp)
                    466: NODE arg;
                    467: LIST *rp;
                    468: {
                    469:        V x,y;
1.8       noro      470:        DCP dc,dct;
                    471:        P t;
                    472:        struct oVL vl1,vl2;
                    473:        VL vl;
1.10      noro      474:        int degbound;
1.7       noro      475:
                    476:        x = VR((P)ARG1(arg));
                    477:        y = VR((P)ARG2(arg));
1.8       noro      478:        vl1.v = x; vl1.next = &vl2;
                    479:        vl2.v = y; vl2.next = 0;
                    480:        vl = &vl1;
1.10      noro      481:        if ( argc(arg) == 4 )
                    482:                degbound = QTOS((Q)ARG3(arg));
                    483:        else
                    484:                degbound = -1;
1.8       noro      485:
1.10      noro      486:        sfbfctr((P)ARG0(arg),x,y,degbound,&dc);
1.8       noro      487:        for ( dct = dc; dct; dct = NEXT(dct) ) {
                    488:                reorderp(CO,vl,COEF(dct),&t); COEF(dct) = t;
1.7       noro      489:        }
1.8       noro      490:        dcptolist(dc,rp);
1.1       noro      491: }
                    492:
1.11      noro      493: void Psfmintdeg(arg,rp)
                    494: NODE arg;
                    495: P *rp;
                    496: {
                    497:        V x,y;
                    498:        P r;
                    499:        struct oVL vl1,vl2;
                    500:        VL vl;
                    501:        int dy,c;
                    502:
                    503:        x = VR((P)ARG1(arg));
                    504:        y = VR((P)ARG2(arg));
                    505:        vl1.v = x; vl1.next = &vl2;
                    506:        vl2.v = y; vl2.next = 0;
                    507:        vl = &vl1;
                    508:        dy = QTOS((Q)ARG3(arg));
                    509:        c = QTOS((Q)ARG4(arg));
                    510:        sfmintdeg(vl,(P)ARG0(arg),dy,c,&r);
                    511:        reorderp(CO,vl,r,rp);
                    512: }
                    513:
1.1       noro      514: void Pmodsqfr(arg,rp)
                    515: NODE arg;
                    516: LIST *rp;
                    517: {
                    518:        DCP dc;
                    519:
1.9       noro      520:        if ( !ARG0(arg) ) {
1.1       noro      521:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9       noro      522:        } else
                    523:                modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),SQFR,&dc);
1.1       noro      524:        dcptolist(dc,rp);
                    525: }
                    526:
                    527: void Pddd(arg,rp)
                    528: NODE arg;
                    529: LIST *rp;
                    530: {
                    531:        DCP dc;
                    532:
1.9       noro      533:        if ( !ARG0(arg) ) {
1.1       noro      534:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9       noro      535:        } else
                    536:                modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),DDD,&dc);
1.1       noro      537:        dcptolist(dc,rp);
                    538: }
                    539:
                    540: void Pnewddd(arg,rp)
                    541: NODE arg;
                    542: LIST *rp;
                    543: {
1.9       noro      544:        DCP dc=0;
1.1       noro      545:
1.9       noro      546:        if ( !ARG0(arg) ) {
1.1       noro      547:                NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; NEXT(dc) = 0;
1.9       noro      548:        } else
                    549:                modfctrp(ARG0(arg),QTOS((Q)ARG1(arg)),NEWDDD,&dc);
1.1       noro      550:        dcptolist(dc,rp);
                    551: }
                    552:
                    553: void Pirred_check(arg,rp)
                    554: NODE arg;
                    555: Q *rp;
                    556: {
                    557:        P p;
                    558:        UM mp;
                    559:        int r,mod;
                    560:
                    561:        p = (P)ARG0(arg);
                    562:        if ( !p ) {
                    563:                *rp = 0; return;
                    564:        }
                    565:        mp = W_UMALLOC(UDEG(p));
                    566:        mod = QTOS((Q)ARG1(arg));
                    567:        ptoum(mod,p,mp);
                    568:        r = irred_check(mp,mod);
                    569:        if ( r )
                    570:                *rp = ONE;
                    571:        else
                    572:                *rp = 0;
                    573: }
                    574:
                    575: void Pnfctr_mod(arg,rp)
                    576: NODE arg;
                    577: Q *rp;
                    578: {
                    579:        P p;
                    580:        UM mp;
                    581:        int r,mod;
                    582:
                    583:        p = (P)ARG0(arg);
                    584:        if ( !p ) {
                    585:                *rp = 0; return;
                    586:        }
                    587:        mp = W_UMALLOC(UDEG(p));
                    588:        mod = QTOS((Q)ARG1(arg));
                    589:        ptoum(mod,p,mp);
                    590:        r = nfctr_mod(mp,mod);
                    591:        STOQ(r,*rp);
                    592: }
                    593:
                    594: void Pddd_tab(arg,rp)
                    595: NODE arg;
                    596: VECT *rp;
                    597: {
                    598:        P p;
                    599:        UM mp,t,q,r1,w,w1;
                    600:        UM *r,*s;
                    601:        int dr,mod,n,i;
                    602:        VECT result;
                    603:        V v;
                    604:
                    605:        p = (P)ARG0(arg); mod = QTOS((Q)ARG1(arg));
                    606:        v = VR(p);
                    607:        n = UDEG(p); mp = W_UMALLOC(n);
                    608:        ptoum(mod,p,mp);
                    609:        r = (UM *)W_ALLOC(n); s = (UM *)W_ALLOC(n);
                    610:        r[0] = UMALLOC(0); DEG(r[0]) = 0; COEF(r[0])[0] = 1;
                    611:        t = W_UMALLOC(mod); bzero(COEF(t),sizeof(int)*(mod+1));
                    612:        DEG(t) = mod; COEF(t)[mod] = 1;
                    613:        q = W_UMALLOC(mod);
                    614:        dr = divum(mod,t,mp,q);
                    615:        DEG(t) = dr; r[1] = r1 = UMALLOC(dr); cpyum(t,r1);
                    616:        s[0] = W_UMALLOC(dr); cpyum(t,s[0]);
                    617:        w = W_UMALLOC(n); bzero(COEF(w),sizeof(int)*(n+1));
                    618:        w1 = W_UMALLOC(2*n); bzero(COEF(w1),sizeof(int)*(2*n+1));
                    619:        for ( i = 1; i < n; i++ ) {
                    620:                DEG(w) = i; COEF(w)[i-1] = 0; COEF(w)[i] = 1;
                    621:                mulum(mod,r1,w,w1);
                    622:                dr = divum(mod,w1,mp,q); DEG(w1) = dr;
                    623:                s[i] = W_UMALLOC(dr); cpyum(w1,s[i]);
                    624:        }
                    625:        for ( i = 2; i < n; i++ ) {
                    626:                mult_mod_tab(r[i-1],mod,s,w,n);
                    627:                r[i] = UMALLOC(DEG(w)); cpyum(w,r[i]);
                    628:        }
                    629:        MKVECT(result,n);
                    630:        for ( i = 0; i < n; i++ )
                    631:                umtop(v,r[i],(P *)&BDY(result)[i]);
                    632:        *rp = result;
1.18      noro      633: }
                    634:
                    635: void reduce_sfdc(DCP sfdc,DCP *dcr)
                    636: {
                    637:        P c,t,s,u,f;
                    638:        DCP dc0,dc,tdc;
                    639:        DCP *a;
                    640:        int i,j,n;
                    641:
                    642:        if ( !current_gfs_ext ) {
                    643:                /* we simply apply sfptop() */
                    644:                for ( dc0 = 0; sfdc; sfdc = NEXT(sfdc) ) {
                    645:                        NEXTDC(dc0,dc);
                    646:                        DEG(dc) = DEG(sfdc);
                    647:                        sfptop(COEF(sfdc),&COEF(dc));
                    648:                }
                    649:                NEXT(dc) = 0;
                    650:                *dcr = dc0;
                    651:                return;
                    652:        }
                    653:
                    654:        if ( NUM(COEF(sfdc)) ) {
                    655:                sfptop(COEF(sfdc),&c);
                    656:                sfdc = NEXT(sfdc);
                    657:        } else
                    658:                c = (P)ONE;
                    659:
                    660:        for ( n = 0, tdc = sfdc; tdc; tdc = NEXT(tdc), n++ );
                    661:        a = (DCP *)ALLOCA(n*sizeof(DCP));
                    662:        for ( i = 0, tdc = sfdc; i < n; tdc = NEXT(tdc), i++ )
                    663:                a[i] = tdc;
                    664:
                    665:        dc0 = 0; NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = c;
                    666:        for ( i = 0; i < n; i++ ) {
                    667:                if ( !a[i] )
                    668:                        continue;
                    669:                t = COEF(a[i]);
                    670:                f = t;
                    671:                while ( 1 ) {
                    672:                        sf_galois_action(t,ONE,&s);
                    673:                        for ( j = i; j < n; j++ )
                    674:                                if ( a[j] && !compp(CO,s,COEF(a[j])) )
                    675:                                        break;
                    676:                        if ( j == n )
                    677:                                error("reduce_sfdc : cannot happen");
                    678:                        if ( j == i ) {
                    679:                                NEXTDC(dc0,dc); DEG(dc) = DEG(a[i]);
                    680:                                sfptop(f,&COEF(dc));
                    681:                                break;
                    682:                        } else {
                    683:                                mulp(CO,f,s,&u); f = u;
                    684:                                t = s;
                    685:                                a[j] = 0;
                    686:                        }
                    687:                }
                    688:        }
                    689:        *dcr = dc0;
1.1       noro      690: }

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