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

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

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