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

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.14    ! ohara      48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.13 2002/07/24 07:38:55 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: {
                     65:        _mp_free_list = 0;
                     66:        _dp_free_list = 0;
                     67:        _dl_free_list = 0;
                     68:        GC_gcollect();
                     69: }
1.1       noro       70:
                     71: void _DL_alloc()
                     72: {
                     73:        int *p;
                     74:        int i,dl_len;
                     75:        static int DL_alloc_count;
                     76:
                     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.8       noro       80:        if ( dl_len & 1 )
                     81:                dl_len += 1;
                     82: #endif
1.1       noro       83:        for ( i = 0; i < 128; i++, p += dl_len ) {
1.13      noro       84:                p = (int *)GC_malloc(dl_len*sizeof(int));
1.1       noro       85:                *(DL *)p = _dl_free_list;
                     86:                _dl_free_list = (DL)p;
                     87:        }
                     88: }
                     89:
                     90: void _MP_alloc()
                     91: {
                     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++ ) {
1.13      noro       98:                p = (MP)GC_malloc(sizeof(struct oMP));
                     99:                p->next = _mp_free_list; _mp_free_list = p;
1.1       noro      100:        }
                    101: }
                    102:
                    103: void _DP_alloc()
                    104: {
                    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++ ) {
1.13      noro      111:                p = (DP)GC_malloc(sizeof(struct oDP));
                    112:                p->body = (MP)_dp_free_list; _dp_free_list = p;
1.1       noro      113:        }
                    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: {
                    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) = STOI(t);
                    138:                                        } else {
1.5       noro      139:                                                _FREEDL(s->dl); _FREEMP(s);
1.1       noro      140:                                        }
1.5       noro      141:                                        s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
1.1       noro      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);
1.5       noro      168:                _FREEDP(p1); _FREEDP(p2);
1.1       noro      169:        }
                    170: }
                    171:
1.11      noro      172: void _mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      173: {
                    174:        if ( !do_weyl )
                    175:                _comm_mulmd_dup(mod,p1,p2,pr);
                    176:        else
                    177:                _weyl_mulmd_dup(mod,p1,p2,pr);
                    178: }
                    179:
1.11      noro      180: void _comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      181: {
                    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 ) GC_free(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:        }
                    210: }
                    211:
1.11      noro      212: void _weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr)
1.1       noro      213: {
                    214:        MP m;
                    215:        DP s,t,u;
1.11      noro      216:        int i,l;
1.1       noro      217:        static MP *w;
                    218:        static int wlen;
                    219:
                    220:        if ( !p1 || !p2 )
                    221:                *pr = 0;
                    222:        else {
1.4       noro      223:                for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.1       noro      224:                if ( l > wlen ) {
                    225:                        if ( w ) GC_free(w);
                    226:                        w = (MP *)MALLOC(l*sizeof(MP));
                    227:                        wlen = l;
                    228:                }
1.4       noro      229:                for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.1       noro      230:                        w[i] = m;
                    231:                for ( s = 0, i = l-1; i >= 0; i-- ) {
1.4       noro      232:                        _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
1.1       noro      233:                }
                    234:                bzero(w,l*sizeof(MP));
                    235:                *pr = s;
                    236:        }
                    237: }
                    238:
1.11      noro      239: void _mulmdm_dup(int mod,DP p,MP m0,DP *pr)
1.1       noro      240: {
                    241:        MP m,mr,mr0;
                    242:        DL d,dt,dm;
1.11      noro      243:        int c,n,i,c1,c2;
1.1       noro      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);
1.10      noro      252:                        c1 = ITOS(C(m));
                    253:                        DMAR(c1,c,0,mod,c2);
                    254:                        C(mr) = (P)STOI(c2);
1.1       noro      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:        }
                    265: }
                    266:
1.11      noro      267: void _weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr)
1.1       noro      268: {
                    269:        DP r,t,t1;
                    270:        MP m;
1.4       noro      271:        DL d0;
                    272:        int n,n2,l,i,j,tlen;
                    273:        static MP *w,*psum;
                    274:        static struct cdlm *tab;
1.1       noro      275:        static int wlen;
1.4       noro      276:        static int rtlen;
1.1       noro      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 ) GC_free(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;
1.4       noro      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 ) GC_free(tab);
                    296:                        if ( psum ) GC_free(psum);
                    297:                        rtlen = tlen;
                    298:                        tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
                    299:                        psum = (MP *)MALLOC(rtlen*sizeof(MP));
1.1       noro      300:                }
1.4       noro      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) = 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:                        }
1.1       noro      317:                if ( r )
                    318:                        r->sugar = p->sugar + m0->dl->td;
                    319:                *pr = r;
                    320:        }
                    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.11      noro      327:        int c,c0,c1;
1.4       noro      328:        DL d,d0,d1,dt;
1.11      noro      329:        int i,j,a,b,k,l,n2,s,min,curlen;
1.4       noro      330:        struct cdlm *p;
                    331:        static int *ctab;
                    332:        static struct cdlm *tab;
1.1       noro      333:        static int tablen;
1.4       noro      334:        static struct cdlm *tmptab;
                    335:        static int tmptablen;
1.1       noro      336:
1.4       noro      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 ) GC_free(tmptab);
                    359:                tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
                    360:                tmptablen = rtablen;
                    361:        }
1.1       noro      362:
1.4       noro      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];
1.12      noro      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:
1.4       noro      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 ) GC_free(tab);
                    386:                        if ( ctab ) GC_free(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 */
1.1       noro      396:                if ( n & 1 )
1.4       noro      397:                        for ( j = 0; j <= min; j++ ) {
                    398:                                _NEWDL(d,n);
1.12      noro      399:                                d->d[i] = a-j; d->d[n2+i] = b-j;
1.4       noro      400:                                d->td = s;
1.12      noro      401:                                d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.4       noro      402:                                tab[j].d = d;
                    403:                                tab[j].c = ctab[j];
                    404:                        }
1.1       noro      405:                else
1.4       noro      406:                        for ( j = 0; j <= min; j++ ) {
                    407:                                _NEWDL(d,n);
1.12      noro      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 */
1.4       noro      410:                                tab[j].d = d;
                    411:                                tab[j].c = ctab[j];
                    412:                        }
                    413: #if 0
                    414:                _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
                    415:                for ( j = 0; j < curlen; j++ )
1.5       noro      416:                        if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1.4       noro      417:                for ( j = 0; j <= min; j++ )
1.5       noro      418:                        if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4       noro      419:                curlen *= k+1;
                    420:                bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
                    421: #else
                    422:                _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
                    423:                for ( j = 0; j <= min; j++ )
1.5       noro      424:                        if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4       noro      425:                curlen *= k+1;
                    426: #endif
                    427:        }
                    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: {
                    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);
1.5       noro      455:                                _adddl_dup(nv,t[i].d,d,&p->d);
1.4       noro      456:                        }
                    457:                }
                    458:        }
                    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: {
                    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);
1.5       noro      476:                                _adddl_dup(nv,t[i].d,d,&p->d);
1.1       noro      477:                        }
                    478:                }
                    479:        }
1.4       noro      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: {
                    492:        int i,n;
1.5       noro      493:        DL t;
1.1       noro      494:
1.5       noro      495:        n = current_dl_length;
                    496:        _NEWDL(t,n); *dr = t;
                    497:        t->td = d->td;
1.1       noro      498:        for ( i = 0; i < n; i++ )
1.5       noro      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: {
                    504:        int i,n;
                    505:        DL t;
                    506:
                    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];
                    512: }
                    513:
1.11      noro      514: void _adddl_dup(int n,DL d1,DL d2,DL *dr)
1.1       noro      515: {
                    516:        DL dt;
                    517:        int i;
                    518:
1.4       noro      519:        _NEWDL(dt,n);
                    520:        *dr = dt;
1.1       noro      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: {
                    528:        int i;
                    529:
1.5       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: {
                    535:        MP m,s;
                    536:
                    537:        if ( !f )
                    538:                return;
                    539:        m = f->body;
                    540:        while ( m ) {
1.5       noro      541:                s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
                    542:        }
                    543:        _FREEDP(f);
                    544: }
                    545:
1.11      noro      546: void dpto_dp(DP p,DP *r)
1.5       noro      547: {
                    548:        MP m,mr0,mr;
1.7       noro      549:        DL t;
1.5       noro      550:
                    551:        if ( !p )
                    552:                *r = 0;
                    553:        else {
1.7       noro      554:                /* XXX : dummy call to set current_dl_length */
                    555:                _NEWDL_NOINIT(t,NV(p));
                    556:
1.5       noro      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;
1.1       noro      565:        }
                    566: }
                    567:
1.11      noro      568: void _dptodp(DP p,DP *r)
1.1       noro      569: {
                    570:        MP m,mr0,mr;
                    571:
                    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:        }
                    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: {
                    595:        NODE top,prev,cur,m,t;
                    596:
                    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;
1.1       noro      611:                                break;
                    612:                }
1.5       noro      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;
1.1       noro      629:                        }
                    630:                }
1.5       noro      631:                if ( !cur )
                    632:                        NEXT(prev) = m;
                    633:                return top;
1.1       noro      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: {
                    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,C(m1),C(m2),&t);
                    654:                                        s = m1; m1 = NEXT(m1);
                    655:                                        if ( t ) {
                    656:                                                _NEXTMP2(mr0,mr,s); C(mr) = 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:        }
                    689: }
                    690:
1.11      noro      691: void _muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      692: {
                    693:        if ( !do_weyl )
                    694:                _comm_muld_dup(vl,p1,p2,pr);
                    695:        else
                    696:                _weyl_muld_dup(vl,p1,p2,pr);
                    697: }
                    698:
1.11      noro      699: void _comm_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      700: {
                    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 ) GC_free(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:        }
                    729: }
                    730:
1.11      noro      731: void _weyl_muld_dup(VL vl,DP p1,DP p2,DP *pr)
1.9       noro      732: {
                    733:        MP m;
                    734:        DP s,t,u;
1.11      noro      735:        int i,l;
1.9       noro      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 ) GC_free(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:        }
                    756: }
                    757:
1.11      noro      758: void _muldm_dup(VL vl,DP p,MP m0,DP *pr)
1.9       noro      759: {
                    760:        MP m,mr,mr0;
                    761:        DL d,dt,dm;
                    762:        P c;
1.11      noro      763:        int n,i;
1.9       noro      764:        int *pt,*p1,*p2;
                    765:
                    766:        if ( !p )
                    767:                *pr = 0;
                    768:        else {
                    769:                for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
                    770:                        m; m = NEXT(m) ) {
                    771:                        _NEXTMP(mr0,mr);
                    772:                        mulp(vl,C(m),c,&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:        }
                    783: }
                    784:
1.11      noro      785: void _weyl_muldm_dup(VL vl,MP m0,DP p,DP *pr)
1.9       noro      786: {
                    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 ) GC_free(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 ) GC_free(tab);
                    814:                        if ( psum ) GC_free(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:        }
                    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: {
                    845:        P c;
                    846:        DL d,d0,d1,dt;
1.11      noro      847:        int i,j,a,b,k,l,n2,s,min,curlen;
1.9       noro      848:        struct cdl *p;
1.11      noro      849:        static Q *ctab;
1.9       noro      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,C(m0),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 = c;
                    872:        rtab[0].d = d;
                    873:
                    874:        if ( rtablen > tmptablen ) {
                    875:                if ( tmptab ) GC_free(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];
1.12      noro      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:
1.9       noro      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 ) GC_free(tab);
                    903:                        if ( ctab ) GC_free(ctab);
                    904:                        tablen = k+1;
                    905:                        tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1.11      noro      906:                        ctab = (Q *)MALLOC(tablen*sizeof(P));
1.9       noro      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);
1.12      noro      916:                                d->d[i] = a-j; d->d[n2+i] = b-j;
1.9       noro      917:                                d->td = s;
1.12      noro      918:                                d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.9       noro      919:                                tab[j].d = d;
1.11      noro      920:                                tab[j].c = (P)ctab[j];
1.9       noro      921:                        }
                    922:                else
                    923:                        for ( j = 0; j <= min; j++ ) {
                    924:                                _NEWDL(d,n);
1.12      noro      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 */
1.9       noro      927:                                tab[j].d = d;
1.11      noro      928:                                tab[j].c = (P)ctab[j];
1.9       noro      929:                        }
                    930: #if 0
                    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));
                    938: #else
                    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;
                    943: #endif
                    944:        }
                    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: {
                    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 = 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,t[i].c,c,&p->c);
                    972:                                _adddl_dup(nv,t[i].d,d,&p->d);
                    973:                        }
                    974:                }
                    975:        }
                    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: {
                    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 = 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,t[i].c,c,&p->c);
                    993:                                _adddl_dup(nv,t[i].d,d,&p->d);
                    994:                        }
                    995:                }
                    996:        }
                    997:        c = 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,t[i].c,c,&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>