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

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.17    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.16 2017/08/31 02:36:21 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "inline.h"
                     52:
                     53: extern int (*cmpdl)();
                     54: extern int do_weyl;
                     55:
                     56: MP _mp_free_list;
                     57: DP _dp_free_list;
                     58: DL _dl_free_list;
                     59: int current_dl_length;
1.6       noro       60:
1.11      noro       61: void GC_gcollect();
                     62:
1.6       noro       63: void _free_private_storage()
                     64: {
1.17    ! noro       65:   _mp_free_list = 0;
        !            66:   _dp_free_list = 0;
        !            67:   _dl_free_list = 0;
        !            68:   GC_gcollect();
1.6       noro       69: }
1.1       noro       70:
                     71: void _DL_alloc()
                     72: {
1.17    ! noro       73:   int *p;
        !            74:   int i,dl_len;
        !            75:   static int DL_alloc_count;
1.1       noro       76:
1.17    ! noro       77: /*  fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
        !            78:   dl_len = (current_dl_length+1);
1.14      ohara      79: #if SIZEOF_LONG == 8
1.17    ! noro       80:   if ( dl_len & 1 )
        !            81:     dl_len += 1;
1.8       noro       82: #endif
1.17    ! noro       83:   for ( i = 0; i < 128; i++, p += dl_len ) {
        !            84:     p = (int *)MALLOC(dl_len*sizeof(int));
        !            85:     *(DL *)p = _dl_free_list;
        !            86:     _dl_free_list = (DL)p;
        !            87:   }
1.1       noro       88: }
                     89:
                     90: void _MP_alloc()
                     91: {
1.17    ! noro       92:   MP p;
        !            93:   int i;
        !            94:   static int MP_alloc_count;
        !            95:
        !            96: /*  fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */
        !            97:   for ( i = 0; i < 1024; i++ ) {
        !            98:     p = (MP)MALLOC(sizeof(struct oMP));
        !            99:     p->next = _mp_free_list; _mp_free_list = p;
        !           100:   }
1.1       noro      101: }
1.17    ! noro      102:
1.1       noro      103: void _DP_alloc()
                    104: {
1.17    ! noro      105:   DP p;
        !           106:   int i;
        !           107:   static int DP_alloc_count;
        !           108:
        !           109: /*  fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */
        !           110:   for ( i = 0; i < 1024; i++ ) {
        !           111:     p = (DP)MALLOC(sizeof(struct oDP));
        !           112:     p->body = (MP)_dp_free_list; _dp_free_list = p;
        !           113:   }
1.1       noro      114: }
                    115:
                    116: /* merge p1 and p2 into pr */
                    117:
1.11      noro      118: void _addmd_destructive(int mod,DP p1,DP p2,DP *pr)
1.1       noro      119: {
1.17    ! noro      120:   int n;
        !           121:   MP m1,m2,mr,mr0,s;
        !           122:   int t;
        !           123:
        !           124:   if ( !p1 )
        !           125:     *pr = p2;
        !           126:   else if ( !p2 )
        !           127:     *pr = p1;
        !           128:   else {
        !           129:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
        !           130:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
        !           131:         case 0:
        !           132:           t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
        !           133:           if ( t < 0 )
        !           134:             t += mod;
        !           135:           s = m1; m1 = NEXT(m1);
        !           136:           if ( t ) {
        !           137:             _NEXTMP2(mr0,mr,s); C(mr) = (Obj)STOI(t);
        !           138:           } else {
        !           139:             _FREEDL(s->dl); _FREEMP(s);
        !           140:           }
        !           141:           s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
        !           142:           break;
        !           143:         case 1:
        !           144:           s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
        !           145:           break;
        !           146:         case -1:
        !           147:           s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
        !           148:           break;
        !           149:       }
        !           150:     if ( !mr0 )
        !           151:       if ( m1 )
        !           152:         mr0 = m1;
        !           153:       else if ( m2 )
        !           154:         mr0 = m2;
        !           155:       else {
        !           156:         *pr = 0;
        !           157:         return;
        !           158:       }
        !           159:     else if ( m1 )
        !           160:       NEXT(mr) = m1;
        !           161:     else if ( m2 )
        !           162:       NEXT(mr) = m2;
        !           163:     else
        !           164:       NEXT(mr) = 0;
        !           165:     _MKDP(NV(p1),mr0,*pr);
        !           166:     if ( *pr )
        !           167:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
        !           168:     _FREEDP(p1); _FREEDP(p2);
        !           169:   }
1.1       noro      170: }
                    171:
1.11      noro      172: void _mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      173: {
1.17    ! noro      174:   if ( !do_weyl )
        !           175:     _comm_mulmd_dup(mod,p1,p2,pr);
        !           176:   else
        !           177:     _weyl_mulmd_dup(mod,p1,p2,pr);
1.1       noro      178: }
                    179:
1.11      noro      180: void _comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      181: {
1.17    ! noro      182:   MP m;
        !           183:   DP s,t,u;
        !           184:   int i,l,l1;
        !           185:   static MP *w;
        !           186:   static int wlen;
        !           187:
        !           188:   if ( !p1 || !p2 )
        !           189:     *pr = 0;
        !           190:   else {
        !           191:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
        !           192:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
        !           193:     if ( l1 < l ) {
        !           194:       t = p1; p1 = p2; p2 = t;
        !           195:       l = l1;
        !           196:     }
        !           197:     if ( l > wlen ) {
        !           198:       if ( w ) GCFREE(w);
        !           199:       w = (MP *)MALLOC(l*sizeof(MP));
        !           200:       wlen = l;
        !           201:     }
        !           202:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
        !           203:       w[i] = m;
        !           204:     for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           205:       _mulmdm_dup(mod,p1,w[i],&t); _addmd_destructive(mod,s,t,&u); s = u;
        !           206:     }
        !           207:     bzero(w,l*sizeof(MP));
        !           208:     *pr = s;
        !           209:   }
1.1       noro      210: }
                    211:
1.11      noro      212: void _weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      213: {
1.17    ! noro      214:   MP m;
        !           215:   DP s,t,u;
        !           216:   int i,l;
        !           217:   static MP *w;
        !           218:   static int wlen;
        !           219:
        !           220:   if ( !p1 || !p2 )
        !           221:     *pr = 0;
        !           222:   else {
        !           223:     for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
        !           224:     if ( l > wlen ) {
        !           225:       if ( w ) GCFREE(w);
        !           226:       w = (MP *)MALLOC(l*sizeof(MP));
        !           227:       wlen = l;
        !           228:     }
        !           229:     for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
        !           230:       w[i] = m;
        !           231:     for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           232:       _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
        !           233:     }
        !           234:     bzero(w,l*sizeof(MP));
        !           235:     *pr = s;
        !           236:   }
1.1       noro      237: }
                    238:
1.11      noro      239: void _mulmdm_dup(int mod,DP p,MP m0,DP *pr)
1.1       noro      240: {
1.17    ! noro      241:   MP m,mr,mr0;
        !           242:   DL d,dt,dm;
        !           243:   int c,n,i,c1,c2;
        !           244:   int *pt,*p1,*p2;
        !           245:
        !           246:   if ( !p )
        !           247:     *pr = 0;
        !           248:   else {
        !           249:     for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
        !           250:       m; m = NEXT(m) ) {
        !           251:       _NEXTMP(mr0,mr);
        !           252:       c1 = ITOS(C(m));
        !           253:       DMAR(c1,c,0,mod,c2);
        !           254:       C(mr) = (Obj)STOI(c2);
        !           255:       _NEWDL_NOINIT(dt,n); mr->dl = dt;
        !           256:       dm = m->dl;
        !           257:       dt->td = d->td + dm->td;
        !           258:       for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
        !           259:         *pt++ = *p1++ + *p2++;
        !           260:     }
        !           261:     NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
        !           262:     if ( *pr )
        !           263:       (*pr)->sugar = p->sugar + m0->dl->td;
        !           264:   }
1.1       noro      265: }
                    266:
1.11      noro      267: void _weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
1.1       noro      268: {
1.17    ! noro      269:   DP r,t,t1;
        !           270:   MP m;
        !           271:   DL d0;
        !           272:   int n,n2,l,i,j,tlen;
        !           273:   static MP *w,*psum;
        !           274:   static struct cdlm *tab;
        !           275:   static int wlen;
        !           276:   static int rtlen;
        !           277:
        !           278:   if ( !p )
        !           279:     *pr = 0;
        !           280:   else {
        !           281:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
        !           282:     if ( l > wlen ) {
        !           283:       if ( w ) GCFREE(w);
        !           284:       w = (MP *)MALLOC(l*sizeof(MP));
        !           285:       wlen = l;
        !           286:     }
        !           287:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
        !           288:       w[i] = m;
        !           289:     n = NV(p); n2 = n>>1;
        !           290:     d0 = m0->dl;
        !           291:
        !           292:     for ( i = 0, tlen = 1; i < n2; i++ )
        !           293:       tlen *= d0->d[n2+i]+1;
        !           294:     if ( tlen > rtlen ) {
        !           295:       if ( tab ) GCFREE(tab);
        !           296:       if ( psum ) GCFREE(psum);
        !           297:       rtlen = tlen;
        !           298:       tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
        !           299:       psum = (MP *)MALLOC(rtlen*sizeof(MP));
        !           300:     }
        !           301:     bzero(psum,tlen*sizeof(MP));
        !           302:     for ( i = l-1; i >= 0; i-- ) {
        !           303:       bzero(tab,tlen*sizeof(struct cdlm));
        !           304:       _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
        !           305:       for ( j = 0; j < tlen; j++ ) {
        !           306:         if ( tab[j].c ) {
        !           307:           _NEWMP(m); m->dl = tab[j].d;
        !           308:           C(m) = (Obj)STOI(tab[j].c); NEXT(m) = psum[j];
        !           309:           psum[j] = m;
        !           310:         }
        !           311:       }
        !           312:     }
        !           313:     for ( j = tlen-1, r = 0; j >= 0; j-- )
        !           314:       if ( psum[j] ) {
        !           315:         _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
        !           316:       }
        !           317:     if ( r )
        !           318:       r->sugar = p->sugar + m0->dl->td;
        !           319:     *pr = r;
        !           320:   }
1.1       noro      321: }
1.4       noro      322:
1.1       noro      323: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    324:
1.11      noro      325: void _weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen)
1.1       noro      326: {
1.17    ! noro      327:   int c,c0,c1;
        !           328:   DL d,d0,d1,dt;
        !           329:   int i,j,a,b,k,l,n2,s,min,curlen;
        !           330:   struct cdlm *p;
        !           331:   static int *ctab;
        !           332:   static struct cdlm *tab;
        !           333:   static int tablen;
        !           334:   static struct cdlm *tmptab;
        !           335:   static int tmptablen;
        !           336:
        !           337:   if ( !m0 || !m1 ) {
        !           338:     rtab[0].c = 0;
        !           339:     rtab[0].d = 0;
        !           340:     return;
        !           341:   }
        !           342:   c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
        !           343:   c = dmar(c0,c1,0,mod);
        !           344:   d0 = m0->dl; d1 = m1->dl;
        !           345:   n2 = n>>1;
        !           346:   curlen = 1;
        !           347:
        !           348:   _NEWDL(d,n);
        !           349:   if ( n & 1 )
        !           350:     /* offset of h-degree */
        !           351:     d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
        !           352:   else
        !           353:     d->td = 0;
        !           354:   rtab[0].c = c;
        !           355:   rtab[0].d = d;
        !           356:
        !           357:   if ( rtablen > tmptablen ) {
        !           358:     if ( tmptab ) GCFREE(tmptab);
        !           359:     tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
        !           360:     tmptablen = rtablen;
        !           361:   }
        !           362:
        !           363:   for ( i = 0; i < n2; i++ ) {
        !           364:     a = d0->d[i]; b = d1->d[n2+i];
        !           365:     k = d0->d[n2+i]; l = d1->d[i];
        !           366:
        !           367:     /* degree of xi^a*(Di^k*xi^l)*Di^b */
        !           368:     a += l;
        !           369:     b += k;
        !           370:     s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
        !           371:
        !           372:     if ( !k || !l ) {
        !           373:       for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
        !           374:         if ( p->c ) {
        !           375:           dt = p->d;
        !           376:           dt->d[i] = a;
        !           377:           dt->d[n2+i] = b;
        !           378:           dt->td += s;
        !           379:         }
        !           380:       }
        !           381:       curlen *= k+1;
        !           382:       continue;
        !           383:     }
        !           384:     if ( k+1 > tablen ) {
        !           385:       if ( tab ) GCFREE(tab);
        !           386:       if ( ctab ) GCFREE(ctab);
        !           387:       tablen = k+1;
        !           388:       tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
        !           389:       ctab = (int *)MALLOC(tablen*sizeof(int));
        !           390:     }
        !           391:     /* compute xi^a*(Di^k*xi^l)*Di^b */
        !           392:     min = MIN(k,l);
        !           393:     mkwcm(k,l,mod,ctab);
        !           394:     bzero(tab,(k+1)*sizeof(struct cdlm));
        !           395:     /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
        !           396:     if ( n & 1 )
        !           397:       for ( j = 0; j <= min; j++ ) {
        !           398:         _NEWDL(d,n);
        !           399:         d->d[i] = a-j; d->d[n2+i] = b-j;
        !           400:         d->td = s;
        !           401:         d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
        !           402:         tab[j].d = d;
        !           403:         tab[j].c = ctab[j];
        !           404:       }
        !           405:     else
        !           406:       for ( j = 0; j <= min; j++ ) {
        !           407:         _NEWDL(d,n);
        !           408:         d->d[i] = a-j; d->d[n2+i] = b-j;
        !           409:         d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i);  /* XXX */
        !           410:         tab[j].d = d;
        !           411:         tab[j].c = ctab[j];
        !           412:       }
1.4       noro      413: #if 0
1.17    ! noro      414:     _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
        !           415:     for ( j = 0; j < curlen; j++ )
        !           416:       if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
        !           417:     for ( j = 0; j <= min; j++ )
        !           418:       if ( tab[j].d ) { _FREEDL(tab[j].d); }
        !           419:     curlen *= k+1;
        !           420:     bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
1.4       noro      421: #else
1.17    ! noro      422:     _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
        !           423:     for ( j = 0; j <= min; j++ )
        !           424:       if ( tab[j].d ) { _FREEDL(tab[j].d); }
        !           425:     curlen *= k+1;
1.4       noro      426: #endif
1.17    ! noro      427:   }
1.4       noro      428: }
                    429:
                    430: /* direct product of two cdlm tables
                    431:   rt[] = [
                    432:     t[0]*t1[0],...,t[n-1]*t1[0],
                    433:     t[0]*t1[1],...,t[n-1]*t1[1],
                    434:     ...
                    435:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                    436:   ]
                    437: */
                    438:
1.11      noro      439: void _comm_mulmd_tab(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1,struct cdlm *rt)
1.4       noro      440: {
1.17    ! noro      441:   int i,j;
        !           442:   struct cdlm *p;
        !           443:   int c;
        !           444:   DL d;
        !           445:
        !           446:   bzero(rt,n*n1*sizeof(struct cdlm));
        !           447:   for ( j = 0, p = rt; j < n1; j++ ) {
        !           448:     c = t1[j].c;
        !           449:     d = t1[j].d;
        !           450:     if ( !c )
        !           451:       break;
        !           452:     for ( i = 0; i < n; i++, p++ ) {
        !           453:       if ( t[i].c ) {
        !           454:         p->c = dmar(t[i].c,c,0,mod);
        !           455:         _adddl_dup(nv,t[i].d,d,&p->d);
        !           456:       }
        !           457:     }
        !           458:   }
1.4       noro      459: }
                    460:
1.11      noro      461: void _comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1)
1.4       noro      462: {
1.17    ! noro      463:   int i,j;
        !           464:   struct cdlm *p;
        !           465:   int c;
        !           466:   DL d;
        !           467:
        !           468:   for ( j = 1, p = t+n; j < n1; j++ ) {
        !           469:     c = t1[j].c;
        !           470:     d = t1[j].d;
        !           471:     if ( !c )
        !           472:       break;
        !           473:     for ( i = 0; i < n; i++, p++ ) {
        !           474:       if ( t[i].c ) {
        !           475:         p->c = dmar(t[i].c,c,0,mod);
        !           476:         _adddl_dup(nv,t[i].d,d,&p->d);
        !           477:       }
        !           478:     }
        !           479:   }
        !           480:   c = t1[0].c;
        !           481:   d = t1[0].d;
        !           482:   for ( i = 0, p = t; i < n; i++, p++ )
        !           483:     if ( t[i].c ) {
        !           484:       p->c = dmar(t[i].c,c,0,mod);
        !           485:       /* t[i].d += d */
        !           486:       adddl_destructive(nv,t[i].d,d);
        !           487:     }
1.1       noro      488: }
                    489:
1.11      noro      490: void dlto_dl(DL d,DL *dr)
1.1       noro      491: {
1.17    ! noro      492:   int i,n;
        !           493:   DL t;
1.1       noro      494:
1.17    ! noro      495:   n = current_dl_length;
        !           496:   _NEWDL(t,n); *dr = t;
        !           497:   t->td = d->td;
        !           498:   for ( i = 0; i < n; i++ )
        !           499:     t->d[i] = d->d[i];
1.1       noro      500: }
                    501:
1.11      noro      502: void _dltodl(DL d,DL *dr)
1.1       noro      503: {
1.17    ! noro      504:   int i,n;
        !           505:   DL t;
1.1       noro      506:
1.17    ! noro      507:   n = current_dl_length;
        !           508:   NEWDL(t,n); *dr = t;
        !           509:   t->td = d->td;
        !           510:   for ( i = 0; i < n; i++ )
        !           511:     t->d[i] = d->d[i];
1.1       noro      512: }
                    513:
1.11      noro      514: void _adddl_dup(int n,DL d1,DL d2,DL *dr)
1.1       noro      515: {
1.17    ! noro      516:   DL dt;
        !           517:   int i;
1.1       noro      518:
1.17    ! noro      519:   _NEWDL(dt,n);
        !           520:   *dr = dt;
        !           521:   dt->td = d1->td + d2->td;
        !           522:   for ( i = 0; i < n; i++ )
        !           523:     dt->d[i] = d1->d[i]+d2->d[i];
1.4       noro      524: }
                    525:
1.11      noro      526: void _free_dlarray(DL *a,int n)
1.4       noro      527: {
1.17    ! noro      528:   int i;
1.4       noro      529:
1.17    ! noro      530:   for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
1.1       noro      531: }
                    532:
1.11      noro      533: void _free_dp(DP f)
1.1       noro      534: {
1.17    ! noro      535:   MP m,s;
1.1       noro      536:
1.17    ! noro      537:   if ( !f )
        !           538:     return;
        !           539:   m = f->body;
        !           540:   while ( m ) {
        !           541:     s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
        !           542:   }
        !           543:   _FREEDP(f);
1.5       noro      544: }
                    545:
1.11      noro      546: void dpto_dp(DP p,DP *r)
1.5       noro      547: {
1.17    ! noro      548:   MP m,mr0,mr;
        !           549:   DL t;
1.5       noro      550:
1.17    ! noro      551:   if ( !p )
        !           552:     *r = 0;
        !           553:   else {
        !           554:     /* XXX : dummy call to set current_dl_length */
        !           555:     _NEWDL_NOINIT(t,NV(p));
        !           556:
        !           557:     for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
        !           558:       _NEXTMP(mr0,mr);
        !           559:       dlto_dl(m->dl,&mr->dl);
        !           560:       mr->c = m->c;
        !           561:     }
        !           562:     NEXT(mr) = 0;
        !           563:     _MKDP(p->nv,mr0,*r);
        !           564:     (*r)->sugar = p->sugar;
        !           565:   }
1.1       noro      566: }
                    567:
1.11      noro      568: void _dptodp(DP p,DP *r)
1.1       noro      569: {
1.17    ! noro      570:   MP m,mr0,mr;
1.1       noro      571:
1.17    ! noro      572:   if ( !p )
        !           573:     *r = 0;
        !           574:   else {
        !           575:     for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
        !           576:       NEXTMP(mr0,mr);
        !           577:       _dltodl(m->dl,&mr->dl);
        !           578:       mr->c = m->c;
        !           579:     }
        !           580:     NEXT(mr) = 0;
        !           581:     MKDP(p->nv,mr0,*r);
        !           582:     (*r)->sugar = p->sugar;
        !           583:   }
1.1       noro      584: }
                    585:
1.5       noro      586: /*
                    587:  * destructive merge of two list
                    588:  *
                    589:  * p1, p2 : list of DL
                    590:  * return : a merged list
                    591:  */
                    592:
1.11      noro      593: NODE _symb_merge(NODE m1,NODE m2,int n)
1.5       noro      594: {
1.17    ! noro      595:   NODE top,prev,cur,m,t;
1.5       noro      596:
1.17    ! noro      597:   if ( !m1 )
        !           598:     return m2;
        !           599:   else if ( !m2 )
        !           600:     return m1;
        !           601:   else {
        !           602:     switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
        !           603:       case 0:
        !           604:         top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
        !           605:         break;
        !           606:       case 1:
        !           607:         top = m1; m = m2;
        !           608:         break;
        !           609:       case -1:
        !           610:         top = m2; m = m1;
        !           611:         break;
        !           612:     }
        !           613:     prev = top; cur = NEXT(top);
        !           614:     /* BDY(prev) > BDY(m) always holds */
        !           615:     while ( cur && m ) {
        !           616:       switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
        !           617:         case 0:
        !           618:           _FREEDL(BDY(m)); m = NEXT(m);
        !           619:           prev = cur; cur = NEXT(cur);
        !           620:           break;
        !           621:         case 1:
        !           622:           t = NEXT(cur); NEXT(cur) = m; m = t;
        !           623:           prev = cur; cur = NEXT(cur);
        !           624:           break;
        !           625:         case -1:
        !           626:           NEXT(prev) = m; m = cur;
        !           627:           prev = NEXT(prev); cur = NEXT(prev);
        !           628:           break;
        !           629:       }
        !           630:     }
        !           631:     if ( !cur )
        !           632:       NEXT(prev) = m;
        !           633:     return top;
        !           634:   }
1.9       noro      635: }
                    636:
                    637: /* merge p1 and p2 into pr */
                    638:
1.11      noro      639: void _addd_destructive(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      640: {
1.17    ! noro      641:   int n;
        !           642:   MP m1,m2,mr,mr0,s;
        !           643:   P t;
        !           644:
        !           645:   if ( !p1 )
        !           646:     *pr = p2;
        !           647:   else if ( !p2 )
        !           648:     *pr = p1;
        !           649:   else {
        !           650:     for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
        !           651:       switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
        !           652:         case 0:
        !           653:           addp(vl,(P)C(m1),(P)C(m2),&t);
        !           654:           s = m1; m1 = NEXT(m1);
        !           655:           if ( t ) {
        !           656:             _NEXTMP2(mr0,mr,s); C(mr) = (Obj)t;
        !           657:           } else {
        !           658:             _FREEDL(s->dl); _FREEMP(s);
        !           659:           }
        !           660:           s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
        !           661:           break;
        !           662:         case 1:
        !           663:           s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
        !           664:           break;
        !           665:         case -1:
        !           666:           s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
        !           667:           break;
        !           668:       }
        !           669:     if ( !mr0 )
        !           670:       if ( m1 )
        !           671:         mr0 = m1;
        !           672:       else if ( m2 )
        !           673:         mr0 = m2;
        !           674:       else {
        !           675:         *pr = 0;
        !           676:         return;
        !           677:       }
        !           678:     else if ( m1 )
        !           679:       NEXT(mr) = m1;
        !           680:     else if ( m2 )
        !           681:       NEXT(mr) = m2;
        !           682:     else
        !           683:       NEXT(mr) = 0;
        !           684:     _MKDP(NV(p1),mr0,*pr);
        !           685:     if ( *pr )
        !           686:       (*pr)->sugar = MAX(p1->sugar,p2->sugar);
        !           687:     _FREEDP(p1); _FREEDP(p2);
        !           688:   }
1.9       noro      689: }
                    690:
1.11      noro      691: void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      692: {
1.17    ! noro      693:   if ( !do_weyl )
        !           694:     _comm_muld_dup(vl,p1,p2,pr);
        !           695:   else
        !           696:     _weyl_muld_dup(vl,p1,p2,pr);
1.9       noro      697: }
                    698:
1.11      noro      699: void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      700: {
1.17    ! noro      701:   MP m;
        !           702:   DP s,t,u;
        !           703:   int i,l,l1;
        !           704:   static MP *w;
        !           705:   static int wlen;
        !           706:
        !           707:   if ( !p1 || !p2 )
        !           708:     *pr = 0;
        !           709:   else {
        !           710:     for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
        !           711:     for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
        !           712:     if ( l1 < l ) {
        !           713:       t = p1; p1 = p2; p2 = t;
        !           714:       l = l1;
        !           715:     }
        !           716:     if ( l > wlen ) {
        !           717:       if ( w ) GCFREE(w);
        !           718:       w = (MP *)MALLOC(l*sizeof(MP));
        !           719:       wlen = l;
        !           720:     }
        !           721:     for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
        !           722:       w[i] = m;
        !           723:     for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           724:       _muldm_dup(vl,p1,w[i],&t); _addd_destructive(vl,s,t,&u); s = u;
        !           725:     }
        !           726:     bzero(w,l*sizeof(MP));
        !           727:     *pr = s;
        !           728:   }
1.9       noro      729: }
                    730:
1.11      noro      731: void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      732: {
1.17    ! noro      733:   MP m;
        !           734:   DP s,t,u;
        !           735:   int i,l;
        !           736:   static MP *w;
        !           737:   static int wlen;
        !           738:
        !           739:   if ( !p1 || !p2 )
        !           740:     *pr = 0;
        !           741:   else {
        !           742:     for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
        !           743:     if ( l > wlen ) {
        !           744:       if ( w ) GCFREE(w);
        !           745:       w = (MP *)MALLOC(l*sizeof(MP));
        !           746:       wlen = l;
        !           747:     }
        !           748:     for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
        !           749:       w[i] = m;
        !           750:     for ( s = 0, i = l-1; i >= 0; i-- ) {
        !           751:       _weyl_muldm_dup(vl,w[i],p2,&t); _addd_destructive(vl,s,t,&u); s = u;
        !           752:     }
        !           753:     bzero(w,l*sizeof(MP));
        !           754:     *pr = s;
        !           755:   }
1.9       noro      756: }
                    757:
1.11      noro      758: void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
1.9       noro      759: {
1.17    ! noro      760:   MP m,mr,mr0;
        !           761:   DL d,dt,dm;
        !           762:   P c;
        !           763:   int n,i;
        !           764:   int *pt,*p1,*p2;
        !           765:
        !           766:   if ( !p )
        !           767:     *pr = 0;
        !           768:   else {
        !           769:     for ( mr0 = 0, m = BDY(p), c = (P)C(m0), d = m0->dl, n = NV(p);
        !           770:       m; m = NEXT(m) ) {
        !           771:       _NEXTMP(mr0,mr);
        !           772:       mulp(vl,(P)C(m),c,(P *)&C(mr));
        !           773:       _NEWDL_NOINIT(dt,n); mr->dl = dt;
        !           774:       dm = m->dl;
        !           775:       dt->td = d->td + dm->td;
        !           776:       for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
        !           777:         *pt++ = *p1++ + *p2++;
        !           778:     }
        !           779:     NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
        !           780:     if ( *pr )
        !           781:       (*pr)->sugar = p->sugar + m0->dl->td;
        !           782:   }
1.9       noro      783: }
                    784:
1.11      noro      785: void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
1.9       noro      786: {
1.17    ! noro      787:   DP r,t,t1;
        !           788:   MP m;
        !           789:   DL d0;
        !           790:   int n,n2,l,i,j,tlen;
        !           791:   static MP *w,*psum;
        !           792:   static struct cdl *tab;
        !           793:   static int wlen;
        !           794:   static int rtlen;
        !           795:
        !           796:   if ( !p )
        !           797:     *pr = 0;
        !           798:   else {
        !           799:     for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
        !           800:     if ( l > wlen ) {
        !           801:       if ( w ) GCFREE(w);
        !           802:       w = (MP *)MALLOC(l*sizeof(MP));
        !           803:       wlen = l;
        !           804:     }
        !           805:     for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
        !           806:       w[i] = m;
        !           807:     n = NV(p); n2 = n>>1;
        !           808:     d0 = m0->dl;
        !           809:
        !           810:     for ( i = 0, tlen = 1; i < n2; i++ )
        !           811:       tlen *= d0->d[n2+i]+1;
        !           812:     if ( tlen > rtlen ) {
        !           813:       if ( tab ) GCFREE(tab);
        !           814:       if ( psum ) GCFREE(psum);
        !           815:       rtlen = tlen;
        !           816:       tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
        !           817:       psum = (MP *)MALLOC(rtlen*sizeof(MP));
        !           818:     }
        !           819:     bzero(psum,tlen*sizeof(MP));
        !           820:     for ( i = l-1; i >= 0; i-- ) {
        !           821:       bzero(tab,tlen*sizeof(struct cdl));
        !           822:       _weyl_mulmm_dup(vl,m0,w[i],n,tab,tlen);
        !           823:       for ( j = 0; j < tlen; j++ ) {
        !           824:         if ( tab[j].c ) {
        !           825:           _NEWMP(m); m->dl = tab[j].d;
        !           826:           C(m) = tab[j].c; NEXT(m) = psum[j];
        !           827:           psum[j] = m;
        !           828:         }
        !           829:       }
        !           830:     }
        !           831:     for ( j = tlen-1, r = 0; j >= 0; j-- )
        !           832:       if ( psum[j] ) {
        !           833:         _MKDP(n,psum[j],t); _addd_destructive(vl,r,t,&t1); r = t1;
        !           834:       }
        !           835:     if ( r )
        !           836:       r->sugar = p->sugar + m0->dl->td;
        !           837:     *pr = r;
        !           838:   }
1.9       noro      839: }
                    840:
                    841: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    842:
1.11      noro      843: void _weyl_mulmm_dup(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.9       noro      844: {
1.17    ! noro      845:   P c;
        !           846:   DL d,d0,d1,dt;
        !           847:   int i,j,a,b,k,l,n2,s,min,curlen;
        !           848:   struct cdl *p;
        !           849:   static Q *ctab;
        !           850:   static struct cdl *tab;
        !           851:   static int tablen;
        !           852:   static struct cdl *tmptab;
        !           853:   static int tmptablen;
        !           854:
        !           855:   if ( !m0 || !m1 ) {
        !           856:     rtab[0].c = 0;
        !           857:     rtab[0].d = 0;
        !           858:     return;
        !           859:   }
        !           860:   mulp(vl,(P)C(m0),(P)C(m1),&c);
        !           861:   d0 = m0->dl; d1 = m1->dl;
        !           862:   n2 = n>>1;
        !           863:   curlen = 1;
        !           864:
        !           865:   _NEWDL(d,n);
        !           866:   if ( n & 1 )
        !           867:     /* offset of h-degree */
        !           868:     d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
        !           869:   else
        !           870:     d->td = 0;
        !           871:   rtab[0].c = (Obj)c;
        !           872:   rtab[0].d = d;
        !           873:
        !           874:   if ( rtablen > tmptablen ) {
        !           875:     if ( tmptab ) GCFREE(tmptab);
        !           876:     tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
        !           877:     tmptablen = rtablen;
        !           878:   }
        !           879:
        !           880:   for ( i = 0; i < n2; i++ ) {
        !           881:     a = d0->d[i]; b = d1->d[n2+i];
        !           882:     k = d0->d[n2+i]; l = d1->d[i];
        !           883:
        !           884:     /* degree of xi^a*(Di^k*xi^l)*Di^b */
        !           885:     a += l;
        !           886:     b += k;
        !           887:     s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
        !           888:
        !           889:     if ( !k || !l ) {
        !           890:       for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
        !           891:         if ( p->c ) {
        !           892:           dt = p->d;
        !           893:           dt->d[i] = a;
        !           894:           dt->d[n2+i] = b;
        !           895:           dt->td += s;
        !           896:         }
        !           897:       }
        !           898:       curlen *= k+1;
        !           899:       continue;
        !           900:     }
        !           901:     if ( k+1 > tablen ) {
        !           902:       if ( tab ) GCFREE(tab);
        !           903:       if ( ctab ) GCFREE(ctab);
        !           904:       tablen = k+1;
        !           905:       tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
        !           906:       ctab = (Q *)MALLOC(tablen*sizeof(P));
        !           907:     }
        !           908:     /* compute xi^a*(Di^k*xi^l)*Di^b */
        !           909:     min = MIN(k,l);
        !           910:     mkwc(k,l,ctab);
        !           911:     bzero(tab,(k+1)*sizeof(struct cdl));
        !           912:     /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
        !           913:     if ( n & 1 )
        !           914:       for ( j = 0; j <= min; j++ ) {
        !           915:         _NEWDL(d,n);
        !           916:         d->d[i] = a-j; d->d[n2+i] = b-j;
        !           917:         d->td = s;
        !           918:         d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
        !           919:         tab[j].d = d;
        !           920:         tab[j].c = (Obj)ctab[j];
        !           921:       }
        !           922:     else
        !           923:       for ( j = 0; j <= min; j++ ) {
        !           924:         _NEWDL(d,n);
        !           925:         d->d[i] = a-j; d->d[n2+i] = b-j;
        !           926:         d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
        !           927:         tab[j].d = d;
        !           928:         tab[j].c = (Obj)ctab[j];
        !           929:       }
1.9       noro      930: #if 0
1.17    ! noro      931:     _comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
        !           932:     for ( j = 0; j < curlen; j++ )
        !           933:       if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
        !           934:     for ( j = 0; j <= min; j++ )
        !           935:       if ( tab[j].d ) { _FREEDL(tab[j].d); }
        !           936:     curlen *= k+1;
        !           937:     bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
1.9       noro      938: #else
1.17    ! noro      939:     _comm_muld_tab_destructive(vl,n,rtab,curlen,tab,k+1);
        !           940:     for ( j = 0; j <= min; j++ )
        !           941:       if ( tab[j].d ) { _FREEDL(tab[j].d); }
        !           942:     curlen *= k+1;
1.9       noro      943: #endif
1.17    ! noro      944:   }
1.9       noro      945: }
                    946:
                    947: /* direct product of two cdl tables
                    948:   rt[] = [
                    949:     t[0]*t1[0],...,t[n-1]*t1[0],
                    950:     t[0]*t1[1],...,t[n-1]*t1[1],
                    951:     ...
                    952:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                    953:   ]
                    954: */
                    955:
1.11      noro      956: void _comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.9       noro      957: {
1.17    ! noro      958:   int i,j;
        !           959:   struct cdl *p;
        !           960:   P c;
        !           961:   DL d;
        !           962:
        !           963:   bzero(rt,n*n1*sizeof(struct cdl));
        !           964:   for ( j = 0, p = rt; j < n1; j++ ) {
        !           965:     c = (P)t1[j].c;
        !           966:     d = t1[j].d;
        !           967:     if ( !c )
        !           968:       break;
        !           969:     for ( i = 0; i < n; i++, p++ ) {
        !           970:       if ( t[i].c ) {
        !           971:         mulp(vl,(P)t[i].c,c,(P *)&p->c);
        !           972:         _adddl_dup(nv,t[i].d,d,&p->d);
        !           973:       }
        !           974:     }
        !           975:   }
1.9       noro      976: }
                    977:
1.11      noro      978: void _comm_muld_tab_destructive(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1)
1.9       noro      979: {
1.17    ! noro      980:   int i,j;
        !           981:   struct cdl *p;
        !           982:   P c;
        !           983:   DL d;
        !           984:
        !           985:   for ( j = 1, p = t+n; j < n1; j++ ) {
        !           986:     c = (P)t1[j].c;
        !           987:     d = t1[j].d;
        !           988:     if ( !c )
        !           989:       break;
        !           990:     for ( i = 0; i < n; i++, p++ ) {
        !           991:       if ( t[i].c ) {
        !           992:         mulp(vl,(P)t[i].c,c,(P *)&p->c);
        !           993:         _adddl_dup(nv,t[i].d,d,&p->d);
        !           994:       }
        !           995:     }
        !           996:   }
        !           997:   c = (P)t1[0].c;
        !           998:   d = t1[0].d;
        !           999:   for ( i = 0, p = t; i < n; i++, p++ )
        !          1000:     if ( t[i].c ) {
        !          1001:       mulp(vl,(P)t[i].c,c,(P *)&p->c);
        !          1002:       /* t[i].d += d */
        !          1003:       adddl_destructive(nv,t[i].d,d);
        !          1004:     }
1.1       noro     1005: }

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