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

Annotation of OpenXM_contrib2/asir2000/engine/distm.c, Revision 1.6

1.5       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.6     ! noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.5       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.6     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.5 2000/08/21 08:31:27 noro Exp $
1.5       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "inline.h"
                     52:
                     53: #define NV(p) ((p)->nv)
                     54: #define C(p) ((p)->c)
                     55: #if 0
                     56: #define ITOS(p) (((unsigned int)(p))>>1)
                     57: #define STOI(i) ((P)((((unsigned int)(i))<<1)|1))
                     58: #else
                     59: #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
                     60: #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
                     61: #endif
                     62:
                     63: extern int (*cmpdl)();
1.2       noro       64: extern int do_weyl;
                     65:
                     66: void comm_mulmd();
                     67: void weyl_mulmd();
                     68: void weyl_mulmdm();
                     69: void weyl_mulmmm();
                     70: void _comm_mulmd();
                     71: void _weyl_mulmd();
                     72: void _weyl_mulmmm();
                     73: void _weyl_mulmdm();
1.1       noro       74:
                     75: void ptomd(vl,mod,dvl,p,pr)
                     76: VL vl,dvl;
                     77: int mod;
                     78: P p;
                     79: DP *pr;
                     80: {
                     81:        P t;
                     82:
                     83:        ptomp(mod,p,&t);
                     84:        mptomd(vl,mod,dvl,t,pr);
                     85: }
                     86:
                     87: void mptomd(vl,mod,dvl,p,pr)
                     88: VL vl,dvl;
                     89: int mod;
                     90: P p;
                     91: DP *pr;
                     92: {
                     93:        int n,i;
                     94:        VL tvl;
                     95:        V v;
                     96:        DL d;
                     97:        MP m;
                     98:        DCP dc;
                     99:        DP r,s,t,u;
                    100:        P x,c;
                    101:
                    102:        if ( !p )
                    103:                *pr = 0;
                    104:        else {
                    105:                for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
                    106:                if ( NUM(p) ) {
                    107:                        NEWDL(d,n);
                    108:                        NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr);
                    109:                } else {
                    110:                        for ( i = 0, tvl = dvl, v = VR(p);
                    111:                                tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
                    112:                        if ( !tvl ) {
                    113:                                for ( dc = DC(p), s = 0, MKV(v,x); dc; dc = NEXT(dc) ) {
                    114:                                        mptomd(vl,mod,dvl,COEF(dc),&t); pwrmp(vl,mod,x,DEG(dc),&c);
                    115:                                        mulmdc(vl,mod,t,c,&r); addmd(vl,mod,r,s,&t); s = t;
                    116:                                }
                    117:                                *pr = s;
                    118:                        } else {
                    119:                                for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
                    120:                                        mptomd(vl,mod,dvl,COEF(dc),&t);
                    121:                                        NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td;
                    122:                                        NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u);
1.2       noro      123:                                        comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t;
1.1       noro      124:                                }
                    125:                                *pr = s;
                    126:                        }
                    127:                }
                    128:        }
                    129: }
                    130:
                    131: void mdtop(vl,mod,dvl,p,pr)
                    132: VL vl,dvl;
                    133: int mod;
                    134: DP p;
                    135: P *pr;
                    136: {
                    137:        int n,i;
                    138:        DL d;
                    139:        MP m;
                    140:        P r,s,t,u,w;
                    141:        Q q;
                    142:        VL tvl;
                    143:
                    144:        if ( !p )
                    145:                *pr = 0;
                    146:        else {
                    147:                for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
                    148:                        for ( i = 0, t = C(m), d = m->dl, tvl = dvl;
                    149:                                i < n; tvl = NEXT(tvl), i++ ) {
                    150:                                MKMV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
                    151:                                mulmp(vl,mod,t,u,&w); t = w;
                    152:                        }
                    153:                        addmp(vl,mod,s,t,&u); s = u;
                    154:                }
                    155:                mptop(s,pr);
                    156:        }
                    157: }
                    158:
                    159: void addmd(vl,mod,p1,p2,pr)
                    160: VL vl;
                    161: int mod;
                    162: DP p1,p2,*pr;
                    163: {
                    164:        int n;
                    165:        MP m1,m2,mr,mr0;
                    166:        P t;
                    167:        int tmp;
                    168:        MQ q;
                    169:
                    170:        if ( !p1 )
                    171:                *pr = p2;
                    172:        else if ( !p2 )
                    173:                *pr = p1;
                    174:        else {
                    175:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    176:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    177:                                case 0:
                    178:                                        if ( NUM(C(m1)) && NUM(C(m2)) ) {
                    179:                                                tmp = (CONT((MQ)C(m1))+CONT((MQ)C(m2))) - mod;
                    180:                                                if ( tmp < 0 )
                    181:                                                        tmp += mod;
                    182:                                                if ( tmp ) {
                    183:                                                        STOMQ(tmp,q); t = (P)q;
                    184:                                                } else
                    185:                                                        t = 0;
                    186:                                        } else
                    187:                                                addmp(vl,mod,C(m1),C(m2),&t);
                    188:                                        if ( t ) {
                    189:                                                NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
                    190:                                        }
                    191:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    192:                                case 1:
                    193:                                        NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    194:                                        m1 = NEXT(m1); break;
                    195:                                case -1:
                    196:                                        NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    197:                                        m2 = NEXT(m2); break;
                    198:                        }
                    199:                if ( !mr0 )
                    200:                        if ( m1 )
                    201:                                mr0 = m1;
                    202:                        else if ( m2 )
                    203:                                mr0 = m2;
                    204:                        else {
                    205:                                *pr = 0;
                    206:                                return;
                    207:                        }
                    208:                else if ( m1 )
                    209:                        NEXT(mr) = m1;
                    210:                else if ( m2 )
                    211:                        NEXT(mr) = m2;
                    212:                else
                    213:                        NEXT(mr) = 0;
                    214:                MKDP(NV(p1),mr0,*pr);
                    215:                if ( *pr )
                    216:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    217:        }
                    218: }
                    219:
                    220: void submd(vl,mod,p1,p2,pr)
                    221: VL vl;
                    222: int mod;
                    223: DP p1,p2,*pr;
                    224: {
                    225:        DP t;
                    226:
                    227:        if ( !p2 )
                    228:                *pr = p1;
                    229:        else {
                    230:                chsgnmd(mod,p2,&t); addmd(vl,mod,p1,t,pr);
                    231:        }
                    232: }
                    233:
                    234: void chsgnmd(mod,p,pr)
                    235: int mod;
                    236: DP p,*pr;
                    237: {
                    238:        MP m,mr,mr0;
                    239:
                    240:        if ( !p )
                    241:                *pr = 0;
                    242:        else {
                    243:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    244:                        NEXTMP(mr0,mr); chsgnmp(mod,C(m),&C(mr)); mr->dl = m->dl;
                    245:                }
                    246:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    247:                if ( *pr )
                    248:                        (*pr)->sugar = p->sugar;
                    249:        }
                    250: }
                    251:
                    252: void mulmd(vl,mod,p1,p2,pr)
                    253: VL vl;
                    254: int mod;
                    255: DP p1,p2,*pr;
                    256: {
1.2       noro      257:        if ( !do_weyl )
                    258:                comm_mulmd(vl,mod,p1,p2,pr);
                    259:        else
                    260:                weyl_mulmd(vl,mod,p1,p2,pr);
                    261: }
                    262:
                    263: void comm_mulmd(vl,mod,p1,p2,pr)
                    264: VL vl;
                    265: int mod;
                    266: DP p1,p2,*pr;
                    267: {
                    268:        MP m;
                    269:        DP s,t,u;
                    270:        int i,l,l1;
                    271:        static MP *w;
                    272:        static int wlen;
                    273:
                    274:        if ( !p1 || !p2 )
                    275:                *pr = 0;
                    276:        else if ( OID(p1) <= O_P )
                    277:                mulmdc(vl,mod,p2,(P)p1,pr);
                    278:        else if ( OID(p2) <= O_P )
                    279:                mulmdc(vl,mod,p1,(P)p2,pr);
                    280:        else {
                    281:                for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    282:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    283:                if ( l1 < l ) {
                    284:                        t = p1; p1 = p2; p2 = t;
                    285:                        l = l1;
                    286:                }
                    287:                if ( l > wlen ) {
                    288:                        if ( w ) GC_free(w);
                    289:                        w = (MP *)MALLOC(l*sizeof(MP));
                    290:                        wlen = l;
                    291:                }
                    292:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    293:                        w[i] = m;
                    294:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    295:                        mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
                    296:                }
                    297:                bzero(w,l*sizeof(MP));
                    298:                *pr = s;
                    299:        }
                    300: }
                    301:
                    302: void weyl_mulmd(vl,mod,p1,p2,pr)
                    303: VL vl;
                    304: int mod;
                    305: DP p1,p2,*pr;
                    306: {
1.1       noro      307:        MP m;
                    308:        DP s,t,u;
1.2       noro      309:        int i,l,l1;
                    310:        static MP *w;
                    311:        static int wlen;
1.1       noro      312:
                    313:        if ( !p1 || !p2 )
                    314:                *pr = 0;
                    315:        else if ( OID(p1) <= O_P )
                    316:                mulmdc(vl,mod,p2,(P)p1,pr);
                    317:        else if ( OID(p2) <= O_P )
                    318:                mulmdc(vl,mod,p1,(P)p2,pr);
                    319:        else {
1.2       noro      320:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    321:                if ( l > wlen ) {
                    322:                        if ( w ) GC_free(w);
                    323:                        w = (MP *)MALLOC(l*sizeof(MP));
                    324:                        wlen = l;
                    325:                }
                    326:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    327:                        w[i] = m;
                    328:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    329:                        weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u;
1.1       noro      330:                }
1.2       noro      331:                bzero(w,l*sizeof(MP));
1.1       noro      332:                *pr = s;
                    333:        }
                    334: }
                    335:
                    336: void mulmdm(vl,mod,p,m0,pr)
                    337: VL vl;
                    338: int mod;
                    339: DP p;
                    340: MP m0;
                    341: DP *pr;
                    342: {
                    343:        MP m,mr,mr0;
                    344:        P c;
                    345:        DL d;
                    346:        int n,t;
                    347:        MQ q;
                    348:
                    349:        if ( !p )
                    350:                *pr = 0;
                    351:        else {
                    352:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    353:                        m; m = NEXT(m) ) {
                    354:                        NEXTMP(mr0,mr);
                    355:                        if ( NUM(C(m)) && NUM(c) ) {
                    356:                                t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
                    357:                                STOMQ(t,q); C(mr) = (P)q;
                    358:                        } else
                    359:                                mulmp(vl,mod,C(m),c,&C(mr));
                    360:                        adddl(n,m->dl,d,&mr->dl);
                    361:                }
                    362:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    363:                if ( *pr )
                    364:                        (*pr)->sugar = p->sugar + m0->dl->td;
                    365:        }
                    366: }
                    367:
1.2       noro      368: void weyl_mulmdm(vl,mod,p,m0,pr)
                    369: VL vl;
                    370: int mod;
                    371: DP p;
                    372: MP m0;
                    373: DP *pr;
                    374: {
                    375:        DP r,t,t1;
                    376:        MP m;
                    377:        int n,l,i;
                    378:        static MP *w;
                    379:        static int wlen;
                    380:
                    381:        if ( !p )
                    382:                *pr = 0;
                    383:        else {
                    384:                for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                    385:                if ( l > wlen ) {
                    386:                        if ( w ) GC_free(w);
                    387:                        w = (MP *)MALLOC(l*sizeof(MP));
                    388:                        wlen = l;
                    389:                }
                    390:                for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                    391:                        w[i] = m;
                    392:                for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
                    393:                        weyl_mulmmm(vl,mod,w[i],m0,n,&t);
                    394:                        addmd(vl,mod,r,t,&t1); r = t1;
                    395:                }
                    396:                bzero(w,l*sizeof(MP));
                    397:                if ( r )
                    398:                        r->sugar = p->sugar + m0->dl->td;
                    399:                *pr = r;
                    400:        }
                    401: }
                    402:
                    403: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    404:
                    405: void weyl_mulmmm(vl,mod,m0,m1,n,pr)
                    406: VL vl;
                    407: int mod;
                    408: MP m0,m1;
                    409: int n;
                    410: DP *pr;
                    411: {
                    412:        MP m,mr,mr0;
                    413:        MQ mq;
                    414:        DP r,t,t1;
                    415:        P c,c0,c1,cc;
                    416:        DL d,d0,d1;
1.3       noro      417:        int i,j,a,b,k,l,n2,s,min;
1.2       noro      418:        static int *tab;
                    419:        static int tablen;
                    420:
                    421:        if ( !m0 || !m1 )
                    422:                *pr = 0;
                    423:        else {
                    424:                c0 = C(m0); c1 = C(m1);
                    425:                mulmp(vl,mod,c0,c1,&c);
                    426:                d0 = m0->dl; d1 = m1->dl;
                    427:                n2 = n>>1;
                    428:                if ( n & 1 ) {
                    429:                        /* homogenized computation; dx-xd=h^2 */
                    430:                        /* offset of h-degree */
1.3       noro      431:                        NEWDL(d,n);
                    432:                        d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1.4       noro      433:                        NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0;
1.3       noro      434:                        MKDP(n,mr,r); r->sugar = d->td;
                    435:                } else
                    436:                        r = (DP)ONEM;
                    437:                for ( i = 0; i < n2; i++ ) {
                    438:                        a = d0->d[i]; b = d1->d[n2+i];
                    439:                        k = d0->d[n2+i]; l = d1->d[i];
                    440:                        /* degree of xi^a*(Di^k*xi^l)*Di^b */
                    441:                        s = a+k+l+b;
                    442:                        /* compute xi^a*(Di^k*xi^l)*Di^b */
                    443:                        min = MIN(k,l);
                    444:
                    445:                        if ( min+1 > tablen ) {
                    446:                                if ( tab ) GC_free(tab);
                    447:                                tab = (int *)MALLOC((min+1)*sizeof(int));
                    448:                                tablen = min+1;
                    449:                        }
                    450:                        mkwcm(k,l,mod,tab);
                    451:                        if ( n & 1 )
1.2       noro      452:                                for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3       noro      453:                                        NEXTMP(mr0,mr); NEWDL(d,n);
1.2       noro      454:                                        d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3       noro      455:                                        d->td = s;
                    456:                                        d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
                    457:                                        STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2       noro      458:                                }
1.3       noro      459:                        else
1.2       noro      460:                                for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3       noro      461:                                        NEXTMP(mr0,mr); NEWDL(d,n);
1.2       noro      462:                                        d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    463:                                        d->td = d->d[i]+d->d[n2+i]; /* XXX */
                    464:                                        s = MAX(s,d->td); /* XXX */
1.3       noro      465:                                        STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d;
1.2       noro      466:                                }
1.3       noro      467:                        bzero(tab,(min+1)*sizeof(int));
                    468:                        if ( mr0 )
                    469:                                NEXT(mr) = 0;
                    470:                        MKDP(n,mr0,t);
                    471:                        if ( t )
                    472:                                t->sugar = s;
                    473:                        comm_mulmd(vl,mod,r,t,&t1); r = t1;
                    474:                }
1.2       noro      475:                mulmdc(vl,mod,r,c,pr);
                    476:        }
                    477: }
                    478:
1.1       noro      479: void mulmdc(vl,mod,p,c,pr)
                    480: VL vl;
                    481: int mod;
                    482: DP p;
                    483: P c;
                    484: DP *pr;
                    485: {
                    486:        MP m,mr,mr0;
                    487:        int t;
                    488:        MQ q;
                    489:
                    490:        if ( !p || !c )
                    491:                *pr = 0;
                    492:        else {
                    493:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    494:                        NEXTMP(mr0,mr);
                    495:                        if ( NUM(C(m)) && NUM(c) ) {
                    496:                                t = dmar(((MQ)(C(m)))->cont,((MQ)c)->cont,0,mod);
                    497:                                STOMQ(t,q); C(mr) = (P)q;
                    498:                        } else
                    499:                                mulmp(vl,mod,C(m),c,&C(mr));
                    500:                        mr->dl = m->dl;
                    501:                }
                    502:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    503:                if ( *pr )
                    504:                        (*pr)->sugar = p->sugar;
                    505:        }
                    506: }
                    507:
                    508: void divsmdc(vl,mod,p,c,pr)
                    509: VL vl;
                    510: int mod;
                    511: DP p;
                    512: P c;
                    513: DP *pr;
                    514: {
                    515:        MP m,mr,mr0;
                    516:
                    517:        if ( !c )
                    518:                error("disvsdc : division by 0");
                    519:        else if ( !p )
                    520:                *pr = 0;
                    521:        else {
                    522:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    523:                        NEXTMP(mr0,mr); divsmp(vl,mod,C(m),c,&C(mr)); mr->dl = m->dl;
                    524:                }
                    525:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
                    526:                if ( *pr )
                    527:                        (*pr)->sugar = p->sugar;
                    528:        }
                    529: }
                    530:
                    531: void _mdtop(vl,mod,dvl,p,pr)
                    532: VL vl,dvl;
                    533: int mod;
                    534: DP p;
                    535: P *pr;
                    536: {
                    537:        int n,i;
                    538:        DL d;
                    539:        MP m;
                    540:        P r,s,t,u,w;
                    541:        Q q;
                    542:        MQ tq;
                    543:        VL tvl;
                    544:
                    545:        if ( !p )
                    546:                *pr = 0;
                    547:        else {
                    548:                for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
                    549:                        STOMQ(ITOS(C(m)),tq); t = (P)tq;
                    550:                        for ( i = 0, d = m->dl, tvl = dvl;
                    551:                                i < n; tvl = NEXT(tvl), i++ ) {
                    552:                                MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u);
                    553:                                mulmp(vl,mod,t,u,&w); t = w;
                    554:                        }
                    555:                        addmp(vl,mod,s,t,&u); s = u;
                    556:                }
                    557:                mptop(s,pr);
                    558:        }
                    559: }
                    560:
                    561: void _addmd(vl,mod,p1,p2,pr)
                    562: VL vl;
                    563: int mod;
                    564: DP p1,p2,*pr;
                    565: {
                    566:        int n;
                    567:        MP m1,m2,mr,mr0;
                    568:        int t;
                    569:
                    570:        if ( !p1 )
                    571:                *pr = p2;
                    572:        else if ( !p2 )
                    573:                *pr = p1;
                    574:        else {
                    575:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    576:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    577:                                case 0:
                    578:                                        t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
                    579:                                        if ( t < 0 )
                    580:                                                t += mod;
                    581:                                        if ( t ) {
                    582:                                                NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t);
                    583:                                        }
                    584:                                        m1 = NEXT(m1); m2 = NEXT(m2); break;
                    585:                                case 1:
                    586:                                        NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
                    587:                                        m1 = NEXT(m1); break;
                    588:                                case -1:
                    589:                                        NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
                    590:                                        m2 = NEXT(m2); break;
                    591:                        }
                    592:                if ( !mr0 )
                    593:                        if ( m1 )
                    594:                                mr0 = m1;
                    595:                        else if ( m2 )
                    596:                                mr0 = m2;
                    597:                        else {
                    598:                                *pr = 0;
                    599:                                return;
                    600:                        }
                    601:                else if ( m1 )
                    602:                        NEXT(mr) = m1;
                    603:                else if ( m2 )
                    604:                        NEXT(mr) = m2;
                    605:                else
                    606:                        NEXT(mr) = 0;
1.3       noro      607:                MKDP(NV(p1),mr0,*pr);
1.1       noro      608:                if ( *pr )
                    609:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
                    610:        }
                    611: }
                    612:
                    613: void _submd(vl,mod,p1,p2,pr)
                    614: VL vl;
                    615: int mod;
                    616: DP p1,p2,*pr;
                    617: {
                    618:        DP t;
                    619:
                    620:        if ( !p2 )
                    621:                *pr = p1;
                    622:        else {
                    623:                _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr);
                    624:        }
                    625: }
                    626:
                    627: void _chsgnmd(mod,p,pr)
                    628: int mod;
                    629: DP p,*pr;
                    630: {
                    631:        MP m,mr,mr0;
                    632:
                    633:        if ( !p )
                    634:                *pr = 0;
                    635:        else {
                    636:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    637:                        NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl;
                    638:                }
1.3       noro      639:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1       noro      640:                if ( *pr )
                    641:                        (*pr)->sugar = p->sugar;
                    642:        }
                    643: }
                    644:
                    645: void _mulmd(vl,mod,p1,p2,pr)
                    646: VL vl;
                    647: int mod;
                    648: DP p1,p2,*pr;
                    649: {
1.2       noro      650:        if ( !do_weyl )
                    651:                _comm_mulmd(vl,mod,p1,p2,pr);
                    652:        else
                    653:                _weyl_mulmd(vl,mod,p1,p2,pr);
                    654: }
                    655:
                    656: void _comm_mulmd(vl,mod,p1,p2,pr)
                    657: VL vl;
                    658: int mod;
                    659: DP p1,p2,*pr;
                    660: {
1.1       noro      661:        MP m;
                    662:        DP s,t,u;
1.2       noro      663:        int i,l,l1;
                    664:        static MP *w;
                    665:        static int wlen;
1.1       noro      666:
                    667:        if ( !p1 || !p2 )
                    668:                *pr = 0;
                    669:        else {
1.2       noro      670:                for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    671:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    672:                if ( l1 < l ) {
                    673:                        t = p1; p1 = p2; p2 = t;
                    674:                        l = l1;
                    675:                }
                    676:                if ( l > wlen ) {
                    677:                        if ( w ) GC_free(w);
                    678:                        w = (MP *)MALLOC(l*sizeof(MP));
                    679:                        wlen = l;
                    680:                }
                    681:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    682:                        w[i] = m;
                    683:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    684:                        _mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
1.1       noro      685:                }
1.2       noro      686:                bzero(w,l*sizeof(MP));
                    687:                *pr = s;
                    688:        }
                    689: }
                    690:
                    691: void _weyl_mulmd(vl,mod,p1,p2,pr)
                    692: VL vl;
                    693: int mod;
                    694: DP p1,p2,*pr;
                    695: {
                    696:        MP m;
                    697:        DP s,t,u;
                    698:        int i,l,l1;
                    699:        static MP *w;
                    700:        static int wlen;
                    701:
                    702:        if ( !p1 || !p2 )
                    703:                *pr = 0;
                    704:        else {
                    705:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    706:                if ( l > wlen ) {
                    707:                        if ( w ) GC_free(w);
                    708:                        w = (MP *)MALLOC(l*sizeof(MP));
                    709:                        wlen = l;
                    710:                }
                    711:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    712:                        w[i] = m;
                    713:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    714:                        _weyl_mulmdm(vl,mod,p1,w[i],&t); _addmd(vl,mod,s,t,&u); s = u;
                    715:                }
                    716:                bzero(w,l*sizeof(MP));
1.1       noro      717:                *pr = s;
                    718:        }
                    719: }
                    720:
                    721: void _mulmdm(vl,mod,p,m0,pr)
                    722: VL vl;
                    723: int mod;
                    724: DP p;
                    725: MP m0;
                    726: DP *pr;
                    727: {
                    728:        MP m,mr,mr0;
                    729:        DL d;
                    730:        int c,n,r;
                    731:
                    732:        if ( !p )
                    733:                *pr = 0;
                    734:        else {
                    735:                for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
                    736:                        m; m = NEXT(m) ) {
                    737:                        NEXTMP(mr0,mr);
                    738:                        C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
                    739:                        adddl(n,m->dl,d,&mr->dl);
                    740:                }
1.3       noro      741:                NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1.1       noro      742:                if ( *pr )
                    743:                        (*pr)->sugar = p->sugar + m0->dl->td;
                    744:        }
                    745: }
                    746:
1.2       noro      747: void _weyl_mulmdm(vl,mod,p,m0,pr)
                    748: VL vl;
                    749: int mod;
                    750: DP p;
                    751: MP m0;
                    752: DP *pr;
                    753: {
                    754:        DP r,t,t1;
                    755:        MP m;
                    756:        int n,l,i;
                    757:        static MP *w;
                    758:        static int wlen;
                    759:
                    760:        if ( !p )
                    761:                *pr = 0;
                    762:        else {
                    763:                for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                    764:                if ( l > wlen ) {
                    765:                        if ( w ) GC_free(w);
                    766:                        w = (MP *)MALLOC(l*sizeof(MP));
                    767:                        wlen = l;
                    768:                }
                    769:                for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                    770:                        w[i] = m;
                    771:                for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) {
                    772:                        _weyl_mulmmm(vl,mod,w[i],m0,n,&t);
                    773:                        _addmd(vl,mod,r,t,&t1); r = t1;
                    774:                }
                    775:                bzero(w,l*sizeof(MP));
                    776:                if ( r )
                    777:                        r->sugar = p->sugar + m0->dl->td;
                    778:                *pr = r;
                    779:        }
                    780: }
                    781:
                    782: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    783:
                    784: void _weyl_mulmmm(vl,mod,m0,m1,n,pr)
                    785: VL vl;
                    786: int mod;
                    787: MP m0,m1;
                    788: int n;
                    789: DP *pr;
                    790: {
                    791:        MP m,mr,mr0;
                    792:        DP r,t,t1;
                    793:        int c,c0,c1,cc;
                    794:        DL d,d0,d1;
                    795:        int i,j,a,b,k,l,n2,s,min,h;
                    796:        static int *tab;
                    797:        static int tablen;
                    798:
                    799:        if ( !m0 || !m1 )
                    800:                *pr = 0;
                    801:        else {
                    802:                c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
                    803:                c = dmar(c0,c1,0,mod);
                    804:                d0 = m0->dl; d1 = m1->dl;
                    805:                n2 = n>>1;
                    806:
1.3       noro      807:                NEWDL(d,n);
                    808:                if ( n & 1 )
1.2       noro      809:                        /* offset of h-degree */
1.3       noro      810:                        d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
                    811:                else
                    812:                        d->td = 0;
1.4       noro      813:                NEWMP(mr); mr->c = STOI(c); mr->dl = d; NEXT(mr) = 0;
1.3       noro      814:                MKDP(n,mr,r); r->sugar = d->td;
                    815:
                    816:                /* homogenized computation; dx-xd=h^2 */
                    817:                for ( i = 0; i < n2; i++ ) {
                    818:                        a = d0->d[i]; b = d1->d[n2+i];
                    819:                        k = d0->d[n2+i]; l = d1->d[i];
                    820:                        /* degree of xi^a*(Di^k*xi^l)*Di^b */
                    821:                        s = a+k+l+b;
                    822:                        /* compute xi^a*(Di^k*xi^l)*Di^b */
                    823:                        min = MIN(k,l);
                    824:
                    825:                        if ( min+1 > tablen ) {
                    826:                                if ( tab ) GC_free(tab);
                    827:                                tab = (int *)MALLOC((min+1)*sizeof(int));
                    828:                                tablen = min+1;
                    829:                        }
                    830:                        mkwcm(k,l,mod,tab);
                    831:                        if ( n & 1 )
1.2       noro      832:                                for ( mr0 = 0, j = 0; j <= min; j++ ) {
1.3       noro      833:                                        NEXTMP(mr0,mr); NEWDL(d,n);
1.2       noro      834:                                        d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
1.3       noro      835:                                        d->td = s;
                    836:                                        d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
                    837:                                        mr->c = STOI(tab[j]); mr->dl = d;
1.2       noro      838:                                }
1.3       noro      839:                        else
1.2       noro      840:                                for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) {
1.3       noro      841:                                        NEXTMP(mr0,mr); NEWDL(d,n);
1.2       noro      842:                                        d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    843:                                        d->td = d->d[i]+d->d[n2+i]; /* XXX */
                    844:                                        s = MAX(s,d->td); /* XXX */
1.3       noro      845:                                        mr->c = STOI(tab[j]); mr->dl = d;
1.2       noro      846:                                }
1.3       noro      847:                        bzero(tab,(min+1)*sizeof(int));
                    848:                        if ( mr0 )
                    849:                                NEXT(mr) = 0;
                    850:                        MKDP(n,mr0,t);
                    851:                        if ( t )
                    852:                                t->sugar = s;
                    853:                        _comm_mulmd(vl,mod,r,t,&t1); r = t1;
1.2       noro      854:                }
                    855:                *pr = r;
                    856:        }
                    857: }
                    858:
1.1       noro      859: void _dtop_mod(vl,dvl,p,pr)
                    860: VL vl,dvl;
                    861: DP p;
                    862: P *pr;
                    863: {
                    864:        int n,i;
                    865:        DL d;
                    866:        MP m;
                    867:        P r,s,t,u,w;
                    868:        Q q;
                    869:        VL tvl;
                    870:
                    871:        if ( !p )
                    872:                *pr = 0;
                    873:        else {
                    874:                for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) {
                    875:                        i = ITOS(m->c); STOQ(i,q); t = (P)q;
                    876:                        for ( i = 0, d = m->dl, tvl = dvl;
                    877:                                i < n; tvl = NEXT(tvl), i++ ) {
                    878:                                MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
                    879:                                mulp(vl,t,u,&w); t = w;
                    880:                        }
                    881:                        addp(vl,s,t,&u); s = u;
                    882:                }
                    883:                *pr = s;
                    884:        }
                    885: }
                    886:
                    887: void _dp_red_mod(p1,p2,mod,rp)
                    888: DP p1,p2;
                    889: int mod;
                    890: DP *rp;
                    891: {
                    892:        int i,n;
                    893:        DL d1,d2,d;
                    894:        MP m;
                    895:        DP t,s;
                    896:        int c,c1;
                    897:
                    898:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    899:        NEWDL(d,n); d->td = d1->td - d2->td;
                    900:        for ( i = 0; i < n; i++ )
                    901:                d->d[i] = d1->d[i]-d2->d[i];
                    902:        c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
                    903:        NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
                    904:        MKDP(n,m,s); s->sugar = d->td;
1.2       noro      905:        _mulmd(CO,mod,s,p2,&t); _addmd(CO,mod,p1,t,rp);
1.1       noro      906: }
                    907:
                    908: void _dp_mod(p,mod,subst,rp)
                    909: DP p;
                    910: int mod;
                    911: NODE subst;
                    912: DP *rp;
                    913: {
                    914:        MP m,mr,mr0;
                    915:        P t,s;
                    916:        NODE tn;
                    917:
                    918:        if ( !p )
                    919:                *rp = 0;
                    920:        else {
                    921:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    922:                        for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) {
                    923:                                substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t);
                    924:                                s = t;
                    925:                        }
                    926:                        ptomp(mod,s,&t);
                    927:                        if ( t ) {
                    928:                                NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl;
                    929:                        }
                    930:                }
                    931:                if ( mr0 ) {
                    932:                        NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    933:                } else
                    934:                        *rp = 0;
                    935:        }
                    936: }
                    937:
1.2       noro      938: void _dp_monic(p,mod,rp)
                    939: DP p;
                    940: int mod;
                    941: DP *rp;
                    942: {
                    943:        MP m,mr,mr0;
                    944:        int c,c1;
                    945:        NODE tn;
                    946:
                    947:        if ( !p )
                    948:                *rp = 0;
                    949:        else {
                    950:                c = invm(ITOS(BDY(p)->c),mod);
                    951:                for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
                    952:                        c1 = dmar(ITOS(m->c),c,0,mod);
                    953:                        NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl;
                    954:                }
                    955:                NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
                    956:        }
                    957: }
                    958:
1.1       noro      959: void _dp_sp_mod(p1,p2,mod,rp)
                    960: DP p1,p2;
                    961: int mod;
                    962: DP *rp;
                    963: {
                    964:        int i,n,td;
                    965:        int *w;
                    966:        DL d1,d2,d;
                    967:        MP m;
                    968:        DP t,s,u;
                    969:
                    970:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                    971:        w = (int *)ALLOCA(n*sizeof(int));
                    972:        for ( i = 0, td = 0; i < n; i++ ) {
                    973:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
                    974:        }
                    975:        NEWDL(d,n); d->td = td - d1->td;
                    976:        for ( i = 0; i < n; i++ )
                    977:                d->d[i] = w[i] - d1->d[i];
                    978:        NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2       noro      979:        MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,&t);
1.1       noro      980:        NEWDL(d,n); d->td = td - d2->td;
                    981:        for ( i = 0; i < n; i++ )
                    982:                d->d[i] = w[i] - d2->d[i];
                    983:        NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
1.2       noro      984:        MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,&u);
1.1       noro      985:        _addmd(CO,mod,t,u,rp);
                    986: }
                    987:
                    988: void _dp_sp_component_mod(p1,p2,mod,f1,f2)
                    989: DP p1,p2;
                    990: int mod;
                    991: DP *f1,*f2;
                    992: {
                    993:        int i,n,td;
                    994:        int *w;
                    995:        DL d1,d2,d;
                    996:        MP m;
                    997:        DP t,s,u;
                    998:
                    999:        n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
                   1000:        w = (int *)ALLOCA(n*sizeof(int));
                   1001:        for ( i = 0, td = 0; i < n; i++ ) {
                   1002:                w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
                   1003:        }
                   1004:        NEWDL(d,n); d->td = td - d1->td;
                   1005:        for ( i = 0; i < n; i++ )
                   1006:                d->d[i] = w[i] - d1->d[i];
                   1007:        NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
1.2       noro     1008:        MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p1,f1);
1.1       noro     1009:        NEWDL(d,n); d->td = td - d2->td;
                   1010:        for ( i = 0; i < n; i++ )
                   1011:                d->d[i] = w[i] - d2->d[i];
                   1012:        NEWMP(m); m->dl = d; m->c = BDY(p1)->c; NEXT(m) = 0;
1.2       noro     1013:        MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,s,p2,f2);
1.1       noro     1014: }
                   1015:
                   1016: void _printdp(d)
                   1017: DP d;
                   1018: {
                   1019:        int n,i;
                   1020:        MP m;
                   1021:        DL dl;
                   1022:
                   1023:        if ( !d ) {
                   1024:                printf("0"); return;
                   1025:        }
                   1026:        for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) {
                   1027:                printf("%d*<<",ITOS(m->c));
                   1028:                for ( i = 0, dl = m->dl; i < n-1; i++ )
                   1029:                        printf("%d,",dl->d[i]);
                   1030:                printf("%d",dl->d[i]);
                   1031:                printf(">>");
                   1032:                if ( NEXT(m) )
                   1033:                        printf("+");
                   1034:        }
                   1035: }

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