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

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.5     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.4 2000/11/07 06:06:39 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:
1.5     ! noro       56: void dpto_dp();
1.1       noro       57: void _dptodp();
1.5     ! noro       58: void _adddl_dup();
1.4       noro       59: void adddl_destructive();
1.1       noro       60: void _mulmdm_dup();
                     61: void _free_dp();
                     62: void _comm_mulmd_dup();
                     63: void _weyl_mulmd_dup();
                     64: void _weyl_mulmmm_dup();
                     65: void _weyl_mulmdm_dup();
1.4       noro       66: void _comm_mulmd_tab();
                     67: void _comm_mulmd_tab_destructive();
1.1       noro       68:
                     69: MP _mp_free_list;
                     70: DP _dp_free_list;
                     71: DL _dl_free_list;
                     72: int current_dl_length;
                     73:
                     74: void _DL_alloc()
                     75: {
                     76:        int *p;
                     77:        int i,dl_len;
                     78:        static int DL_alloc_count;
                     79:
                     80: /*     fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
                     81:        dl_len = (current_dl_length+1);
                     82:        p = (int *)GC_malloc(128*dl_len*sizeof(int));
                     83:        for ( i = 0; i < 128; i++, p += dl_len ) {
                     84:                *(DL *)p = _dl_free_list;
                     85:                _dl_free_list = (DL)p;
                     86:        }
                     87: }
                     88:
                     89: void _MP_alloc()
                     90: {
                     91:        MP p;
                     92:        int i;
                     93:        static int MP_alloc_count;
                     94:
                     95: /*     fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */
                     96:        p = (MP)GC_malloc(1024*sizeof(struct oMP));
                     97:        for ( i = 0; i < 1024; i++ ) {
                     98:                p[i].next = _mp_free_list; _mp_free_list = &p[i];
                     99:        }
                    100: }
                    101:
                    102: void _DP_alloc()
                    103: {
                    104:        DP p;
                    105:        int i;
                    106:        static int DP_alloc_count;
                    107:
                    108: /*     fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */
                    109:        p = (DP)GC_malloc(1024*sizeof(struct oDP));
                    110:        for ( i = 0; i < 1024; i++ ) {
                    111:                p[i].body = (MP)_dp_free_list; _dp_free_list = &p[i];
                    112:        }
                    113: }
                    114:
                    115: /* merge p1 and p2 into pr */
                    116:
                    117: void _addmd_destructive(mod,p1,p2,pr)
                    118: int mod;
                    119: DP p1,p2,*pr;
                    120: {
                    121:        int n;
                    122:        MP m1,m2,mr,mr0,s;
                    123:        int t;
                    124:
                    125:        if ( !p1 )
                    126:                *pr = p2;
                    127:        else if ( !p2 )
                    128:                *pr = p1;
                    129:        else {
                    130:                for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
                    131:                        switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
                    132:                                case 0:
                    133:                                        t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
                    134:                                        if ( t < 0 )
                    135:                                                t += mod;
                    136:                                        s = m1; m1 = NEXT(m1);
                    137:                                        if ( t ) {
                    138:                                                _NEXTMP2(mr0,mr,s); C(mr) = STOI(t);
                    139:                                        } else {
1.5     ! noro      140:                                                _FREEDL(s->dl); _FREEMP(s);
1.1       noro      141:                                        }
1.5     ! noro      142:                                        s = m2; m2 = NEXT(m2); _FREEDL(s->dl); _FREEMP(s);
1.1       noro      143:                                        break;
                    144:                                case 1:
                    145:                                        s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
                    146:                                        break;
                    147:                                case -1:
                    148:                                        s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
                    149:                                        break;
                    150:                        }
                    151:                if ( !mr0 )
                    152:                        if ( m1 )
                    153:                                mr0 = m1;
                    154:                        else if ( m2 )
                    155:                                mr0 = m2;
                    156:                        else {
                    157:                                *pr = 0;
                    158:                                return;
                    159:                        }
                    160:                else if ( m1 )
                    161:                        NEXT(mr) = m1;
                    162:                else if ( m2 )
                    163:                        NEXT(mr) = m2;
                    164:                else
                    165:                        NEXT(mr) = 0;
                    166:                _MKDP(NV(p1),mr0,*pr);
                    167:                if ( *pr )
                    168:                        (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.5     ! noro      169:                _FREEDP(p1); _FREEDP(p2);
1.1       noro      170:        }
                    171: }
                    172:
                    173: void _mulmd_dup(mod,p1,p2,pr)
                    174: int mod;
                    175: DP p1,p2,*pr;
                    176: {
                    177:        if ( !do_weyl )
                    178:                _comm_mulmd_dup(mod,p1,p2,pr);
                    179:        else
                    180:                _weyl_mulmd_dup(mod,p1,p2,pr);
                    181: }
                    182:
                    183: void _comm_mulmd_dup(mod,p1,p2,pr)
                    184: int mod;
                    185: DP p1,p2,*pr;
                    186: {
                    187:        MP m;
                    188:        DP s,t,u;
                    189:        int i,l,l1;
                    190:        static MP *w;
                    191:        static int wlen;
                    192:
                    193:        if ( !p1 || !p2 )
                    194:                *pr = 0;
                    195:        else {
                    196:                for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
                    197:                for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
                    198:                if ( l1 < l ) {
                    199:                        t = p1; p1 = p2; p2 = t;
                    200:                        l = l1;
                    201:                }
                    202:                if ( l > wlen ) {
                    203:                        if ( w ) GC_free(w);
                    204:                        w = (MP *)MALLOC(l*sizeof(MP));
                    205:                        wlen = l;
                    206:                }
                    207:                for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
                    208:                        w[i] = m;
                    209:                for ( s = 0, i = l-1; i >= 0; i-- ) {
                    210:                        _mulmdm_dup(mod,p1,w[i],&t); _addmd_destructive(mod,s,t,&u); s = u;
                    211:                }
                    212:                bzero(w,l*sizeof(MP));
                    213:                *pr = s;
                    214:        }
                    215: }
                    216:
                    217: void _weyl_mulmd_dup(mod,p1,p2,pr)
                    218: int mod;
                    219: DP p1,p2,*pr;
                    220: {
                    221:        MP m;
                    222:        DP s,t,u;
                    223:        int i,l,l1;
                    224:        static MP *w;
                    225:        static int wlen;
                    226:
                    227:        if ( !p1 || !p2 )
                    228:                *pr = 0;
                    229:        else {
1.4       noro      230:                for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.1       noro      231:                if ( l > wlen ) {
                    232:                        if ( w ) GC_free(w);
                    233:                        w = (MP *)MALLOC(l*sizeof(MP));
                    234:                        wlen = l;
                    235:                }
1.4       noro      236:                for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.1       noro      237:                        w[i] = m;
                    238:                for ( s = 0, i = l-1; i >= 0; i-- ) {
1.4       noro      239:                        _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
1.1       noro      240:                }
                    241:                bzero(w,l*sizeof(MP));
                    242:                *pr = s;
                    243:        }
                    244: }
                    245:
                    246: void _mulmdm_dup(mod,p,m0,pr)
                    247: int mod;
                    248: DP p;
                    249: MP m0;
                    250: DP *pr;
                    251: {
                    252:        MP m,mr,mr0;
                    253:        DL d,dt,dm;
                    254:        int c,n,r,i;
                    255:        int *pt,*p1,*p2;
                    256:
                    257:        if ( !p )
                    258:                *pr = 0;
                    259:        else {
                    260:                for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
                    261:                        m; m = NEXT(m) ) {
                    262:                        _NEXTMP(mr0,mr);
                    263:                        C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
                    264:                        _NEWDL_NOINIT(dt,n); mr->dl = dt;
                    265:                        dm = m->dl;
                    266:                        dt->td = d->td + dm->td;
                    267:                        for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
                    268:                                *pt++ = *p1++ + *p2++;
                    269:                }
                    270:                NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
                    271:                if ( *pr )
                    272:                        (*pr)->sugar = p->sugar + m0->dl->td;
                    273:        }
                    274: }
                    275:
1.4       noro      276: void _weyl_mulmmm_dup_bug();
                    277:
                    278: void _weyl_mulmdm_dup(mod,m0,p,pr)
1.1       noro      279: int mod;
1.4       noro      280: MP m0;
1.1       noro      281: DP p;
                    282: DP *pr;
                    283: {
                    284:        DP r,t,t1;
                    285:        MP m;
1.4       noro      286:        DL d0;
                    287:        int n,n2,l,i,j,tlen;
                    288:        static MP *w,*psum;
                    289:        static struct cdlm *tab;
1.1       noro      290:        static int wlen;
1.4       noro      291:        static int rtlen;
1.1       noro      292:
                    293:        if ( !p )
                    294:                *pr = 0;
                    295:        else {
                    296:                for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
                    297:                if ( l > wlen ) {
                    298:                        if ( w ) GC_free(w);
                    299:                        w = (MP *)MALLOC(l*sizeof(MP));
                    300:                        wlen = l;
                    301:                }
                    302:                for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
                    303:                        w[i] = m;
1.4       noro      304:                n = NV(p); n2 = n>>1;
                    305:                d0 = m0->dl;
                    306:
                    307:                for ( i = 0, tlen = 1; i < n2; i++ )
                    308:                        tlen *= d0->d[n2+i]+1;
                    309:                if ( tlen > rtlen ) {
                    310:                        if ( tab ) GC_free(tab);
                    311:                        if ( psum ) GC_free(psum);
                    312:                        rtlen = tlen;
                    313:                        tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
                    314:                        psum = (MP *)MALLOC(rtlen*sizeof(MP));
1.1       noro      315:                }
1.4       noro      316:                bzero(psum,tlen*sizeof(MP));
                    317:                for ( i = l-1; i >= 0; i-- ) {
                    318:                        bzero(tab,tlen*sizeof(struct cdlm));
                    319:                        _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
                    320:                        for ( j = 0; j < tlen; j++ ) {
                    321:                                if ( tab[j].c ) {
                    322:                                        _NEWMP(m); m->dl = tab[j].d;
                    323:                                        C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
                    324:                                        psum[j] = m;
                    325:                                }
                    326:                        }
                    327:                }
                    328:                for ( j = tlen-1, r = 0; j >= 0; j-- )
                    329:                        if ( psum[j] ) {
                    330:                                _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
                    331:                        }
1.1       noro      332:                if ( r )
                    333:                        r->sugar = p->sugar + m0->dl->td;
                    334:                *pr = r;
                    335:        }
                    336: }
1.4       noro      337:
1.1       noro      338: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
                    339:
1.4       noro      340: void _weyl_mulmmm_dup(mod,m0,m1,n,rtab,rtablen)
1.1       noro      341: int mod;
                    342: MP m0,m1;
                    343: int n;
1.4       noro      344: struct cdlm *rtab;
                    345: int rtablen;
1.1       noro      346: {
                    347:        MP m,mr,mr0;
                    348:        DP r,t,t1;
                    349:        int c,c0,c1,cc;
1.4       noro      350:        DL d,d0,d1,dt;
                    351:        int i,j,a,b,k,l,n2,s,min,h,curlen;
                    352:        struct cdlm *p;
                    353:        static int *ctab;
                    354:        static struct cdlm *tab;
1.1       noro      355:        static int tablen;
1.4       noro      356:        static struct cdlm *tmptab;
                    357:        static int tmptablen;
1.1       noro      358:
1.4       noro      359:        if ( !m0 || !m1 ) {
                    360:                rtab[0].c = 0;
                    361:                rtab[0].d = 0;
                    362:                return;
                    363:        }
                    364:        c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
                    365:        c = dmar(c0,c1,0,mod);
                    366:        d0 = m0->dl; d1 = m1->dl;
                    367:        n2 = n>>1;
                    368:        curlen = 1;
                    369:
                    370:        _NEWDL(d,n);
                    371:        if ( n & 1 )
                    372:                /* offset of h-degree */
                    373:                d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
                    374:        else
                    375:                d->td = 0;
                    376:        rtab[0].c = c;
                    377:        rtab[0].d = d;
                    378:
                    379:        if ( rtablen > tmptablen ) {
                    380:                if ( tmptab ) GC_free(tmptab);
                    381:                tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
                    382:                tmptablen = rtablen;
                    383:        }
1.1       noro      384:
1.4       noro      385:        for ( i = 0; i < n2; i++ ) {
                    386:                a = d0->d[i]; b = d1->d[n2+i];
                    387:                k = d0->d[n2+i]; l = d1->d[i];
                    388:                if ( !k || !l ) {
                    389:                        a += l;
                    390:                        b += k;
                    391:                        s = a+b;
                    392:                        for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
                    393:                                if ( p->c ) {
                    394:                                        dt = p->d;
                    395:                                        dt->d[i] = a;
                    396:                                        dt->d[n2+i] = b;
                    397:                                        dt->td += s;
                    398:                                }
                    399:                        }
                    400:                        curlen *= k+1;
                    401:                        continue;
                    402:                }
                    403:                if ( k+1 > tablen ) {
                    404:                        if ( tab ) GC_free(tab);
                    405:                        if ( ctab ) GC_free(ctab);
                    406:                        tablen = k+1;
                    407:                        tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
                    408:                        ctab = (int *)MALLOC(tablen*sizeof(int));
                    409:                }
                    410:                /* degree of xi^a*(Di^k*xi^l)*Di^b */
                    411:                s = a+k+l+b;
                    412:                /* compute xi^a*(Di^k*xi^l)*Di^b */
                    413:                min = MIN(k,l);
                    414:                mkwcm(k,l,mod,ctab);
                    415:                bzero(tab,(k+1)*sizeof(struct cdlm));
                    416:                /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1.1       noro      417:                if ( n & 1 )
1.4       noro      418:                        for ( j = 0; j <= min; j++ ) {
                    419:                                _NEWDL(d,n);
                    420:                                d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    421:                                d->td = s;
                    422:                                d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
                    423:                                tab[j].d = d;
                    424:                                tab[j].c = ctab[j];
                    425:                        }
1.1       noro      426:                else
1.4       noro      427:                        for ( j = 0; j <= min; j++ ) {
                    428:                                _NEWDL(d,n);
                    429:                                d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
                    430:                                d->td = d->d[i]+d->d[n2+i]; /* XXX */
                    431:                                tab[j].d = d;
                    432:                                tab[j].c = ctab[j];
                    433:                        }
                    434: #if 0
                    435:                _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
                    436:                for ( j = 0; j < curlen; j++ )
1.5     ! noro      437:                        if ( rtab[j].d ) { _FREEDL(rtab[j].d); }
1.4       noro      438:                for ( j = 0; j <= min; j++ )
1.5     ! noro      439:                        if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4       noro      440:                curlen *= k+1;
                    441:                bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
                    442: #else
                    443:                _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
                    444:                for ( j = 0; j <= min; j++ )
1.5     ! noro      445:                        if ( tab[j].d ) { _FREEDL(tab[j].d); }
1.4       noro      446:                curlen *= k+1;
                    447: #endif
                    448:        }
                    449: }
                    450:
                    451: /* direct product of two cdlm tables
                    452:   rt[] = [
                    453:     t[0]*t1[0],...,t[n-1]*t1[0],
                    454:     t[0]*t1[1],...,t[n-1]*t1[1],
                    455:     ...
                    456:     t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
                    457:   ]
                    458: */
                    459:
                    460: void _comm_mulmd_tab(mod,nv,t,n,t1,n1,rt)
                    461: int mod;
                    462: int nv;
                    463: struct cdlm *t;
                    464: int n;
                    465: struct cdlm *t1;
                    466: int n1;
                    467: struct cdlm *rt;
                    468: {
                    469:        int i,j;
                    470:        struct cdlm *p;
                    471:        int c;
                    472:        DL d;
                    473:
                    474:        bzero(rt,n*n1*sizeof(struct cdlm));
                    475:        for ( j = 0, p = rt; j < n1; j++ ) {
                    476:                c = t1[j].c;
                    477:                d = t1[j].d;
                    478:                if ( !c )
                    479:                        break;
                    480:                for ( i = 0; i < n; i++, p++ ) {
                    481:                        if ( t[i].c ) {
                    482:                                p->c = dmar(t[i].c,c,0,mod);
1.5     ! noro      483:                                _adddl_dup(nv,t[i].d,d,&p->d);
1.4       noro      484:                        }
                    485:                }
                    486:        }
                    487: }
                    488:
                    489: void _comm_mulmd_tab_destructive(mod,nv,t,n,t1,n1)
                    490: int mod;
                    491: int nv;
                    492: struct cdlm *t;
                    493: int n;
                    494: struct cdlm *t1;
                    495: int n1;
                    496: {
                    497:        int i,j;
                    498:        struct cdlm *p;
                    499:        int c;
                    500:        DL d;
                    501:
                    502:        for ( j = 1, p = t+n; j < n1; j++ ) {
                    503:                c = t1[j].c;
                    504:                d = t1[j].d;
                    505:                if ( !c )
                    506:                        break;
                    507:                for ( i = 0; i < n; i++, p++ ) {
                    508:                        if ( t[i].c ) {
                    509:                                p->c = dmar(t[i].c,c,0,mod);
1.5     ! noro      510:                                _adddl_dup(nv,t[i].d,d,&p->d);
1.1       noro      511:                        }
                    512:                }
                    513:        }
1.4       noro      514:        c = t1[0].c;
                    515:        d = t1[0].d;
                    516:        for ( i = 0, p = t; i < n; i++, p++ )
                    517:                if ( t[i].c ) {
                    518:                        p->c = dmar(t[i].c,c,0,mod);
                    519:                        /* t[i].d += d */
                    520:                        adddl_destructive(nv,t[i].d,d);
                    521:                }
1.1       noro      522: }
                    523:
1.5     ! noro      524: void dlto_dl(d,dr)
        !           525: DL d;
        !           526: DL *dr;
1.1       noro      527: {
                    528:        int i,n;
1.5     ! noro      529:        DL t;
1.1       noro      530:
1.5     ! noro      531:        n = current_dl_length;
        !           532:        _NEWDL(t,n); *dr = t;
        !           533:        t->td = d->td;
1.1       noro      534:        for ( i = 0; i < n; i++ )
1.5     ! noro      535:                t->d[i] = d->d[i];
1.1       noro      536: }
                    537:
                    538: void _dltodl(d,dr)
                    539: DL d;
                    540: DL *dr;
                    541: {
                    542:        int i,n;
                    543:        DL t;
                    544:
                    545:        n = current_dl_length;
                    546:        NEWDL(t,n); *dr = t;
                    547:        t->td = d->td;
                    548:        for ( i = 0; i < n; i++ )
                    549:                t->d[i] = d->d[i];
                    550: }
                    551:
1.5     ! noro      552: void _adddl_dup(n,d1,d2,dr)
1.1       noro      553: int n;
                    554: DL d1,d2;
                    555: DL *dr;
                    556: {
                    557:        DL dt;
                    558:        int i;
                    559:
1.4       noro      560:        _NEWDL(dt,n);
                    561:        *dr = dt;
1.1       noro      562:        dt->td = d1->td + d2->td;
                    563:        for ( i = 0; i < n; i++ )
                    564:                dt->d[i] = d1->d[i]+d2->d[i];
1.4       noro      565: }
                    566:
                    567: void _free_dlarray(a,n)
                    568: DL *a;
                    569: int n;
                    570: {
                    571:        int i;
                    572:
1.5     ! noro      573:        for ( i = 0; i < n; i++ ) { _FREEDL(a[i]); }
1.1       noro      574: }
                    575:
                    576: void _free_dp(f)
                    577: DP f;
                    578: {
                    579:        MP m,s;
                    580:
                    581:        if ( !f )
                    582:                return;
                    583:        m = f->body;
                    584:        while ( m ) {
1.5     ! noro      585:                s = m; m = NEXT(m); _FREEDL(s->dl); _FREEMP(s);
        !           586:        }
        !           587:        _FREEDP(f);
        !           588: }
        !           589:
        !           590: void dpto_dp(p,r)
        !           591: DP p;
        !           592: DP *r;
        !           593: {
        !           594:        MP m,mr0,mr;
        !           595:
        !           596:        if ( !p )
        !           597:                *r = 0;
        !           598:        else {
        !           599:                current_dl_length = NV(p);
        !           600:                for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
        !           601:                        _NEXTMP(mr0,mr);
        !           602:                        dlto_dl(m->dl,&mr->dl);
        !           603:                        mr->c = m->c;
        !           604:                }
        !           605:                NEXT(mr) = 0;
        !           606:                _MKDP(p->nv,mr0,*r);
        !           607:                (*r)->sugar = p->sugar;
1.1       noro      608:        }
                    609: }
                    610:
                    611: void _dptodp(p,r)
                    612: DP p;
                    613: DP *r;
                    614: {
                    615:        MP m,mr0,mr;
                    616:
                    617:        if ( !p )
                    618:                *r = 0;
                    619:        else {
                    620:                for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
                    621:                        NEXTMP(mr0,mr);
                    622:                        _dltodl(m->dl,&mr->dl);
                    623:                        mr->c = m->c;
                    624:                }
                    625:                NEXT(mr) = 0;
                    626:                MKDP(p->nv,mr0,*r);
                    627:                (*r)->sugar = p->sugar;
                    628:        }
                    629: }
                    630:
1.5     ! noro      631: /*
        !           632:  * destructive merge of two list
        !           633:  *
        !           634:  * p1, p2 : list of DL
        !           635:  * return : a merged list
        !           636:  */
        !           637:
        !           638: NODE _symb_merge(m1,m2,n)
        !           639: NODE m1,m2;
        !           640: int n;
        !           641: {
        !           642:        NODE top,prev,cur,m,t;
        !           643:
        !           644:        if ( !m1 )
        !           645:                return m2;
        !           646:        else if ( !m2 )
        !           647:                return m1;
        !           648:        else {
        !           649:                switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
        !           650:                        case 0:
        !           651:                                top = m1; _FREEDL((DL)BDY(m2)); m = NEXT(m2);
        !           652:                                break;
        !           653:                        case 1:
        !           654:                                top = m1; m = m2;
        !           655:                                break;
        !           656:                        case -1:
        !           657:                                top = m2; m = m1;
1.1       noro      658:                                break;
                    659:                }
1.5     ! noro      660:                prev = top; cur = NEXT(top);
        !           661:                /* BDY(prev) > BDY(m) always holds */
        !           662:                while ( cur && m ) {
        !           663:                        switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
        !           664:                                case 0:
        !           665:                                        _FREEDL(BDY(m)); m = NEXT(m);
        !           666:                                        prev = cur; cur = NEXT(cur);
        !           667:                                        break;
        !           668:                                case 1:
        !           669:                                        t = NEXT(cur); NEXT(cur) = m; m = t;
        !           670:                                        prev = cur; cur = NEXT(cur);
        !           671:                                        break;
        !           672:                                case -1:
        !           673:                                        NEXT(prev) = m; m = cur;
        !           674:                                        prev = NEXT(prev); cur = NEXT(prev);
        !           675:                                        break;
1.1       noro      676:                        }
                    677:                }
1.5     ! noro      678:                if ( !cur )
        !           679:                        NEXT(prev) = m;
        !           680:                return top;
1.1       noro      681:        }
                    682: }

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