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

Annotation of OpenXM_contrib2/asir2000/builtin/gr.c, Revision 1.54

1.8       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.9       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.8       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.54    ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.53 2003/07/18 10:13:12 noro Exp $
1.8       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "parse.h"
                     52: #include "base.h"
                     53: #include "ox.h"
                     54:
1.27      noro       55: #if defined(__GNUC__)
                     56: #define INLINE inline
                     57: #elif defined(VISUAL)
                     58: #define INLINE __inline
                     59: #else
                     60: #define INLINE
                     61: #endif
                     62:
1.1       noro       63: #define HMAG(p) (p_mag(BDY(p)->c))
                     64:
1.37      noro       65: #define NEWDP_pairs ((DP_pairs)MALLOC(sizeof(struct dp_pairs)))
1.1       noro       66:
1.37      noro       67: static DP_pairs collect_pairs_of_hdlcm( DP_pairs d1, DP_pairs *prest );
                     68: double get_rtime();
1.1       noro       69:
1.14      noro       70: struct oEGT eg_nf,eg_nfm;
                     71: struct oEGT eg_znfm,eg_pz,eg_np,eg_ra,eg_mc,eg_gc;
1.1       noro       72: int TP,NBP,NMP,NFP,NDP,ZR,NZR;
                     73:
                     74: extern int (*cmpdl)();
1.5       noro       75: extern int do_weyl;
1.1       noro       76:
1.13      noro       77: extern DP_Print;
                     78:
1.1       noro       79: extern int dp_nelim;
                     80: extern int dp_fcoeffs;
                     81: static DP *ps,*psm;
                     82: static DL *psh;
                     83: static P *psc;
                     84:
                     85: static int *pss;
                     86: static int psn,pslen;
1.16      noro       87: static int NVars,CNVars;
1.1       noro       88: static VL VC;
1.14      noro       89:
1.16      noro       90: int PCoeffs;
1.14      noro       91: int DP_Print = 0;
1.46      noro       92: int DP_PrintShort = 0;
1.14      noro       93: int DP_Multiple = 0;
1.16      noro       94: int DP_NFStat = 0;
1.14      noro       95: LIST Dist = 0;
                     96: int NoGCD = 0;
                     97: int GenTrace = 0;
                     98: int OXCheck = -1;
                     99:
1.54    ! noro      100: int NoSugar = 0;
1.1       noro      101: static int NoCriB = 0;
                    102: static int NoGC = 0;
                    103: static int NoMC = 0;
                    104: static int NoRA = 0;
                    105: static int ShowMag = 0;
                    106: static int Stat = 0;
                    107: static int Denominator = 1;
1.53      noro      108: int Top = 0;
                    109: int Reverse = 0;
1.1       noro      110: static int Max_mag = 0;
1.42      noro      111: static int Max_coef = 0;
1.54    ! noro      112: char *Demand = 0;
1.1       noro      113: static int PtozpRA = 0;
1.14      noro      114:
1.1       noro      115: int doing_f4;
1.7       noro      116: NODE TraceList;
1.23      noro      117: NODE AllTraceList;
1.1       noro      118:
1.37      noro      119: void Pox_cmo_rpc(NODE,Obj *);
                    120: void Pox_rpc(NODE,Obj *);
                    121: void Pox_pop_local(NODE,Obj *);
                    122:
                    123: INLINE int eqdl(int nv,DL dl1,DL dl2)
1.1       noro      124: {
                    125:        int i;
1.51      noro      126:        int *p1,*p2;
1.1       noro      127:
                    128:        if ( dl1->td != dl2->td )
                    129:                return 0;
1.51      noro      130:        i = nv-1;
                    131:        p1 = dl1->d;
                    132:        p2 = dl2->d;
                    133:        while ( i >= 7 ) {
                    134:                if ( *p1++ != *p2++ ) return 0;
                    135:                if ( *p1++ != *p2++ ) return 0;
                    136:                if ( *p1++ != *p2++ ) return 0;
                    137:                if ( *p1++ != *p2++ ) return 0;
                    138:                if ( *p1++ != *p2++ ) return 0;
                    139:                if ( *p1++ != *p2++ ) return 0;
                    140:                if ( *p1++ != *p2++ ) return 0;
                    141:                if ( *p1++ != *p2++ ) return 0;
                    142:                i -= 8;
                    143:        }
                    144:        switch ( i ) {
                    145:                case 6:
                    146:                        if ( *p1++ != *p2++ ) return 0;
                    147:                        if ( *p1++ != *p2++ ) return 0;
                    148:                        if ( *p1++ != *p2++ ) return 0;
                    149:                        if ( *p1++ != *p2++ ) return 0;
                    150:                        if ( *p1++ != *p2++ ) return 0;
                    151:                        if ( *p1++ != *p2++ ) return 0;
                    152:                        if ( *p1++ != *p2++ ) return 0;
                    153:                        return 1;
                    154:                case 5:
                    155:                        if ( *p1++ != *p2++ ) return 0;
                    156:                        if ( *p1++ != *p2++ ) return 0;
                    157:                        if ( *p1++ != *p2++ ) return 0;
                    158:                        if ( *p1++ != *p2++ ) return 0;
                    159:                        if ( *p1++ != *p2++ ) return 0;
                    160:                        if ( *p1++ != *p2++ ) return 0;
                    161:                        return 1;
                    162:                case 4:
                    163:                        if ( *p1++ != *p2++ ) return 0;
                    164:                        if ( *p1++ != *p2++ ) return 0;
                    165:                        if ( *p1++ != *p2++ ) return 0;
                    166:                        if ( *p1++ != *p2++ ) return 0;
                    167:                        if ( *p1++ != *p2++ ) return 0;
                    168:                        return 1;
                    169:                case 3:
                    170:                        if ( *p1++ != *p2++ ) return 0;
                    171:                        if ( *p1++ != *p2++ ) return 0;
                    172:                        if ( *p1++ != *p2++ ) return 0;
                    173:                        if ( *p1++ != *p2++ ) return 0;
                    174:                        return 1;
                    175:                case 2:
                    176:                        if ( *p1++ != *p2++ ) return 0;
                    177:                        if ( *p1++ != *p2++ ) return 0;
                    178:                        if ( *p1++ != *p2++ ) return 0;
                    179:                        return 1;
                    180:                case 1:
                    181:                        if ( *p1++ != *p2++ ) return 0;
                    182:                        if ( *p1++ != *p2++ ) return 0;
                    183:                        return 1;
                    184:                case 0:
                    185:                        if ( *p1++ != *p2++ ) return 0;
                    186:                        return 1;
                    187:                default:
                    188:                        return 1;
                    189:        }
1.1       noro      190: }
                    191:
1.5       noro      192: /* b[] should be cleared */
                    193:
1.37      noro      194: void _dpmod_to_vect(DP f,DL *at,int *b)
1.1       noro      195: {
                    196:        int i,nv;
                    197:        MP m;
                    198:
                    199:        nv = f->nv;
                    200:        for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
                    201:                for ( ; !eqdl(nv,m->dl,at[i]); i++ );
                    202:                b[i] = ITOS(m->c);
                    203:        }
                    204: }
                    205:
1.30      noro      206: /* [t,findex] -> tf -> compressed vector */
                    207:
1.37      noro      208: void _tf_to_vect_compress(NODE tf,DL *at,CDP *b)
1.30      noro      209: {
                    210:        int i,j,k,nv,len;
                    211:        DL t,s,d1;
                    212:        DP f;
                    213:        MP m;
                    214:        CDP r;
                    215:
                    216:        t = (DL)BDY(tf);
                    217:        f = ps[(int)BDY(NEXT(tf))];
                    218:
                    219:        nv = f->nv;
                    220:        for ( m = BDY(f), len = 0; m; m = NEXT(m), len++ );
                    221:        r = (CDP)MALLOC(sizeof(struct oCDP));
                    222:        r->len = len;
1.32      noro      223:        r->psindex = (int)BDY(NEXT(tf));
1.33      noro      224:        r->body = (unsigned int *)MALLOC_ATOMIC(sizeof(unsigned int)*len);
1.30      noro      225:
1.34      noro      226:        NEWDL_NOINIT(s,nv);
1.30      noro      227:        for ( m = BDY(f), i = j = 0; m; m = NEXT(m), j++ ) {
                    228:                d1 = m->dl;
                    229:                s->td = t->td+d1->td;
                    230:                for ( k = 0; k < nv; k++ )
                    231:                        s->d[k] = t->d[k]+d1->d[k];
                    232:                for ( ; !eqdl(nv,s,at[i]); i++ );
1.32      noro      233:                r->body[j] = i;
1.30      noro      234:        }
                    235:        *b = r;
                    236: }
                    237:
1.37      noro      238: void dp_to_vect(DP f,DL *at,Q *b)
1.1       noro      239: {
                    240:        int i,nv;
                    241:        MP m;
                    242:
                    243:        nv = f->nv;
                    244:        for ( m = BDY(f), i = 0; m; m = NEXT(m), i++ ) {
                    245:                for ( ; !eqdl(nv,m->dl,at[i]); i++ );
                    246:                b[i] =(Q)m->c;
                    247:        }
                    248: }
                    249:
1.37      noro      250: NODE dp_dllist(DP f)
1.1       noro      251: {
1.2       noro      252:        MP m;
                    253:        NODE mp,mp0;
1.1       noro      254:
1.2       noro      255:        if ( !f )
                    256:                return 0;
1.1       noro      257:        mp0 = 0;
                    258:        for ( m = BDY(f); m; m = NEXT(m) ) {
1.2       noro      259:                NEXTNODE(mp0,mp); BDY(mp) = (pointer)m->dl;
1.1       noro      260:        }
                    261:        NEXT(mp) = 0;
1.2       noro      262:        return mp0;
                    263: }
                    264:
1.37      noro      265: NODE mul_dllist(DL d,DP f)
1.30      noro      266: {
                    267:        MP m;
                    268:        NODE mp,mp0;
                    269:        DL t,d1;
                    270:        int i,nv;
                    271:
                    272:        if ( !f )
                    273:                return 0;
                    274:        nv = NV(f);
                    275:        mp0 = 0;
                    276:        for ( m = BDY(f); m; m = NEXT(m) ) {
                    277:                NEXTNODE(mp0,mp);
1.34      noro      278:                NEWDL_NOINIT(t,nv);
1.30      noro      279:                d1 = m->dl;
                    280:                t->td = d->td+d1->td;
                    281:                for ( i = 0; i < nv; i++ )
                    282:                        t->d[i] = d->d[i]+d1->d[i];
                    283:                BDY(mp) = (pointer)t;
                    284:        }
                    285:        NEXT(mp) = 0;
                    286:        return mp0;
                    287: }
                    288:
1.37      noro      289: void pdl(NODE f)
1.2       noro      290: {
                    291:        while ( f ) {
                    292:                printdl(BDY(f)); f = NEXT(f);
                    293:        }
                    294:        fflush(stdout);
                    295:        printf("\n");
1.1       noro      296: }
                    297:
1.37      noro      298: void dp_gr_main(LIST f,LIST v,Num homo,int modular,int field,struct order_spec *ord,LIST *rp)
1.1       noro      299: {
                    300:        int i,mindex,m,nochk;
                    301:        struct order_spec ord1;
1.23      noro      302:        Q q;
1.1       noro      303:        VL fv,vv,vc;
                    304:        NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
1.23      noro      305:        NODE ind,ind0;
                    306:        LIST trace,gbindex;
1.1       noro      307:
1.21      noro      308:        mindex = 0; nochk = 0; dp_fcoeffs = field;
1.1       noro      309:        get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
                    310:        NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
                    311:        CNVars = homo ? NVars+1 : NVars;
                    312:        if ( ord->id && NVars != ord->nv )
                    313:                error("dp_gr_main : invalid order specification");
                    314:        initd(ord);
                    315:        if ( homo ) {
                    316:                homogenize_order(ord,NVars,&ord1);
                    317:                for ( fd0 = fi0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                    318:                        NEXTNODE(fd0,fd); NEXTNODE(fi0,fi);
                    319:                        ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); dp_homo((DP)BDY(fi),(DP *)&BDY(fd));
                    320:                }
                    321:                if ( fd0 ) NEXT(fd) = 0;
                    322:                if ( fi0 ) NEXT(fi) = 0;
                    323:                initd(&ord1);
                    324:        } else {
                    325:                for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                    326:                        NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
                    327:                }
                    328:                if ( fd0 ) NEXT(fd) = 0;
                    329:                fi0 = fd0;
                    330:        }
                    331:        if ( modular < 0 ) {
                    332:                modular = -modular; nochk = 1;
                    333:        }
                    334:        if ( modular )
1.19      noro      335:                m = modular > 1 ? modular : get_lprime(mindex);
1.1       noro      336:        else
                    337:                m = 0;
                    338:        makesubst(vc,&subst);
                    339:        setup_arrays(fd0,0,&s);
                    340:        init_stat();
                    341:        while ( 1 ) {
                    342:                if ( homo ) {
                    343:                        initd(&ord1); CNVars = NVars+1;
                    344:                }
1.12      noro      345:                if ( DP_Print && modular ) {
1.1       noro      346:                        fprintf(asir_out,"mod= %d, eval = ",m); printsubst(subst);
                    347:                }
                    348:                x = gb(s,m,subst);
                    349:                if ( x ) {
                    350:                        if ( homo ) {
                    351:                                reducebase_dehomo(x,&xx); x = xx;
                    352:                                initd(ord); CNVars = NVars;
                    353:                        }
                    354:                        reduceall(x,&xx); x = xx;
                    355:                        if ( modular ) {
                    356:                                if ( nochk || (membercheck(fi0,x) && gbcheck(x)) )
                    357:                                        break;
                    358:                        } else
                    359:                                break;
                    360:                }
                    361:                if ( modular )
                    362:                        if ( modular > 1 ) {
                    363:                                *rp = 0; return;
                    364:                        } else
1.19      noro      365:                                m = get_lprime(++mindex);
1.1       noro      366:                makesubst(vc,&subst);
                    367:                psn = length(s);
                    368:                for ( i = psn; i < pslen; i++ ) {
                    369:                        pss[i] = 0; psh[i] = 0; psc[i] = 0; ps[i] = 0;
                    370:                }
                    371:        }
1.23      noro      372:        for ( r0 = 0, ind0 = 0; x; x = NEXT(x) ) {
1.1       noro      373:                NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
                    374:                dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
1.23      noro      375:                NEXTNODE(ind0,ind);
                    376:                STOQ((int)BDY(x),q); BDY(ind) = q;
1.1       noro      377:        }
                    378:        if ( r0 ) NEXT(r) = 0;
1.23      noro      379:        if ( ind0 ) NEXT(ind) = 0;
1.1       noro      380:        MKLIST(*rp,r0);
1.23      noro      381:        MKLIST(gbindex,ind0);
                    382:
                    383:        if ( GenTrace && OXCheck < 0 ) {
                    384:
                    385:                x = AllTraceList;
                    386:                for ( r = 0; x; x = NEXT(x) ) {
                    387:                        MKNODE(r0,BDY(x),r); r = r0;
                    388:                }
                    389:                MKLIST(trace,r);
                    390:                r0 = mknode(3,*rp,gbindex,trace);
                    391:                MKLIST(*rp,r0);
                    392:        }
1.1       noro      393:        print_stat();
                    394:        if ( ShowMag )
1.42      noro      395:                fprintf(asir_out,"\nMax_mag=%d, Max_coef=%d\n",Max_mag, Max_coef);
1.1       noro      396: }
                    397:
1.37      noro      398: void dp_gr_mod_main(LIST f,LIST v,Num homo,int m,struct order_spec *ord,LIST *rp)
1.1       noro      399: {
                    400:        struct order_spec ord1;
                    401:        VL fv,vv,vc;
                    402:        NODE fd,fd0,r,r0,t,x,s,xx;
                    403:        DP a,b,c;
1.17      noro      404: extern struct oEGT eg_red_mod;
1.1       noro      405:
                    406:        get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
                    407:        NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
                    408:        CNVars = homo ? NVars+1 : NVars;
                    409:        if ( ord->id && NVars != ord->nv )
                    410:                error("dp_gr_mod_main : invalid order specification");
                    411:        initd(ord);
                    412:        if ( homo ) {
                    413:                for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                    414:                        ptod(CO,vv,(P)BDY(t),&a); dp_homo(a,&b);
                    415:                        if ( PCoeffs )
                    416:                                dp_mod(b,m,0,&c);
                    417:                        else
                    418:                                _dp_mod(b,m,(NODE)0,&c);
                    419:                        if ( c ) {
                    420:                                NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
                    421:                        }
                    422:                }
                    423:                homogenize_order(ord,NVars,&ord1); initd(&ord1);
                    424:        } else {
                    425:                for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
                    426:                        ptod(CO,vv,(P)BDY(t),&b);
                    427:                        if ( PCoeffs )
                    428:                                dp_mod(b,m,0,&c);
                    429:                        else
                    430:                                _dp_mod(b,m,0,&c);
                    431:                        if ( c ) {
                    432:                                NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
                    433:                        }
                    434:                }
                    435:        }
                    436:        if ( fd0 ) NEXT(fd) = 0;
                    437:        setup_arrays(fd0,m,&s);
                    438:        init_stat();
                    439:        if ( homo ) {
                    440:                initd(&ord1); CNVars = NVars+1;
                    441:        }
1.17      noro      442: /* init_eg(&eg_red_mod); */
1.1       noro      443:        x = gb_mod(s,m);
1.17      noro      444: /* print_eg("Red_mod",&eg_red_mod); */
1.1       noro      445:        if ( homo ) {
                    446:                reducebase_dehomo(x,&xx); x = xx;
                    447:                initd(ord); CNVars = NVars;
                    448:        }
                    449:        reduceall_mod(x,m,&xx); x = xx;
                    450:        if ( PCoeffs )
                    451:                for ( r0 = 0; x; x = NEXT(x) ) {
                    452:                        NEXTNODE(r0,r); mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r));
                    453:                }
                    454:        else
                    455:                for ( r0 = 0; x; x = NEXT(x) ) {
                    456:                        NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
                    457:                }
                    458:        print_stat();
                    459:        if ( r0 ) NEXT(r) = 0;
                    460:        MKLIST(*rp,r0);
                    461: }
                    462:
1.37      noro      463: void dp_f4_main(LIST f,LIST v,struct order_spec *ord,LIST *rp)
1.1       noro      464: {
1.37      noro      465:        int homogen;
1.1       noro      466:        VL fv,vv,vc;
1.37      noro      467:        NODE fd,fd0,r,r0,t,x,s,xx;
1.1       noro      468:
1.22      noro      469:        dp_fcoeffs = 0;
1.1       noro      470:        get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
                    471:        NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
                    472:        CNVars = NVars;
                    473:        if ( ord->id && NVars != ord->nv )
                    474:                error("dp_f4_main : invalid order specification");
                    475:        initd(ord);
1.27      noro      476:        for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
1.1       noro      477:                NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
1.27      noro      478:                if ( homogen )
                    479:                        homogen = dp_homogeneous(BDY(fd));
1.1       noro      480:        }
                    481:        if ( fd0 ) NEXT(fd) = 0;
                    482:        setup_arrays(fd0,0,&s);
                    483:        x = gb_f4(s);
1.27      noro      484:        if ( !homogen ) {
                    485:                reduceall(x,&xx); x = xx;
                    486:        }
1.1       noro      487:        for ( r0 = 0; x; x = NEXT(x) ) {
                    488:                NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
                    489:                dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
                    490:        }
                    491:        if ( r0 ) NEXT(r) = 0;
                    492:        MKLIST(*rp,r0);
                    493: }
                    494:
1.37      noro      495: void dp_f4_mod_main(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp)
1.1       noro      496: {
1.37      noro      497:        int homogen;
1.1       noro      498:        VL fv,vv,vc;
1.5       noro      499:        DP b,c,c1;
1.37      noro      500:        NODE fd,fd0,r,r0,t,x,s,xx;
1.1       noro      501:
                    502:        dp_fcoeffs = 0;
                    503:        get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
                    504:        NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
                    505:        CNVars = NVars;
                    506:        if ( ord->id && NVars != ord->nv )
                    507:                error("dp_f4_mod_main : invalid order specification");
                    508:        initd(ord);
1.27      noro      509:        for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
1.1       noro      510:                ptod(CO,vv,(P)BDY(t),&b);
1.27      noro      511:                if ( homogen )
                    512:                        homogen = dp_homogeneous(b);
1.1       noro      513:                _dp_mod(b,m,0,&c);
1.5       noro      514:                _dp_monic(c,m,&c1);
1.1       noro      515:                if ( c ) {
1.5       noro      516:                        NEXTNODE(fd0,fd); BDY(fd) = (pointer)c1;
1.1       noro      517:                }
                    518:        }
                    519:        if ( fd0 ) NEXT(fd) = 0;
                    520:        setup_arrays(fd0,m,&s);
1.35      noro      521:        init_stat();
1.36      noro      522:        if ( do_weyl )
                    523:                x = gb_f4_mod_old(s,m);
                    524:        else
                    525:                x = gb_f4_mod(s,m);
1.27      noro      526:        if ( !homogen ) {
                    527:                reduceall_mod(x,m,&xx); x = xx;
                    528:        }
1.1       noro      529:        for ( r0 = 0; x; x = NEXT(x) ) {
                    530:                NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
                    531:        }
                    532:        if ( r0 ) NEXT(r) = 0;
                    533:        MKLIST(*rp,r0);
1.35      noro      534:        print_stat();
1.1       noro      535: }
                    536:
1.37      noro      537: NODE gb_f4(NODE f)
1.1       noro      538: {
1.37      noro      539:        int i,k,nh,row,col,nv;
1.1       noro      540:        NODE r,g,gall;
1.2       noro      541:        NODE s,s0;
1.1       noro      542:        DP_pairs d,dm,dr,t;
1.37      noro      543:        DP nf,nf1,f2,sp,sd,tdp;
1.1       noro      544:        MP mp,mp0;
1.37      noro      545:        NODE blist,bt;
1.1       noro      546:        DL *ht,*at;
                    547:        MAT mat,nm;
                    548:        int *rind,*cind;
                    549:        int rank,nred;
                    550:        Q dn;
1.37      noro      551:        struct oEGT tmp0,tmp1,eg_split_symb;
1.1       noro      552:        extern struct oEGT eg_mod,eg_elim,eg_chrem,eg_gschk,eg_intrat,eg_symb;
                    553:
                    554:        init_eg(&eg_mod); init_eg(&eg_elim); init_eg(&eg_chrem);
                    555:        init_eg(&eg_gschk); init_eg(&eg_intrat); init_eg(&eg_symb);
                    556:
                    557:        doing_f4 = 1;
                    558:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                    559:                i = (int)BDY(r);
                    560:                d = updpairs(d,g,i);
                    561:                g = updbase(g,i);
                    562:                gall = append_one(gall,i);
                    563:        }
                    564:        if ( gall )
                    565:                nv = ((DP)ps[(int)BDY(gall)])->nv;
                    566:        while ( d ) {
                    567:                get_eg(&tmp0);
                    568:                minsugar(d,&dm,&dr); d = dr;
1.12      noro      569:                if ( DP_Print )
1.1       noro      570:                        fprintf(asir_out,"sugar=%d\n",dm->sugar);
                    571:                blist = 0; s0 = 0;
                    572:                /* asph : sum of all head terms of spoly */
                    573:                for ( t = dm; t; t = NEXT(t) ) {
                    574:                        dp_sp(ps[t->dp1],ps[t->dp2],&sp);
                    575:                        if ( sp ) {
                    576:                                MKNODE(bt,sp,blist); blist = bt;
1.2       noro      577:                                s0 = symb_merge(s0,dp_dllist(sp),nv);
1.1       noro      578:                        }
                    579:                }
                    580:                /* s0 : all the terms appeared in symbolic redunction */
1.2       noro      581:                for ( s = s0, nred = 0; s; s = NEXT(s) ) {
1.1       noro      582:                        for ( r = gall; r; r = NEXT(r) )
1.2       noro      583:                                if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
1.1       noro      584:                                        break;
                    585:                        if ( r ) {
1.2       noro      586:                                dltod(BDY(s),nv,&tdp);
                    587:                                dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.5       noro      588:                                muld(CO,sd,ps[(int)BDY(r)],&f2);
1.1       noro      589:                                MKNODE(bt,f2,blist); blist = bt;
1.2       noro      590:                                s = symb_merge(s,dp_dllist(f2),nv);
1.1       noro      591:                                nred++;
                    592:                        }
                    593:                }
                    594:
                    595:                /* the first nred polys in blist are reducers */
                    596:                /* row = the number of all the polys */
                    597:                for ( r = blist, row = 0; r; r = NEXT(r), row++ );
                    598:                ht = (DL *)MALLOC(nred*sizeof(DL));
                    599:                for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
                    600:                        ht[i] = BDY((DP)BDY(r))->dl;
1.2       noro      601:                for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.1       noro      602:                at = (DL *)MALLOC(col*sizeof(DL));
1.2       noro      603:                for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
                    604:                        at[i] = (DL)BDY(s);
1.1       noro      605:                MKMAT(mat,row,col);
                    606:                for ( i = 0, r = blist; i < row; r = NEXT(r), i++ )
                    607:                        dp_to_vect(BDY(r),at,(Q *)mat->body[i]);
                    608:                get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
                    609:                init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
1.12      noro      610:                if ( DP_Print ) {
1.1       noro      611:                        print_eg("Symb",&eg_split_symb);
                    612:                        fprintf(asir_out,"mat : %d x %d",row,col);
                    613:                        fflush(asir_out);
                    614:                }
1.3       noro      615: #if 0
                    616:                rank = generic_gauss_elim_hensel(mat,&nm,&dn,&rind,&cind);
                    617: #else
1.1       noro      618:                rank = generic_gauss_elim(mat,&nm,&dn,&rind,&cind);
1.3       noro      619: #endif
1.12      noro      620:                if ( DP_Print )
1.1       noro      621:                        fprintf(asir_out,"done rank = %d\n",rank,row,col);
                    622:                for ( i = 0; i < rank; i++ ) {
                    623:                        for ( k = 0; k < nred; k++ )
                    624:                                if ( !cmpdl(nv,at[rind[i]],ht[k]) )
                    625:                                        break;
                    626:                        if ( k == nred ) {
                    627:                                /* this is a new base */
                    628:                                mp0 = 0;
                    629:                                NEXTMP(mp0,mp); mp->dl = at[rind[i]]; mp->c = (P)dn;
                    630:                                for ( k = 0; k < col-rank; k++ )
                    631:                                        if ( nm->body[i][k] ) {
                    632:                                                NEXTMP(mp0,mp); mp->dl = at[cind[k]];
                    633:                                                mp->c = (P)nm->body[i][k];
                    634:                                        }
                    635:                                NEXT(mp) = 0;
                    636:                                MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
                    637:                                dp_ptozp(nf,&nf1);
                    638:                                nh = newps(nf1,0,0);
                    639:                                d = updpairs(d,g,nh);
                    640:                                g = updbase(g,nh);
                    641:                                gall = append_one(gall,nh);
                    642:                        }
                    643:                }
                    644:        }
1.12      noro      645:        if ( DP_Print ) {
1.1       noro      646:                print_eg("Symb",&eg_symb);
                    647:                print_eg("Mod",&eg_mod); print_eg("GaussElim",&eg_elim);
                    648:                print_eg("ChRem",&eg_chrem); print_eg("IntToRat",&eg_intrat);
                    649:                print_eg("Check",&eg_gschk);
                    650:        }
                    651:        return g;
                    652: }
                    653:
1.5       noro      654: /* initial bases are monic */
                    655:
1.32      noro      656: unsigned int **psca;
1.51      noro      657: GeoBucket create_bucket();
                    658: DL remove_head_bucket(GeoBucket,int);
1.32      noro      659:
1.37      noro      660: NODE gb_f4_mod(NODE f,int m)
1.1       noro      661: {
                    662:        int i,j,k,nh,row,col,nv;
                    663:        NODE r,g,gall;
1.29      noro      664:        NODE s,s0;
1.1       noro      665:        DP_pairs d,dm,dr,t;
1.37      noro      666:        DP nf,sp,sd,tdp;
1.1       noro      667:        MP mp,mp0;
1.37      noro      668:        NODE blist,bt,bt1,dt;
                    669:        DL *at,*st;
1.24      noro      670:        int **spmat;
                    671:        CDP *redmat;
1.25      noro      672:        int *colstat,*w,*w1;
1.35      noro      673:        int rank,nred,nsp,nsp0,nonzero,spcol;
1.24      noro      674:        int *indred,*isred;
                    675:        CDP ri;
1.32      noro      676:        int pscalen;
1.51      noro      677:        GeoBucket bucket;
                    678:        DL head;
1.50      noro      679:        struct oEGT tmp0,tmp1,eg_split_symb,eg_split_conv,eg_split_elim1,eg_split_elim2;
                    680:        extern struct oEGT eg_symb,eg_conv,eg_elim1,eg_elim2;
1.1       noro      681:
1.32      noro      682:        /* initialize coeffcient array list of ps[] */
                    683:        pscalen = pslen;
                    684:        psca = (unsigned int **)MALLOC(pscalen*sizeof(unsigned int *));
                    685:
1.50      noro      686:        init_eg(&eg_symb); init_eg(&eg_conv); init_eg(&eg_elim1); init_eg(&eg_elim2);
1.1       noro      687:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                    688:                i = (int)BDY(r);
                    689:                d = updpairs(d,g,i);
                    690:                g = updbase(g,i);
                    691:                gall = append_one(gall,i);
1.32      noro      692:                dptoca(ps[i],&psca[i]);
1.1       noro      693:        }
                    694:        if ( gall )
                    695:                nv = ((DP)ps[(int)BDY(gall)])->nv;
                    696:        while ( d ) {
                    697:                get_eg(&tmp0);
                    698:                minsugar(d,&dm,&dr); d = dr;
1.12      noro      699:                if ( DP_Print )
1.1       noro      700:                        fprintf(asir_out,"sugar=%d\n",dm->sugar);
1.51      noro      701:                blist = 0;
                    702:                bucket = create_bucket();
1.1       noro      703:                /* asph : sum of all head terms of spoly */
                    704:                for ( t = dm; t; t = NEXT(t) ) {
1.13      noro      705:                        _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
1.36      noro      706: /*                     fprintf(stderr,"splen=%d-",dp_nt(sp)); */
1.1       noro      707:                        if ( sp ) {
                    708:                                MKNODE(bt,sp,blist); blist = bt;
1.51      noro      709:                                add_bucket(bucket,dp_dllist(sp),nv);
1.36      noro      710: /*                             fprintf(stderr,"%d-",length(s0)); */
1.1       noro      711:                        }
                    712:                }
1.51      noro      713: #if 0
1.35      noro      714:                if ( DP_Print )
                    715:                        fprintf(asir_out,"initial spmat : %d x %d ",length(blist),length(s0));
1.51      noro      716: #endif
1.24      noro      717:                /* s0 : all the terms appeared in symbolic reduction */
1.51      noro      718:                nred = 0;
                    719:                s0 = 0;
                    720:                while ( 1 ) {
                    721:                        head = remove_head_bucket(bucket,nv);
                    722:                        if ( !head ) break;
                    723:                        else {
                    724:                                NEXTNODE(s0,s);
                    725:                                BDY(s) = (pointer)head;
                    726:                        }
1.29      noro      727:                        for ( r = gall; r; r = NEXT(r) )
1.51      noro      728:                                if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,head,nv) )
1.29      noro      729:                                        break;
                    730:                        if ( r ) {
1.51      noro      731:                                dltod(head,nv,&tdp);
1.29      noro      732:                                dp_subd(tdp,ps[(int)BDY(r)],&sd);
1.30      noro      733:                                dt = mul_dllist(BDY(sd)->dl,ps[(int)BDY(r)]);
1.51      noro      734:                                add_bucket(bucket,NEXT(dt),nv);
1.48      noro      735: /*                             fprintf(stderr,"[%d]",length(dt)); */
1.30      noro      736:                                /* list of [t,f] */
                    737:                                bt1 = mknode(2,BDY(sd)->dl,BDY(r));
                    738:                                MKNODE(bt,bt1,blist); blist = bt;
1.48      noro      739: /*                             fprintf(stderr,"%d-",length(s0));  */
1.29      noro      740:                                nred++;
1.1       noro      741:                        }
                    742:                }
1.51      noro      743:                if ( s0 ) NEXT(s) = 0;
1.36      noro      744: /*             fprintf(stderr,"\n"); */
1.50      noro      745:                get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
                    746:                init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
                    747:
1.35      noro      748:                if ( DP_Print )
                    749:                        fprintf(asir_out,"number of reducers : %d\n",nred);
1.1       noro      750:
1.50      noro      751:                get_eg(&tmp0);
1.1       noro      752:                /* the first nred polys in blist are reducers */
                    753:                /* row = the number of all the polys */
                    754:                for ( r = blist, row = 0; r; r = NEXT(r), row++ );
1.5       noro      755:
                    756:                /* col = number of all terms */
1.29      noro      757:                for ( s = s0, col = 0; s; s = NEXT(s), col++ );
1.5       noro      758:
                    759:                /* head terms of all terms */
1.1       noro      760:                at = (DL *)MALLOC(col*sizeof(DL));
1.29      noro      761:                for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
                    762:                        at[i] = (DL)BDY(s);
1.5       noro      763:
                    764:                /* store coefficients separately in spmat and redmat */
                    765:                nsp = row-nred;
                    766:
                    767:                /* reducer matrix */
1.25      noro      768:                /* indred : register the position of the head term */
1.24      noro      769:                redmat = (CDP *)MALLOC(nred*sizeof(CDP));
1.5       noro      770:                for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
1.30      noro      771:                        _tf_to_vect_compress(BDY(r),at,&redmat[i]);
1.32      noro      772:
1.5       noro      773:                /* register the position of the head term */
1.31      noro      774:                indred = (int *)MALLOC_ATOMIC(nred*sizeof(int));
1.5       noro      775:                bzero(indred,nred*sizeof(int));
1.31      noro      776:                isred = (int *)MALLOC_ATOMIC(col*sizeof(int));
1.5       noro      777:                bzero(isred,col*sizeof(int));
                    778:                for ( i = 0; i < nred; i++ ) {
                    779:                        ri = redmat[i];
1.32      noro      780:                        indred[i] = ri->body[0];
1.24      noro      781:                        isred[indred[i]] = 1;
1.5       noro      782:                }
                    783:
                    784:                spcol = col-nred;
                    785:                /* head terms not in ht */
                    786:                st = (DL *)MALLOC(spcol*sizeof(DL));
                    787:                for ( j = 0, k = 0; j < col; j++ )
                    788:                        if ( !isred[j] )
                    789:                                st[k++] = at[j];
1.50      noro      790:                get_eg(&tmp1); add_eg(&eg_conv,&tmp0,&tmp1);
                    791:                init_eg(&eg_split_conv); add_eg(&eg_split_conv,&tmp0,&tmp1);
1.5       noro      792:
1.24      noro      793:                get_eg(&tmp1);
1.5       noro      794:                /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
1.25      noro      795:                spmat = (int **)MALLOC(nsp*sizeof(int *));
1.31      noro      796:                w = (int *)MALLOC_ATOMIC(col*sizeof(int));
1.25      noro      797:
                    798:                /* skip reducers in blist */
                    799:                for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ );
                    800:                for ( i = 0; r; r = NEXT(r) ) {
1.5       noro      801:                        bzero(w,col*sizeof(int));
                    802:                        _dpmod_to_vect(BDY(r),at,w);
1.24      noro      803:                        reduce_sp_by_red_mod_compress(w,redmat,indred,nred,col,m);
1.25      noro      804:                        for ( j = 0; j < col; j++ )
                    805:                                if ( w[j] )
                    806:                                        break;
                    807:                        if ( j < col ) {
                    808:                                w1 = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
                    809:                                for ( j = 0, k = 0; j < col; j++ )
                    810:                                        if ( !isred[j] )
                    811:                                                w1[k++] = w[j];
                    812:                                spmat[i] = w1;
                    813:                                i++;
                    814:                        }
1.5       noro      815:                }
1.25      noro      816:                /* update nsp */
1.35      noro      817:                nsp0 = nsp;
1.25      noro      818:                nsp = i;
1.5       noro      819:
1.30      noro      820:                /* XXX free redmat explicitly */
                    821:                for ( k = 0; k < nred; k++ ) {
                    822:                        GC_free(BDY(redmat[k]));
                    823:                        GC_free(redmat[k]);
                    824:                }
                    825:
1.5       noro      826:                get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
                    827:                init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
                    828:
                    829:                colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
1.25      noro      830:                bzero(colstat,spcol*sizeof(int));
1.5       noro      831:                for ( i = 0, nonzero=0; i < nsp; i++ )
                    832:                        for ( j = 0; j < spcol; j++ )
                    833:                                if ( spmat[i][j] )
1.1       noro      834:                                        nonzero++;
1.17      noro      835:                if ( DP_Print && nsp )
1.5       noro      836:                        fprintf(asir_out,"spmat : %d x %d (nonzero=%f%%)...",
                    837:                                nsp,spcol,((double)nonzero*100)/(nsp*spcol));
1.17      noro      838:                if ( nsp )
                    839:                        rank = generic_gauss_elim_mod(spmat,nsp,spcol,m,colstat);
                    840:                else
                    841:                        rank = 0;
1.5       noro      842:                get_eg(&tmp1); add_eg(&eg_elim2,&tmp0,&tmp1);
                    843:                init_eg(&eg_split_elim2); add_eg(&eg_split_elim2,&tmp0,&tmp1);
                    844:
1.12      noro      845:                if ( DP_Print ) {
1.1       noro      846:                        fprintf(asir_out,"done rank = %d\n",rank,row,col);
                    847:                        print_eg("Symb",&eg_split_symb);
1.50      noro      848:                        print_eg("Conv",&eg_split_conv);
1.5       noro      849:                        print_eg("Elim1",&eg_split_elim1);
                    850:                        print_eg("Elim2",&eg_split_elim2);
1.1       noro      851:                        fprintf(asir_out,"\n");
                    852:                }
1.25      noro      853:
1.35      noro      854:                NZR += rank;
                    855:                ZR += nsp0-rank;
                    856:
1.25      noro      857:                if ( !rank )
                    858:                        continue;
                    859:
1.5       noro      860:                for ( j = 0, i = 0; j < spcol; j++ )
1.1       noro      861:                        if ( colstat[j] ) {
1.5       noro      862:                                mp0 = 0;
1.13      noro      863:                                NEXTMP(mp0,mp); mp->dl = st[j]; mp->c = STOI(1);
1.5       noro      864:                                for ( k = j+1; k < spcol; k++ )
                    865:                                        if ( !colstat[k] && spmat[i][k] ) {
1.13      noro      866:                                                NEXTMP(mp0,mp); mp->dl = st[k];
1.5       noro      867:                                                mp->c = STOI(spmat[i][k]);
1.1       noro      868:                                }
1.5       noro      869:                                NEXT(mp) = 0;
                    870:                                MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
                    871:                                nh = newps_mod(nf,m);
1.32      noro      872:                                if ( nh == pscalen ) {
                    873:                                        psca = (unsigned int **)
                    874:                                                REALLOC(psca,2*pscalen*sizeof(unsigned int *));
                    875:                                        pscalen *= 2;
                    876:                                }
                    877:                                dptoca(ps[nh],&psca[nh]);
1.5       noro      878:                                d = updpairs(d,g,nh);
                    879:                                g = updbase(g,nh);
                    880:                                gall = append_one(gall,nh);
1.1       noro      881:                                i++;
                    882:                        }
1.30      noro      883:
                    884:                /* XXX free spmat[] explicitly */
                    885:                for ( j = 0; j < nsp; j++ ) {
                    886:                        GC_free(spmat[j]);
                    887:                }
1.36      noro      888:        }
                    889:        if ( DP_Print ) {
                    890:                print_eg("Symb",&eg_symb);
1.50      noro      891:                print_eg("Conv",&eg_conv);
1.36      noro      892:                print_eg("Elim1",&eg_elim1);
                    893:                print_eg("Elim2",&eg_elim2);
                    894:                fflush(asir_out);
                    895:        }
                    896:        return g;
                    897: }
                    898:
1.37      noro      899: NODE gb_f4_mod_old(NODE f,int m)
1.36      noro      900: {
                    901:        int i,j,k,nh,row,col,nv;
                    902:        NODE r,g,gall;
                    903:        NODE s,s0;
                    904:        DP_pairs d,dm,dr,t;
1.37      noro      905:        DP nf,f2,sp,sd,sdm,tdp;
1.36      noro      906:        MP mp,mp0;
1.37      noro      907:        NODE blist,bt;
1.36      noro      908:        DL *ht,*at,*st;
                    909:        int **spmat,**redmat;
                    910:        int *colstat,*w;
                    911:        int rank,nred,nsp,nonzero,spcol;
                    912:        int *indred,*isred,*ri;
1.37      noro      913:        struct oEGT tmp0,tmp1,eg_split_symb,eg_split_elim1,eg_split_elim2;
1.36      noro      914:        extern struct oEGT eg_symb,eg_elim1,eg_elim2;
                    915:
                    916:        init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);
                    917:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                    918:                i = (int)BDY(r);
                    919:                d = updpairs(d,g,i);
                    920:                g = updbase(g,i);
                    921:                gall = append_one(gall,i);
                    922:        }
                    923:        if ( gall )
                    924:                nv = ((DP)ps[(int)BDY(gall)])->nv;
                    925:        while ( d ) {
                    926:                get_eg(&tmp0);
                    927:                minsugar(d,&dm,&dr); d = dr;
                    928:                if ( DP_Print )
                    929:                        fprintf(asir_out,"sugar=%d\n",dm->sugar);
                    930:                blist = 0; s0 = 0;
                    931:                /* asph : sum of all head terms of spoly */
                    932:                for ( t = dm; t; t = NEXT(t) ) {
                    933:                        _dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
                    934:                        if ( sp ) {
                    935:                                MKNODE(bt,sp,blist); blist = bt;
                    936:                                s0 = symb_merge(s0,dp_dllist(sp),nv);
                    937:                        }
                    938:                }
                    939:                /* s0 : all the terms appeared in symbolic redunction */
                    940:                for ( s = s0, nred = 0; s; s = NEXT(s) ) {
                    941:                        for ( r = gall; r; r = NEXT(r) )
                    942:                                if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
                    943:                                        break;
                    944:                        if ( r ) {
                    945:                                dltod(BDY(s),nv,&tdp);
                    946:                                dp_subd(tdp,ps[(int)BDY(r)],&sd);
                    947:                                _dp_mod(sd,m,0,&sdm);
                    948:                                mulmd_dup(m,sdm,ps[(int)BDY(r)],&f2);
                    949:                                MKNODE(bt,f2,blist); blist = bt;
                    950:                                s = symb_merge(s,dp_dllist(f2),nv);
                    951:                                nred++;
                    952:                        }
                    953:                }
                    954:
                    955:                get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
                    956:                init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
                    957:
                    958:                /* the first nred polys in blist are reducers */
                    959:                /* row = the number of all the polys */
                    960:                for ( r = blist, row = 0; r; r = NEXT(r), row++ );
                    961:
                    962:                /* head terms of reducers */
                    963:                ht = (DL *)MALLOC(nred*sizeof(DL));
                    964:                for ( r = blist, i = 0; i < nred; r = NEXT(r), i++ )
                    965:                        ht[i] = BDY((DP)BDY(r))->dl;
                    966:
                    967:                /* col = number of all terms */
                    968:                for ( s = s0, col = 0; s; s = NEXT(s), col++ );
                    969:
                    970:                /* head terms of all terms */
                    971:                at = (DL *)MALLOC(col*sizeof(DL));
                    972:                for ( s = s0, i = 0; i < col; s = NEXT(s), i++ )
                    973:                        at[i] = (DL)BDY(s);
                    974:
                    975:                /* store coefficients separately in spmat and redmat */
                    976:                nsp = row-nred;
                    977:
                    978:                /* reducer matrix */
                    979:                redmat = (int **)almat(nred,col);
                    980:                for ( i = 0, r = blist; i < nred; r = NEXT(r), i++ )
                    981:                        _dpmod_to_vect(BDY(r),at,redmat[i]);
                    982:                /* XXX */
                    983: /*             reduce_reducers_mod(redmat,nred,col,m); */
                    984:                /* register the position of the head term */
                    985:                indred = (int *)MALLOC(nred*sizeof(int));
                    986:                bzero(indred,nred*sizeof(int));
                    987:                isred = (int *)MALLOC(col*sizeof(int));
                    988:                bzero(isred,col*sizeof(int));
                    989:                for ( i = 0; i < nred; i++ ) {
                    990:                        ri = redmat[i];
                    991:                        for ( j = 0; j < col && !ri[j]; j++ );
                    992:                        indred[i] = j;
                    993:                        isred[j] = 1;
                    994:                }
                    995:
                    996:                spcol = col-nred;
                    997:                /* head terms not in ht */
                    998:                st = (DL *)MALLOC(spcol*sizeof(DL));
                    999:                for ( j = 0, k = 0; j < col; j++ )
                   1000:                        if ( !isred[j] )
                   1001:                                st[k++] = at[j];
                   1002:
                   1003:                /* spoly matrix; stored in reduced form; terms in ht[] are omitted */
                   1004:                spmat = almat(nsp,spcol);
                   1005:                w = (int *)MALLOC(col*sizeof(int));
                   1006:                for ( ; i < row; r = NEXT(r), i++ ) {
                   1007:                        bzero(w,col*sizeof(int));
                   1008:                        _dpmod_to_vect(BDY(r),at,w);
                   1009:                        reduce_sp_by_red_mod(w,redmat,indred,nred,col,m);
                   1010:                        for ( j = 0, k = 0; j < col; j++ )
                   1011:                                if ( !isred[j] )
                   1012:                                        spmat[i-nred][k++] = w[j];
                   1013:                }
                   1014:
                   1015:                get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
                   1016:                init_eg(&eg_split_elim1); add_eg(&eg_split_elim1,&tmp1,&tmp0);
                   1017:
                   1018:                colstat = (int *)MALLOC_ATOMIC(spcol*sizeof(int));
                   1019:                for ( i = 0, nonzero=0; i < nsp; i++ )
                   1020:                        for ( j = 0; j < spcol; j++ )
                   1021:                                if ( spmat[i][j] )
                   1022:                                        nonzero++;
                   1023:                if ( DP_Print && nsp )
                   1024:                        fprintf(asir_out,"spmat : %d x %d (nonzero=%f%%)...",
                   1025:                                nsp,spcol,((double)nonzero*100)/(nsp*spcol));
                   1026:                if ( nsp )
                   1027:                        rank = generic_gauss_elim_mod(spmat,nsp,spcol,m,colstat);
                   1028:                else
                   1029:                        rank = 0;
                   1030:                get_eg(&tmp1); add_eg(&eg_elim2,&tmp0,&tmp1);
                   1031:                init_eg(&eg_split_elim2); add_eg(&eg_split_elim2,&tmp0,&tmp1);
                   1032:
                   1033:                if ( DP_Print ) {
                   1034:                        fprintf(asir_out,"done rank = %d\n",rank,row,col);
                   1035:                        print_eg("Symb",&eg_split_symb);
                   1036:                        print_eg("Elim1",&eg_split_elim1);
                   1037:                        print_eg("Elim2",&eg_split_elim2);
                   1038:                        fprintf(asir_out,"\n");
                   1039:                }
                   1040:                for ( j = 0, i = 0; j < spcol; j++ )
                   1041:                        if ( colstat[j] ) {
                   1042:                                mp0 = 0;
                   1043:                                NEXTMP(mp0,mp); mp->dl = st[j]; mp->c = STOI(1);
                   1044:                                for ( k = j+1; k < spcol; k++ )
                   1045:                                        if ( !colstat[k] && spmat[i][k] ) {
                   1046:                                                NEXTMP(mp0,mp); mp->dl = st[k];
                   1047:                                                mp->c = STOI(spmat[i][k]);
                   1048:                                }
                   1049:                                NEXT(mp) = 0;
                   1050:                                MKDP(nv,mp0,nf); nf->sugar = dm->sugar;
                   1051:                                nh = newps_mod(nf,m);
                   1052:                                d = updpairs(d,g,nh);
                   1053:                                g = updbase(g,nh);
                   1054:                                gall = append_one(gall,nh);
                   1055:                                i++;
                   1056:                        }
1.1       noro     1057:        }
1.12      noro     1058:        if ( DP_Print ) {
1.1       noro     1059:                print_eg("Symb",&eg_symb);
1.5       noro     1060:                print_eg("Elim1",&eg_elim1);
                   1061:                print_eg("Elim2",&eg_elim2);
1.1       noro     1062:                fflush(asir_out);
                   1063:        }
                   1064:        return g;
                   1065: }
                   1066:
1.37      noro     1067: int DPPlength(DP_pairs n)
1.1       noro     1068: {
                   1069:        int i;
                   1070:
                   1071:        for ( i = 0; n; n = NEXT(n), i++ );
                   1072:        return i;
                   1073: }
                   1074:
1.37      noro     1075: void printdl(DL dl)
1.1       noro     1076: {
                   1077:        int i;
                   1078:
                   1079:        fprintf(asir_out,"<<");
                   1080:        for ( i = 0; i < CNVars-1; i++ )
                   1081:                fprintf(asir_out,"%d,",dl->d[i]);
                   1082:        fprintf(asir_out,"%d>>",dl->d[i]);
                   1083: }
                   1084:
1.37      noro     1085: void pltovl(LIST l,VL *vl)
1.1       noro     1086: {
                   1087:        NODE n;
                   1088:        VL r,r0;
                   1089:
                   1090:        n = BDY(l);
                   1091:        for ( r0 = 0; n; n = NEXT(n) ) {
                   1092:                NEXTVL(r0,r); r->v = VR((P)BDY(n));
                   1093:        }
                   1094:        if ( r0 ) NEXT(r) = 0;
                   1095:        *vl = r0;
                   1096: }
                   1097:
1.37      noro     1098: void makesubst(VL v,NODE *s)
1.1       noro     1099: {
                   1100:        NODE r,r0;
                   1101:        Q q;
                   1102:        unsigned int n;
                   1103:
                   1104:        for ( r0 = 0; v; v = NEXT(v) ) {
                   1105:                NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
                   1106: #if defined(_PA_RISC1_1)
                   1107:                n = mrand48()&BMASK; UTOQ(n,q);
                   1108: #else
                   1109:                n = random(); UTOQ(n,q);
                   1110: #endif
                   1111:                NEXTNODE(r0,r); BDY(r) = (pointer)q;
                   1112:        }
                   1113:        if ( r0 ) NEXT(r) = 0;
                   1114:        *s = r0;
                   1115: }
                   1116:
1.37      noro     1117: void printsubst(NODE s)
1.1       noro     1118: {
                   1119:        fputc('[',asir_out);
                   1120:        while ( s ) {
                   1121:                printv(CO,(V)BDY(s)); s = NEXT(s);
                   1122:                fprintf(asir_out,"->%d",QTOS((Q)BDY(s)));
                   1123:                if ( NEXT(s) ) {
                   1124:                        fputc(',',asir_out); s = NEXT(s);
                   1125:                } else
                   1126:                        break;
                   1127:        }
                   1128:        fprintf(asir_out,"]\n"); return;
                   1129: }
                   1130:
1.37      noro     1131: void vlminus(VL v,VL w,VL *d)
1.1       noro     1132: {
                   1133:        int i,j,n,m;
                   1134:        V *va,*wa;
                   1135:        V a;
                   1136:        VL r,r0;
                   1137:        VL t;
                   1138:
                   1139:        for ( n = 0, t = v; t; t = NEXT(t), n++ );
                   1140:        va = (V *)ALLOCA(n*sizeof(V));
                   1141:        for ( i = 0, t = v; t; t = NEXT(t), i++ )
                   1142:                va[i] = t->v;
                   1143:        for ( m = 0, t = w; t; t = NEXT(t), m++ );
                   1144:        wa = (V *)ALLOCA(m*sizeof(V));
                   1145:        for ( i = 0, t = w; t; t = NEXT(t), i++ )
                   1146:                wa[i] = t->v;
                   1147:        for ( i = 0; i < n; i++ ) {
                   1148:                a = va[i];
                   1149:                for ( j = 0; j < m; j++ )
                   1150:                        if ( a == wa[j] )
                   1151:                                break;
                   1152:                if ( j < m )
                   1153:                        va[i] = 0;
                   1154:        }
                   1155:        for ( r0 = 0, i = 0; i < n; i++ )
                   1156:                if ( va[i] ) { NEXTVL(r0,r); r->v = va[i];      }
                   1157:        if ( r0 ) NEXT(r) = 0;
                   1158:        *d = r0;
                   1159: }
                   1160:
1.37      noro     1161: int validhc(P a,int m,NODE s)
1.1       noro     1162: {
                   1163:        P c,c1;
                   1164:        V v;
                   1165:
                   1166:        if ( !a )
                   1167:                return 0;
                   1168:        for ( c = a; s; s = NEXT(s) ) {
                   1169:                v = (V)BDY(s); s = NEXT(s);
                   1170:                substp(CO,c,v,(P)BDY(s),&c1); c = c1;
                   1171:        }
                   1172:        ptomp(m,c,&c1);
                   1173:        return c1 ? 1 : 0;
                   1174: }
                   1175:
1.37      noro     1176: void setup_arrays(NODE f,int m,NODE *r)
1.1       noro     1177: {
                   1178:        int i;
1.7       noro     1179:        NODE s,s0,f0;
1.1       noro     1180:
1.17      noro     1181: #if 1
1.7       noro     1182:        f0 = f = NODE_sortb(f,1);
1.17      noro     1183: #else
                   1184:        f0 = f;
                   1185: #endif
1.1       noro     1186:        psn = length(f); pslen = 2*psn;
                   1187:        ps = (DP *)MALLOC(pslen*sizeof(DP));
                   1188:        psh = (DL *)MALLOC(pslen*sizeof(DL));
                   1189:        pss = (int *)MALLOC(pslen*sizeof(int));
                   1190:        psc = (P *)MALLOC(pslen*sizeof(P));
                   1191:        for ( i = 0; i < psn; i++, f = NEXT(f) ) {
                   1192:                prim_part((DP)BDY(f),m,&ps[i]);
                   1193:                if ( Demand )
1.7       noro     1194:                        dp_save(i,(Obj)ps[i],0);
1.1       noro     1195:                psh[i] = BDY(ps[i])->dl;
                   1196:                pss[i] = ps[i]->sugar;
                   1197:                psc[i] = BDY(ps[i])->c;
                   1198:        }
1.23      noro     1199:        if ( GenTrace ) {
1.7       noro     1200:                Q q;
                   1201:                STRING fname;
                   1202:                LIST input;
1.23      noro     1203:                NODE arg,t,t1;
1.37      noro     1204:                Obj obj;
1.23      noro     1205:
                   1206:                t = 0;
                   1207:                for ( i = psn-1; i >= 0; i-- ) {
                   1208:                        MKNODE(t1,ps[i],t);
                   1209:                        t = t1;
                   1210:                }
                   1211:                MKLIST(input,t);
1.7       noro     1212:
1.23      noro     1213:                if ( OXCheck >= 0 ) {
                   1214:                        STOQ(OXCheck,q);
                   1215:                        MKSTR(fname,"register_input");
                   1216:                        arg = mknode(3,q,fname,input);
1.37      noro     1217:                        Pox_cmo_rpc(arg,&obj);
1.23      noro     1218:                } else if ( OXCheck < 0 ) {
                   1219:                        MKNODE(AllTraceList,input,0);
                   1220:                }
1.7       noro     1221:        }
1.1       noro     1222:        for ( s0 = 0, i = 0; i < psn; i++ ) {
                   1223:                NEXTNODE(s0,s); BDY(s) = (pointer)i;
                   1224:        }
                   1225:        if ( s0 ) NEXT(s) = 0;
                   1226:        *r = s0;
                   1227: }
                   1228:
1.37      noro     1229: void prim_part(DP f,int m,DP *r)
1.1       noro     1230: {
1.7       noro     1231:        P d,t;
                   1232:
1.1       noro     1233:        if ( m > 0 ) {
                   1234:                if ( PCoeffs )
                   1235:                        dp_prim_mod(f,m,r);
                   1236:                else
1.49      noro     1237:                        _dp_monic(f,m,r);
1.1       noro     1238:        } else {
1.44      noro     1239:                if ( dp_fcoeffs || PCoeffs )
1.1       noro     1240:                        dp_prim(f,r);
                   1241:                else
                   1242:                        dp_ptozp(f,r);
1.7       noro     1243:                if ( GenTrace && TraceList ) {
1.23      noro     1244:                        /* adust the denominator according to the final
                   1245:                           content reduction */
1.7       noro     1246:                        divsp(CO,BDY(f)->c,BDY(*r)->c,&d);
                   1247:                        mulp(CO,(P)ARG3(BDY((LIST)BDY(TraceList))),d,&t);
                   1248:                        ARG3(BDY((LIST)BDY(TraceList))) = t;
                   1249:                }
1.1       noro     1250:        }
                   1251: }
                   1252:
1.37      noro     1253: NODE /* of DP */ NODE_sortb_insert( DP newdp, NODE /* of DP */ nd, int dec )
1.1       noro     1254: {
                   1255:        register NODE last, p;
                   1256:        register DL newdl = BDY(newdp)->dl;
                   1257:        register int (*cmpfun)() = cmpdl, nv = CNVars;
                   1258:        NODE newnd;
                   1259:        int sgn = dec ? 1 : -1;
                   1260:        MKNODE( newnd, newdp, 0 );
                   1261:        if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(last))->dl ) > 0 ) {
                   1262:                NEXT(newnd) = last;
                   1263:                return newnd;
                   1264:        }
                   1265:        for ( ; p = NEXT(last); last = p )
                   1266:                if ( sgn*(*cmpfun)( nv, newdl, BDY((DP) BDY(p))->dl ) > 0 ) break;
                   1267:        if ( p ) NEXT(NEXT(last) = newnd) = p;
                   1268:        else NEXT(last) = newnd;
                   1269:        return nd;
                   1270: }
                   1271:
1.37      noro     1272: NODE NODE_sortb( NODE node, int dec )
1.1       noro     1273: {
                   1274:        register NODE nd, ans;
                   1275:
                   1276:        for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
                   1277:                ans = NODE_sortb_insert( (DP) BDY(nd), ans, dec );
                   1278:        return ans;
                   1279: }
                   1280:
1.37      noro     1281: NODE /* of index */ NODE_sortbi_insert( int newdpi, NODE /* of index */ nd, int dec )
1.1       noro     1282: {
                   1283:        register NODE last, p;
                   1284:        register DL newdl = psh[newdpi];
                   1285:        register int (*cmpfun)() = cmpdl, nv = CNVars;
                   1286:        NODE newnd;
                   1287:        int sgn = dec ? 1 : -1;
                   1288:        MKNODE( newnd, newdpi, 0 );
                   1289:        if ( !(last = nd) || sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(last)] ) > 0 ) {
                   1290:                NEXT(newnd) = last;
                   1291:                return newnd;
                   1292:        }
                   1293:        for ( ; p = NEXT(last); last = p )
                   1294:                if ( sgn*(*cmpfun)( nv, newdl, psh[(int)BDY(p)] ) > 0 ) break;
                   1295:        if ( p ) NEXT(NEXT(last) = newnd) = p;
                   1296:        else NEXT(last) = newnd;
                   1297:        return nd;
                   1298: }
                   1299:
1.37      noro     1300: NODE NODE_sortbi( NODE node, int dec )
1.1       noro     1301: {
                   1302:        register NODE nd, ans;
                   1303:
                   1304:        for ( ans = 0, nd = node; nd; nd = NEXT(nd) )
                   1305:                ans = NODE_sortbi_insert( (int) BDY(nd), ans, dec );
                   1306:        return ans;
                   1307: }
                   1308:
1.37      noro     1309: void reduceall(NODE in,NODE *h)
1.1       noro     1310: {
                   1311:        NODE r,t,top;
                   1312:        int n,i,j;
                   1313:        int *w;
                   1314:        DP g,g1;
                   1315:        struct oEGT tmp0,tmp1;
                   1316:
                   1317:        if ( NoRA ) {
                   1318:                *h = in; return;
                   1319:        }
1.12      noro     1320:        if ( DP_Print || DP_PrintShort ) {
1.1       noro     1321:                fprintf(asir_out,"reduceall\n"); fflush(asir_out);
                   1322:        }
                   1323:        r = NODE_sortbi(in,0);
                   1324:        n = length(r);
                   1325:        w = (int *)ALLOCA(n*sizeof(int));
                   1326:        for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
                   1327:                w[i] = (int)BDY(t);
                   1328:        for ( i = 0; i < n; i++ ) {
                   1329:                for ( top = 0, j = n-1; j >= 0; j-- )
                   1330:                        if ( j != i ) {
                   1331:                                MKNODE(t,(pointer)w[j],top); top = t;
                   1332:                        }
                   1333:                get_eg(&tmp0);
                   1334:                dp_load(w[i],&ps[w[i]]);
                   1335:
1.7       noro     1336:                if ( GenTrace ) {
                   1337:                        Q q;
                   1338:                        NODE node;
                   1339:                        LIST hist;
                   1340:
                   1341:                        STOQ(w[i],q);
                   1342:                        node = mknode(4,ONE,q,ONE,ONE);
                   1343:                        MKLIST(hist,node);
                   1344:                        MKNODE(TraceList,hist,0);
                   1345:                }
1.16      noro     1346:                _dp_nf(top,ps[w[i]],ps,1,&g);
1.1       noro     1347:                prim_part(g,0,&g1);
                   1348:                get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12      noro     1349:                if ( DP_Print || DP_PrintShort ) {
1.1       noro     1350:                        fprintf(asir_out,"."); fflush(asir_out);
                   1351:                }
                   1352:                w[i] = newps(g1,0,(NODE)0);
                   1353:        }
                   1354:        for ( top = 0, j = n-1; j >= 0; j-- ) {
                   1355:                MKNODE(t,(pointer)w[j],top); top = t;
                   1356:        }
                   1357:        *h = top;
1.12      noro     1358:        if ( DP_Print || DP_PrintShort )
1.1       noro     1359:                fprintf(asir_out,"\n");
                   1360: }
                   1361:
1.37      noro     1362: void reduceall_mod(NODE in,int m,NODE *h)
1.1       noro     1363: {
                   1364:        NODE r,t,top;
                   1365:        int n,i,j;
                   1366:        int *w;
1.12      noro     1367:        DP g,p;
1.1       noro     1368:        struct oEGT tmp0,tmp1;
                   1369:
                   1370:        if ( NoRA ) {
                   1371:                *h = in; return;
                   1372:        }
1.12      noro     1373:        if ( DP_Print || DP_PrintShort ) {
1.1       noro     1374:                fprintf(asir_out,"reduceall\n"); fflush(asir_out);
                   1375:        }
                   1376:        r = NODE_sortbi(in,0);
                   1377:        n = length(r);
                   1378:        w = (int *)ALLOCA(n*sizeof(int));
                   1379:        for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
                   1380:                w[i] = (int)BDY(t);
                   1381:        for ( i = 0; i < n; i++ ) {
                   1382:                for ( top = 0, j = n-1; j >= 0; j-- )
                   1383:                        if ( j != i ) {
                   1384:                                MKNODE(t,(pointer)w[j],top); top = t;
                   1385:                        }
                   1386:                get_eg(&tmp0);
                   1387:                if ( PCoeffs )
                   1388:                        dp_nf_mod(top,ps[w[i]],ps,m,1,&g);
1.12      noro     1389:                else {
                   1390:                        dpto_dp(ps[w[i]],&p);
                   1391:                        _dp_nf_mod_destructive(top,p,ps,m,1,&g);
                   1392:                }
1.1       noro     1393:                get_eg(&tmp1); add_eg(&eg_ra,&tmp0,&tmp1);
1.12      noro     1394:                if ( DP_Print || DP_PrintShort ) {
1.1       noro     1395:                        fprintf(asir_out,"."); fflush(asir_out);
                   1396:                }
                   1397:                w[i] = newps_mod(g,m);
                   1398:        }
                   1399:        for ( top = 0, j = n-1; j >= 0; j-- ) {
                   1400:                MKNODE(t,(pointer)w[j],top); top = t;
                   1401:        }
                   1402:        *h = top;
1.12      noro     1403:        if ( DP_Print || DP_PrintShort )
1.1       noro     1404:                fprintf(asir_out,"\n");
                   1405: }
                   1406:
1.37      noro     1407: int newps(DP a,int m,NODE subst)
1.1       noro     1408: {
                   1409:        if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
                   1410:                return -1;
                   1411:        if ( psn == pslen ) {
                   1412:                pslen *= 2;
                   1413:                ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
                   1414:                psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
                   1415:                pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
                   1416:                psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
                   1417:                if ( m )
                   1418:                        psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
                   1419:        }
                   1420:        if ( Demand ) {
                   1421:                if ( doing_f4 )
                   1422:                        ps[psn] = a;
                   1423:                else
                   1424:                        ps[psn] = 0;
1.7       noro     1425:                dp_save(psn,(Obj)a,0);
1.1       noro     1426:        } else
                   1427:                ps[psn] = a;
                   1428:        psh[psn] = BDY(a)->dl;
                   1429:        pss[psn] = a->sugar;
                   1430:        psc[psn] = BDY(a)->c;
                   1431:        if ( m )
                   1432:                _dp_mod(a,m,subst,&psm[psn]);
1.7       noro     1433:        if ( GenTrace ) {
                   1434:                NODE tn,tr,tr1;
1.23      noro     1435:                LIST trace,trace1;
                   1436:                NODE arg;
                   1437:                Q q1,q2;
                   1438:                STRING fname;
1.37      noro     1439:                Obj obj;
1.7       noro     1440:
                   1441:                /* reverse the TraceList */
                   1442:                tn = TraceList;
                   1443:                for ( tr = 0; tn; tn = NEXT(tn) ) {
                   1444:                        MKNODE(tr1,BDY(tn),tr); tr = tr1;
                   1445:                }
                   1446:                MKLIST(trace,tr);
                   1447:                if ( OXCheck >= 0 ) {
                   1448:                        STOQ(OXCheck,q1);
                   1449:                        MKSTR(fname,"check_trace");
                   1450:                        STOQ(psn,q2);
                   1451:                        arg = mknode(5,q1,fname,a,q2,trace);
1.37      noro     1452:                        Pox_cmo_rpc(arg,&obj);
1.23      noro     1453:                } else if ( OXCheck < 0 ) {
                   1454:                        STOQ(psn,q1);
                   1455:                        tn = mknode(2,q1,trace);
                   1456:                        MKLIST(trace1,tn);
                   1457:                        MKNODE(tr,trace1,AllTraceList);
                   1458:                        AllTraceList = tr;
1.7       noro     1459:                } else
                   1460:                        dp_save(psn,(Obj)trace,"t");
                   1461:                TraceList = 0;
                   1462:        }
1.1       noro     1463:        return psn++;
                   1464: }
                   1465:
1.37      noro     1466: int newps_nosave(DP a,int m,NODE subst)
1.1       noro     1467: {
                   1468:        if ( m && !validhc(!a?0:BDY(a)->c,m,subst) )
                   1469:                return -1;
                   1470:        if ( psn == pslen ) {
                   1471:                pslen *= 2;
                   1472:                ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
                   1473:                psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
                   1474:                pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
                   1475:                psc = (P *)REALLOC((char *)psc,pslen*sizeof(P));
                   1476:                if ( m )
                   1477:                        psm = (DP *)REALLOC((char *)psm,pslen*sizeof(DP));
                   1478:        }
                   1479:        ps[psn] = 0;
                   1480:        psh[psn] = BDY(a)->dl;
                   1481:        pss[psn] = a->sugar;
                   1482:        psc[psn] = BDY(a)->c;
                   1483:        if ( m )
                   1484:                _dp_mod(a,m,subst,&psm[psn]);
                   1485:        return psn++;
                   1486: }
                   1487:
1.37      noro     1488: int newps_mod(DP a,int m)
1.1       noro     1489: {
                   1490:        if ( psn == pslen ) {
                   1491:                pslen *= 2;
                   1492:                ps = (DP *)REALLOC((char *)ps,pslen*sizeof(DP));
                   1493:                psh = (DL *)REALLOC((char *)psh,pslen*sizeof(DL));
                   1494:                pss = (int *)REALLOC((char *)pss,pslen*sizeof(int));
                   1495:                psc = (P *)REALLOC((char *)psc,pslen*sizeof(P)); /* XXX */
                   1496:        }
                   1497:        ps[psn] = a;
                   1498:        psh[psn] = BDY(ps[psn])->dl;
                   1499:        pss[psn] = ps[psn]->sugar;
                   1500:        return psn++;
                   1501: }
                   1502:
1.37      noro     1503: void reducebase_dehomo(NODE f,NODE *g)
1.1       noro     1504: {
                   1505:        int n,i,j,k;
                   1506:        int *r;
                   1507:        DL *w,d;
                   1508:        DP u;
                   1509:        NODE t,top;
                   1510:
                   1511:        n = length(f);
                   1512:        w = (DL *)ALLOCA(n*sizeof(DL));
                   1513:        r = (int *)ALLOCA(n*sizeof(int));
                   1514:        for ( i = 0, t = f; i < n; i++, t = NEXT(t) ) {
                   1515:                r[i] = (int)BDY(t); w[i] = psh[r[i]];
                   1516:        }
                   1517:        for ( i = 0; i < n; i++ ) {
                   1518:                for ( j = 0, d = w[i]; j < n; j++ ) {
                   1519:                        if ( j != i ) {
                   1520:                                for ( k = 0; k < NVars; k++ )
                   1521:                                        if ( d->d[k] < w[j]->d[k] )
                   1522:                                                break;
                   1523:                                if ( k == NVars )
                   1524:                                        break;
                   1525:                        }
                   1526:                }
                   1527:                if ( j != n )
                   1528:                        r[i] = -1;
                   1529:        }
                   1530:        for ( top = 0, i = n-1; i >= 0; i-- )
                   1531:                if ( r[i] >= 0 ) {
1.7       noro     1532:                        dp_load(r[i],&ps[r[i]]); dp_dehomo(ps[r[i]],&u);
                   1533:                        if ( GenTrace ) {
                   1534:                                Q q;
                   1535:                                LIST hist;
                   1536:                                NODE node;
                   1537:
                   1538:                                STOQ(r[i],q);
                   1539:                                node = mknode(4,0,q,0,0);
                   1540:                                MKLIST(hist,node);
                   1541:                                MKNODE(TraceList,hist,0);
                   1542:                        }
                   1543:                        j = newps(u,0,0);
1.1       noro     1544:                        MKNODE(t,j,top); top = t;
                   1545:                }
                   1546:        *g = top;
                   1547: }
                   1548:
1.37      noro     1549: NODE append_one(NODE f,int n)
1.1       noro     1550: {
                   1551:        NODE t;
                   1552:
                   1553:        if ( Reverse || !f ) {
                   1554:                MKNODE(t,(pointer)n,f); return t;
                   1555:        } else {
                   1556:                for ( t = f; NEXT(t); t = NEXT(t) );
                   1557:                MKNODE(NEXT(t),(pointer)n,0);
                   1558:                return f;
                   1559:        }
                   1560: }
                   1561:
1.37      noro     1562: DP_pairs minp( DP_pairs d, DP_pairs *prest )
1.1       noro     1563: {
                   1564:        register DP_pairs m, ml, p, l;
                   1565:        register DL lcm;
                   1566:        register int s, nv = CNVars;
                   1567:        register int (*cmpfun)() = cmpdl;
                   1568:
                   1569:        if ( !(p = NEXT(m = d)) ) {
                   1570:                *prest = p;
                   1571:                NEXT(m) = 0;
                   1572:                return m;
                   1573:        }
                   1574:        for ( lcm = m->lcm, s = m->sugar, ml = 0, l = m; p; p = NEXT(l = p) )
                   1575:                if ( NoSugar ? (*cmpfun)( nv, lcm, p->lcm ) >= 0 :
                   1576:                     (s > p->sugar || s == p->sugar && (*cmpfun)( nv, lcm, p->lcm ) >= 0) )
                   1577:                        ml = l,  lcm = (m = p)->lcm,  s = p->sugar;
                   1578:        if ( !ml ) *prest = NEXT(m);
                   1579:        else {
                   1580:                NEXT(ml) = NEXT(m);
                   1581:                *prest = d;
                   1582:        }
                   1583:        NEXT(m) = 0;
                   1584:        return m;
                   1585: }
                   1586:
1.37      noro     1587: void minsugar(DP_pairs d,DP_pairs *dm,DP_pairs *dr)
1.1       noro     1588: {
                   1589:        int msugar;
                   1590:        DP_pairs t,dm0,dr0,dmt,drt;
                   1591:
                   1592:        for ( msugar = d->sugar, t = NEXT(d); t; t = NEXT(t) )
                   1593:                if ( t->sugar < msugar )
                   1594:                        msugar = t->sugar;
                   1595:        dm0 = 0; dr0 = 0;
                   1596:        for ( t = d; t; t = NEXT(t) ) {
                   1597:                if ( t->sugar == msugar ) {
                   1598:                        NEXTDPP(dm0,dmt);
                   1599:                        dmt->dp1 = t->dp1; dmt->dp2 = t->dp2;
                   1600:                        dmt->lcm = t->lcm; dmt->sugar = t->sugar;
                   1601:                } else {
                   1602:                        NEXTDPP(dr0,drt);
                   1603:                        drt->dp1 = t->dp1; drt->dp2 = t->dp2;
                   1604:                        drt->lcm = t->lcm; drt->sugar = t->sugar;
                   1605:                }
                   1606:        }
                   1607:        if ( dm0 ) NEXT(dmt) = 0;
                   1608:        if ( dr0 ) NEXT(drt) = 0;
                   1609:        *dm = dm0; *dr = dr0;
                   1610: }
                   1611:
1.37      noro     1612: NODE gb(NODE f,int m,NODE subst)
1.1       noro     1613: {
1.42      noro     1614:        int i,nh,prev,mag,mag0,magt;
1.1       noro     1615:        NODE r,g,gall;
1.37      noro     1616:        DP_pairs d;
1.1       noro     1617:        DP_pairs l;
                   1618:        DP h,nf,nfm,dp1,dp2;
                   1619:        MP mp;
1.37      noro     1620:        struct oEGT tnf0,tnf1,tnfm0,tnfm1,tpz0,tpz1,tnp0,tnp1;
1.1       noro     1621:        int skip_nf_flag;
                   1622:        double t_0;
1.7       noro     1623:        Q q;
1.6       noro     1624:        int new_sugar;
1.1       noro     1625:        static prev_sugar = -1;
                   1626:
                   1627:        Max_mag = 0;
1.42      noro     1628:        Max_coef = 0;
1.1       noro     1629:        prev = 1;
                   1630:        doing_f4 = 0;
                   1631:        if ( m ) {
                   1632:                psm = (DP *)MALLOC(pslen*sizeof(DP));
                   1633:                for ( i = 0; i < psn; i++ )
                   1634:                        if ( psh[i] && !validhc(psc[i],m,subst) )
                   1635:                                return 0;
                   1636:                        else
                   1637:                                _dp_mod(ps[i],m,subst,&psm[i]);
                   1638:        }
                   1639:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                   1640:                i = (int)BDY(r);
                   1641:                d = updpairs(d,g,i);
                   1642:                g = updbase(g,i);
                   1643:                gall = append_one(gall,i);
                   1644:        }
                   1645:        while ( d ) {
1.14      noro     1646:                l = minp(d,&d);
1.1       noro     1647:                if ( m ) {
1.6       noro     1648:                        _dp_sp_mod_dup(psm[l->dp1],psm[l->dp2],m,&h);
1.10      noro     1649:                        if ( h )
                   1650:                                new_sugar = h->sugar;
1.1       noro     1651:                        get_eg(&tnfm0);
1.6       noro     1652:                        _dp_nf_mod_destructive(gall,h,psm,m,0,&nfm);
1.1       noro     1653:                        get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
                   1654:                } else
                   1655:                        nfm = (DP)1;
                   1656:                if ( nfm ) {
                   1657:                        if ( Demand ) {
                   1658:                                if ( dp_load_t(psn,&nf) ) {
                   1659:                                        skip_nf_flag = 1;
                   1660:                                        goto skip_nf;
                   1661:                                } else {
                   1662:                                        skip_nf_flag = 0;
                   1663:                                        dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
                   1664:                                        dp_sp(dp1,dp2,&h);
                   1665:                                }
                   1666:                        } else
                   1667:                                dp_sp(ps[l->dp1],ps[l->dp2],&h);
1.7       noro     1668:                        if ( GenTrace ) {
                   1669:                                STOQ(l->dp1,q); ARG1(BDY((LIST)BDY(NEXT(TraceList)))) = q;
                   1670:                                STOQ(l->dp2,q); ARG1(BDY((LIST)BDY(TraceList))) = q;
                   1671:                        }
1.10      noro     1672:                        if ( h )
                   1673:                                new_sugar = h->sugar;
1.1       noro     1674:                        get_eg(&tnf0);
                   1675:                        t_0 = get_rtime();
1.18      noro     1676:                        if ( PCoeffs || dp_fcoeffs )
1.16      noro     1677:                                _dp_nf(gall,h,ps,!Top,&nf);
                   1678:                        else
1.21      noro     1679:                                _dp_nf_z(gall,h,ps,!Top,DP_Multiple,&nf);
1.50      noro     1680:                        if ( DP_Print && nf )
1.1       noro     1681:                                fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
                   1682:                        get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
                   1683:                } else
                   1684:                        nf = 0;
                   1685: skip_nf:
                   1686:                if ( nf ) {
                   1687:                        NZR++;
                   1688:                        get_eg(&tpz0);
                   1689:                        prim_part(nf,0,&h);
                   1690:                        get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
                   1691:                        get_eg(&tnp0);
                   1692:                        if ( Demand && skip_nf_flag )
                   1693:                                nh = newps_nosave(h,m,subst);
                   1694:                        else
                   1695:                                nh = newps(h,m,subst);
                   1696:                        get_eg(&tnp1); add_eg(&eg_np,&tnp0,&tnp1);
                   1697:                        if ( nh < 0 )
                   1698:                                return 0;
                   1699:                        d = updpairs(d,g,nh);
                   1700:                        g = updbase(g,nh);
                   1701:                        gall = append_one(gall,nh);
                   1702:                        if ( !dp_fcoeffs && ShowMag ) {
1.42      noro     1703:                                for ( mag = 0, mag0 = 0, mp = BDY(h); mp; mp = NEXT(mp) ) {
                   1704:                                        magt = p_mag((P)mp->c);
                   1705:                                        mag0 = MAX(mag0,magt);
                   1706:                                        mag += magt;
                   1707:                                }
                   1708:                                Max_coef = MAX(Max_coef,mag0);
1.1       noro     1709:                                Max_mag = MAX(Max_mag,mag);
                   1710:                        }
1.12      noro     1711:                        if ( DP_Print ) {
1.1       noro     1712:                                if ( !prev )
                   1713:                                        fprintf(asir_out,"\n");
                   1714:                                print_split_e(&tnf0,&tnf1); print_split_e(&tpz0,&tpz1);
                   1715:                                printdl(psh[nh]);
                   1716:                                fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",
                   1717:                                        l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
                   1718:                                        pss[nh]);
                   1719:                                if ( ShowMag )
1.42      noro     1720:                                        fprintf(asir_out,",mag=(%d,%d)",mag,mag0);
1.1       noro     1721:                                fprintf(asir_out,"\n"); fflush(asir_out);
1.12      noro     1722:                        } else if ( DP_PrintShort ) {
1.1       noro     1723:                                fprintf(asir_out,"+"); fflush(asir_out);
                   1724:                        }
                   1725:                        prev = 1;
                   1726:                } else {
                   1727:                        if ( m )
                   1728:                                add_eg(&eg_znfm,&tnfm0,&tnfm1);
                   1729:                        ZR++;
1.12      noro     1730:                        if ( DP_Print || DP_PrintShort ) {
1.6       noro     1731:                                if ( new_sugar != prev_sugar ) {
                   1732:                                        fprintf(asir_out,"[%d]",new_sugar);
                   1733:                                        prev_sugar = new_sugar;
1.1       noro     1734:                                }
                   1735:                                fprintf(asir_out,"."); fflush(asir_out); prev = 0;
                   1736:                        }
                   1737:                }
                   1738:        }
1.12      noro     1739:        if ( DP_Print || DP_PrintShort )
1.1       noro     1740:                fprintf(asir_out,"gb done\n");
                   1741:        return g;
                   1742: }
                   1743:
1.37      noro     1744: NODE gb_mod(NODE f,int m)
1.1       noro     1745: {
                   1746:        int i,nh,prev;
                   1747:        NODE r,g,gall;
1.37      noro     1748:        DP_pairs d;
1.1       noro     1749:        DP_pairs l;
                   1750:        DP h,nf;
1.37      noro     1751:        struct oEGT tnfm0,tnfm1,tpz0,tpz1;
1.1       noro     1752:
                   1753:        prev = 1;
                   1754:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                   1755:                i = (int)BDY(r);
                   1756:                d = updpairs(d,g,i);
                   1757:                g = updbase(g,i);
                   1758:                gall = append_one(gall,i);
                   1759:        }
                   1760:        while ( d ) {
1.14      noro     1761:                l = minp(d,&d);
1.1       noro     1762:                if ( PCoeffs ) {
                   1763:                        dp_sp_mod(ps[l->dp1],ps[l->dp2],m,&h);
1.14      noro     1764:                        get_eg(&tnfm0);
1.1       noro     1765:                        dp_nf_mod(gall,h,ps,m,!Top,&nf);
                   1766:                } else {
1.6       noro     1767:                        _dp_sp_mod_dup(ps[l->dp1],ps[l->dp2],m,&h);
1.14      noro     1768:                        get_eg(&tnfm0);
1.6       noro     1769:                        _dp_nf_mod_destructive(gall,h,ps,m,!Top,&nf);
1.1       noro     1770:                }
                   1771:                get_eg(&tnfm1); add_eg(&eg_nfm,&tnfm0,&tnfm1);
                   1772:                if ( nf ) {
                   1773:                        NZR++;
                   1774:                        get_eg(&tpz0);
                   1775:                        prim_part(nf,m,&h);
                   1776:                        get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
                   1777:                        nh = newps_mod(h,m);
                   1778:                        if ( nh < 0 )
                   1779:                                return 0;
                   1780:                        d = updpairs(d,g,nh);
                   1781:                        g = updbase(g,nh);
                   1782:                        gall = append_one(gall,nh);
1.12      noro     1783:                        if ( DP_Print ) {
1.1       noro     1784:                                if ( !prev )
                   1785:                                        fprintf(asir_out,"\n");
                   1786:                                print_split_eg(&tnfm0,&tnfm1); fflush(asir_out);
                   1787:                                fprintf(asir_out,"(%d,%d),nb=%d,nab=%d,rp=%d,sugar=%d",l->dp1,l->dp2,length(g),length(gall),DPPlength(d),pss[nh]);
                   1788:                                printdl(psh[nh]); fprintf(asir_out,"\n"); fflush(asir_out);
1.12      noro     1789:                        } else if ( DP_PrintShort ) {
1.6       noro     1790:                                fprintf(asir_out,"+"); fflush(asir_out);
1.1       noro     1791:                        }
                   1792:                        prev = 1;
                   1793:                } else {
                   1794:                        add_eg(&eg_znfm,&tnfm0,&tnfm1);
                   1795:                        ZR++;
1.12      noro     1796:                        if ( DP_Print || DP_PrintShort ) {
1.1       noro     1797:                                fprintf(asir_out,"."); fflush(asir_out); prev = 0;
                   1798:                        }
                   1799:                }
                   1800:        }
1.12      noro     1801:        if ( DP_Print || DP_PrintShort )
1.1       noro     1802:                fprintf(asir_out,"gb_mod done\n");
                   1803:        return g;
                   1804: }
                   1805:
1.37      noro     1806: DP_pairs updpairs( DP_pairs d, NODE /* of index */ g, int t)
1.1       noro     1807: {
1.14      noro     1808:        register DP_pairs d1, dd, nd;
                   1809:        int dl,dl1;
1.1       noro     1810:
1.14      noro     1811:        if ( !g ) return d;
                   1812:        if ( !NoCriB && d ) {
                   1813:                dl = DPPlength(d);
                   1814:                d = criterion_B( d, t );
                   1815:                dl -= DPPlength(d); NBP += dl;
1.1       noro     1816:        }
1.14      noro     1817:        d1 = newpairs( g, t );
                   1818:        if ( NEXT(d1) ) {
                   1819:                dl = DPPlength(d1); TP += dl;
                   1820:                d1 = criterion_M( d1 );
                   1821:                dl1 = DPPlength(d1); NMP += (dl-dl1); dl = dl1;
                   1822:                d1 = criterion_F( d1 );
                   1823:                dl1 = DPPlength(d1); NFP += (dl-dl1); dl = dl1;
                   1824:        } else
                   1825:                dl = 1;
                   1826:        if ( !do_weyl )
                   1827:                for ( dd = 0; d1; d1 = nd ) {
                   1828:                        nd = NEXT(d1);
                   1829:                        if ( !criterion_2( d1->dp1, d1->dp2 ) ) {
                   1830:                                NEXT(d1) = dd;
                   1831:                                dd = d1;
1.1       noro     1832:                        }
                   1833:                }
1.14      noro     1834:        else
                   1835:                dd = d1;
                   1836:        dl1 = DPPlength(dd); NDP += (dl-dl1);
                   1837:        if ( !(nd = d) ) return dd;
                   1838:        while ( nd = NEXT(d1 = nd) ) ;
                   1839:        NEXT(d1) = dd;
                   1840:        return d;
1.1       noro     1841: }
                   1842:
1.37      noro     1843: DP_pairs newpairs( NODE /* of index */ g, int t )
1.14      noro     1844: {
                   1845:        register NODE r;
                   1846:        register DL tdl = psh[t];
                   1847:        register int ts;
                   1848:        register DP_pairs p, last;
                   1849:        int dp;
                   1850:        register DL dl;
                   1851:        register int s;
1.1       noro     1852:
1.14      noro     1853:        ts = pss[t] - tdl->td;
                   1854:        for ( last = 0, r = g; r; r = NEXT(r) ) {
                   1855:                NEXT(p = NEWDP_pairs) = last;
                   1856:                last = p;
                   1857:                dp = p->dp1 = (int)BDY(r);  p->dp2 = t;
                   1858:                p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
                   1859: #if 0
                   1860:                if ( do_weyl )
1.41      noro     1861:                        p->sugar = dl_weyl_weight(p->lcm);
1.14      noro     1862:                else
                   1863: #endif
                   1864:                        p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
                   1865:        }
                   1866:        return last;
                   1867: }
1.1       noro     1868:
1.37      noro     1869: DP_pairs criterion_B( DP_pairs d, int s )
1.1       noro     1870: {
1.14      noro     1871:        register DP_pairs dd, p;
                   1872:        register DL tij, t = psh[s], dltmp;
1.1       noro     1873:
1.14      noro     1874:        if ( !d ) return 0;
                   1875:        NEWDL( dltmp, CNVars );
                   1876:        for ( dd = 0; d; d = p ) {
                   1877:                p = NEXT(d),
                   1878:                tij = d->lcm;
                   1879:                if ( tij->td != lcm_of_DL(CNVars, tij, t, dltmp )->td
                   1880:                        || !dl_equal(CNVars, tij, dltmp )
                   1881:                        || (tij->td == lcm_of_DL(CNVars, psh[d->dp1], t, dltmp )->td
                   1882:                            && dl_equal(CNVars, dltmp, tij ))
                   1883:                        || (tij->td == lcm_of_DL(CNVars, psh[d->dp2], t, dltmp )->td
                   1884:                            && dl_equal(CNVars, dltmp, tij )) ) {
                   1885:                        NEXT(d) = dd;
                   1886:                        dd = d;
                   1887:                }
1.1       noro     1888:        }
1.14      noro     1889:        return dd;
                   1890: }
1.1       noro     1891:
1.37      noro     1892: DP_pairs criterion_M( DP_pairs d1 )
1.14      noro     1893: {
                   1894:        register DP_pairs dd, e, d3, d2, p;
                   1895:        register DL itdl, jtdl;
                   1896:        register int itdltd, jtdltd;
1.1       noro     1897:
1.14      noro     1898:        for ( dd = 0, e = d1; e; e = d3 ) {
                   1899:                if ( !(d2 = NEXT(e)) ) {
                   1900:                        NEXT(e) = dd;
                   1901:                        return e;
                   1902:                }
                   1903:                itdltd = (itdl = e->lcm)->td;
                   1904:                for ( d3 = 0; d2; d2 = p ) {
                   1905:                        p = NEXT(d2),
                   1906:                        jtdltd = (jtdl = d2->lcm)->td;
                   1907:                        if ( jtdltd == itdltd  )
                   1908:                                if ( dl_equal(CNVars, itdl, jtdl ) ) ;
                   1909:                                else if ( dl_redble( jtdl, itdl ) ) continue;
                   1910:                                else if ( dl_redble( itdl, jtdl ) ) goto delit;
                   1911:                                else ;
                   1912:                        else if ( jtdltd > itdltd )
                   1913:                                if ( dl_redble( jtdl, itdl ) ) continue;
1.1       noro     1914:                                else ;
                   1915:                        else if ( dl_redble( itdl, jtdl ) ) goto delit;
                   1916:                        NEXT(d2) = d3;
                   1917:                        d3 = d2;
                   1918:                }
                   1919:                NEXT(e) = dd;
                   1920:                dd = e;
                   1921:                continue;
                   1922:                /**/
                   1923:        delit:  NEXT(d2) = d3;
                   1924:                d3 = d2;
                   1925:                for ( ; p; p = d2 ) {
                   1926:                        d2 = NEXT(p);
                   1927:                        NEXT(p) = d3;
                   1928:                        d3 = p;
                   1929:                }
                   1930:        }
                   1931:        return dd;
                   1932: }
                   1933:
1.37      noro     1934: static DP_pairs collect_pairs_of_hdlcm( DP_pairs d1, DP_pairs *prest )
1.1       noro     1935: {
                   1936:        register DP_pairs w, p, r, s;
                   1937:        register DL ti;
                   1938:        register int td;
                   1939:
                   1940:        td = (ti = (w = d1)->lcm)->td;
                   1941:        s = NEXT(w);
                   1942:        NEXT(w) = 0;
                   1943:        for ( r = 0; s; s = p ) {
                   1944:                p = NEXT(s);
                   1945:                if ( td == s->lcm->td && dl_equal(CNVars, ti, s->lcm ) )
                   1946:                {
                   1947:                        NEXT(s) = w;
                   1948:                        w = s;
                   1949:                } else {
                   1950:                        NEXT(s) = r;
                   1951:                        r = s;
                   1952:                }
                   1953:        }
                   1954:        *prest = r;
                   1955:        return w;
                   1956: }
                   1957:
1.37      noro     1958: int criterion_2( int dp1, int dp2 )
1.1       noro     1959: {
                   1960:        register int i, *d1, *d2;
                   1961:
                   1962:        d1 = psh[dp1]->d,  d2 = psh[dp2]->d;
                   1963:        for ( i = CNVars; --i >= 0; d1++, d2++ )
                   1964:                if ( (*d1 <= *d2 ? *d1 : *d2) > 0  ) return 0;
                   1965:        return 1;
                   1966: }
                   1967:
1.37      noro     1968: DP_pairs criterion_F( DP_pairs d1 )
1.1       noro     1969: {
                   1970:        DP_pairs rest, head;
                   1971:        register DP_pairs last, p, r, w;
                   1972:        register int s;
                   1973:
                   1974:        for ( head = last = 0, p = d1; NEXT(p); ) {
                   1975:                s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
                   1976:                while ( w = NEXT(w) )
1.52      noro     1977:                        if ( !do_weyl && criterion_2( w->dp1, w->dp2 ) ) {
1.1       noro     1978:                                r = w;
                   1979:                                break;
                   1980:                        } else if ( w->sugar < s ) s = (r = w)->sugar;
                   1981:                if ( last ) NEXT(last) = r;
                   1982:                else head = r;
                   1983:                NEXT(last = r) = 0;
                   1984:                if ( !(p = rest) ) return head;
                   1985:        }
                   1986:        if ( !last ) return p;
                   1987:        NEXT(last) = p;
                   1988:        return head;
                   1989: }
                   1990:
1.37      noro     1991: NODE updbase(NODE g,int t)
1.1       noro     1992: {
                   1993:        g = remove_reducibles(g,t);
                   1994:        g = append_one(g,t);
                   1995:        return g;
                   1996: }
                   1997:
1.37      noro     1998: NODE /* of index */ remove_reducibles(NODE /* of index */ nd, int newdp )
1.1       noro     1999: {
                   2000:        register DL dl, dln;
                   2001:        register NODE last, p, head;
                   2002:        register int td;
                   2003:
                   2004:        dl = psh[newdp];
                   2005:        td = dl->td;
                   2006:        for ( head = last = 0, p = nd; p; ) {
                   2007:                dln = psh[(int)BDY(p)];
                   2008:                if ( dln->td >= td && dl_redble( dln, dl ) ) {
                   2009:                        p = NEXT(p);
                   2010:                        if ( last ) NEXT(last) = p;
                   2011:                } else {
                   2012:                        if ( !last ) head = p;
                   2013:                        p = NEXT(last = p);
                   2014:                }
                   2015:        }
                   2016:        return head;
                   2017: }
                   2018:
1.37      noro     2019: int dl_redble(DL dl1,DL dl2)
1.1       noro     2020: {
                   2021:        register int n, *d1, *d2;
                   2022:
                   2023:        for ( d1 = dl1->d, d2 = dl2->d, n = CNVars; --n >= 0; d1++, d2++ )
                   2024:                if ( *d1 < *d2 ) return 0;
                   2025:        return 1;
1.5       noro     2026: }
                   2027:
1.41      noro     2028: #if 0
                   2029: int dl_weyl_weight(DL dl)
1.5       noro     2030: {
                   2031:        int n,w,i;
                   2032:
                   2033:        n = CNVars/2;
                   2034:        for ( i = 0, w = 0; i < n; i++ )
                   2035:                w += (-dl->d[i]+dl->d[n+i]);
                   2036:        return w;
1.1       noro     2037: }
1.41      noro     2038: #endif
1.1       noro     2039:
1.37      noro     2040: int gbcheck(NODE f)
1.1       noro     2041: {
                   2042:        int i;
                   2043:        NODE r,g,gall;
                   2044:        DP_pairs d,l;
                   2045:        DP h,nf,dp1,dp2;
                   2046:        struct oEGT tmp0,tmp1;
                   2047:
                   2048:        if ( NoGC )
                   2049:                return 1;
                   2050:        for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
                   2051:                i = (int)BDY(r);
                   2052:                d = updpairs(d,g,i);
                   2053:                g = updbase(g,i);
                   2054:                gall = append_one(gall,i);
                   2055:        }
1.12      noro     2056:        if ( DP_Print || DP_PrintShort ) {
1.1       noro     2057:                fprintf(asir_out,"gbcheck total %d pairs\n",DPPlength(d)); fflush(asir_out);
                   2058:        }
                   2059:        while ( d ) {
                   2060:                l = d; d = NEXT(d);
                   2061:                get_eg(&tmp0);
1.17      noro     2062:                dp_load(l->dp1,&dp1); dp_load(l->dp2,&dp2);
                   2063:                dp_sp(dp1,dp2,&h);
                   2064: /* fprintf(stderr,"{%d,%d}",l->dp1,l->dp2); */
1.16      noro     2065:                _dp_nf(gall,h,ps,1,&nf);
1.1       noro     2066:                get_eg(&tmp1); add_eg(&eg_gc,&tmp0,&tmp1);
1.12      noro     2067:                if ( DP_Print || DP_PrintShort ) {
1.1       noro     2068:                        fprintf(asir_out,"."); fflush(asir_out);
                   2069:                }
                   2070:                if ( nf )
                   2071:                        return 0;
                   2072:        }
1.12      noro     2073:        if ( DP_Print || DP_PrintShort )
1.1       noro     2074:                fprintf(asir_out,"\n");
                   2075:        return 1;
1.38      noro     2076: }
                   2077:
1.40      noro     2078: void gbcheck_list(NODE f,int n,VECT *gp,LIST *pp)
1.38      noro     2079: {
                   2080:        int i;
                   2081:        NODE r,g,gall,u,u0,t;
1.39      noro     2082:        VECT vect;
1.38      noro     2083:        LIST pair;
                   2084:        DP_pairs d,l;
                   2085:        Q q1,q2;
                   2086:
1.40      noro     2087:        /* we need the following settings */
                   2088:        NVars = CNVars = n;
1.39      noro     2089:        setup_arrays(f,0,&r);
                   2090:        for ( gall = g = 0, d = 0; r; r = NEXT(r) ) {
1.38      noro     2091:                i = (int)BDY(r);
                   2092:                d = updpairs(d,g,i);
                   2093:                g = updbase(g,i);
                   2094:                gall = append_one(gall,i);
                   2095:        }
1.39      noro     2096:        NEWVECT(vect); vect->len = psn; vect->body = (pointer)ps;
                   2097:        *gp = vect;
                   2098:
1.38      noro     2099:        for ( u0 = 0, l = d; l; l = NEXT(l) ) {
                   2100:                NEXTNODE(u0,u);
                   2101:                STOQ(l->dp1,q1);
                   2102:                STOQ(l->dp2,q2);
                   2103:                t = mknode(2,q1,q2);
                   2104:                MKLIST(pair,t);
                   2105:                BDY(u) = (pointer)pair;
                   2106:        }
                   2107:        if ( u0 )
                   2108:                NEXT(u) = 0;
1.39      noro     2109:        MKLIST(*pp,u0);
1.1       noro     2110: }
                   2111:
1.37      noro     2112: int membercheck(NODE f,NODE x)
1.1       noro     2113: {
                   2114:        DP g;
                   2115:        struct oEGT tmp0,tmp1;
                   2116:
                   2117:        if ( NoMC )
                   2118:                return 1;
1.12      noro     2119:        if ( DP_Print || DP_PrintShort ) {
1.1       noro     2120:                fprintf(asir_out,"membercheck\n"); fflush(asir_out);
                   2121:        }
                   2122:        for ( ; f; f = NEXT(f) ) {
                   2123:                get_eg(&tmp0);
1.16      noro     2124:                _dp_nf(x,(DP)BDY(f),ps,1,&g);
1.1       noro     2125:                get_eg(&tmp1); add_eg(&eg_mc,&tmp0,&tmp1);
1.12      noro     2126:                if ( DP_Print ) {
1.1       noro     2127:                        print_split_eg(&tmp0,&tmp1); fflush(asir_out);
1.12      noro     2128:                } else if ( DP_PrintShort ) {
1.1       noro     2129:                        fprintf(asir_out,"."); fflush(asir_out);
                   2130:                }
                   2131:                if ( g )
                   2132:                        return 0;
                   2133:        }
1.12      noro     2134:        if ( DP_Print || DP_PrintShort )
1.1       noro     2135:                fprintf(asir_out,"\n");
                   2136:        return 1;
                   2137: }
                   2138:
1.37      noro     2139: void dp_set_flag(Obj name,Obj value)
1.1       noro     2140: {
                   2141:        char *n;
                   2142:        int v;
1.43      noro     2143:        Q ratio;
1.1       noro     2144:
                   2145:        if ( OID(name) != O_STR )
                   2146:                return;
                   2147:        n = BDY((STRING)name);
                   2148:        if ( !strcmp(n,"Demand") ) {
                   2149:                Demand = value ? BDY((STRING)value) : 0; return;
                   2150:        }
                   2151:        if ( !strcmp(n,"Dist") ) {
                   2152:                Dist = (LIST)value; return;
                   2153:        }
1.43      noro     2154:        if ( !strcmp(n,"Content") ) {
                   2155:                ratio = (Q)value;
                   2156:                if ( ratio ) {
                   2157:                        DP_Multiple = BD(NM(ratio))[0];
                   2158:                        Denominator = INT(ratio) ? 1 : BD(DN(ratio))[0];
                   2159:                } else {
                   2160:                        DP_Multiple = 0;
                   2161:                        Denominator = 1;
                   2162:                }
                   2163:        }
1.1       noro     2164:        if ( value && OID(value) != O_N )
                   2165:                return;
                   2166:        v = QTOS((Q)value);
                   2167:        if ( !strcmp(n,"NoSugar") )
                   2168:                NoSugar = v;
                   2169:        else if ( !strcmp(n,"NoCriB") )
                   2170:                NoCriB = v;
                   2171:        else if ( !strcmp(n,"NoGC") )
                   2172:                NoGC = v;
                   2173:        else if ( !strcmp(n,"NoMC") )
                   2174:                NoMC = v;
                   2175:        else if ( !strcmp(n,"NoRA") )
                   2176:                NoRA = v;
                   2177:        else if ( !strcmp(n,"NoGCD") )
                   2178:                NoGCD = v;
                   2179:        else if ( !strcmp(n,"Top") )
                   2180:                Top = v;
                   2181:        else if ( !strcmp(n,"ShowMag") )
                   2182:                ShowMag = v;
1.14      noro     2183:        else if ( !strcmp(n,"PrintShort") )
1.12      noro     2184:                DP_PrintShort = v;
1.14      noro     2185:        else if ( !strcmp(n,"Print") )
1.12      noro     2186:                DP_Print = v;
1.14      noro     2187:        else if ( !strcmp(n,"NFStat") )
                   2188:                DP_NFStat = v;
1.1       noro     2189:        else if ( !strcmp(n,"Stat") )
                   2190:                Stat = v;
                   2191:        else if ( !strcmp(n,"Reverse") )
                   2192:                Reverse = v;
                   2193:        else if ( !strcmp(n,"Multiple") )
1.14      noro     2194:                DP_Multiple = v;
1.1       noro     2195:        else if ( !strcmp(n,"Denominator") )
                   2196:                Denominator = v;
                   2197:        else if ( !strcmp(n,"PtozpRA") )
                   2198:                PtozpRA = v;
1.7       noro     2199:        else if ( !strcmp(n,"GenTrace") )
                   2200:                GenTrace = v;
                   2201:        else if ( !strcmp(n,"OXCheck") )
                   2202:                OXCheck = v;
1.1       noro     2203: }
                   2204:
1.37      noro     2205: void dp_make_flaglist(LIST *list)
1.1       noro     2206: {
1.43      noro     2207:        Q v,nm,dn;
1.1       noro     2208:        STRING name,path;
                   2209:        NODE n,n1;
                   2210:
1.43      noro     2211: #if 0
1.14      noro     2212:        STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1;
1.1       noro     2213:        STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
1.43      noro     2214: #else
                   2215:        if ( DP_Multiple ) {
                   2216:                STOQ(DP_Multiple,nm); STOQ(Denominator,dn); divq(nm,dn,&v);
                   2217:        } else
                   2218:                v = 0;
                   2219:        MKNODE(n,v,0); MKSTR(name,"Content"); MKNODE(n1,name,n); n = n1;
                   2220: #endif
1.1       noro     2221:        MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
                   2222:        STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
                   2223:        STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
1.14      noro     2224:        STOQ(DP_Print,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Print"); MKNODE(n1,name,n); n = n1;
1.20      noro     2225:        STOQ(DP_PrintShort,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PrintShort"); MKNODE(n1,name,n); n = n1;
1.14      noro     2226:        STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1;
1.7       noro     2227:        STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1;
                   2228:        STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1;
1.1       noro     2229:        STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
                   2230:        STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
                   2231:        STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
                   2232:        STOQ(NoGCD,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGCD"); MKNODE(n1,name,n); n = n1;
                   2233:        STOQ(NoRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoRA"); MKNODE(n1,name,n); n = n1;
                   2234:        STOQ(NoMC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoMC"); MKNODE(n1,name,n); n = n1;
                   2235:        STOQ(NoGC,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoGC"); MKNODE(n1,name,n); n = n1;
                   2236:        STOQ(NoCriB,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoCriB"); MKNODE(n1,name,n); n = n1;
                   2237:        STOQ(NoSugar,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NoSugar"); MKNODE(n1,name,n); n = n1;
                   2238:        if ( Demand )
                   2239:                MKSTR(path,Demand);
                   2240:        else
                   2241:                path = 0;
                   2242:        MKNODE(n1,path,n); n = n1; MKSTR(name,"Demand"); MKNODE(n1,name,n); n = n1;
                   2243:        MKLIST(*list,n);
                   2244: }
                   2245:
                   2246: #define DELIM '/'
                   2247:
1.37      noro     2248: void dp_save(int index,Obj p,char *prefix)
1.1       noro     2249: {
                   2250:        FILE *fp;
                   2251:        char path[BUFSIZ];
                   2252:
1.7       noro     2253:        if ( prefix )
                   2254:                sprintf(path,"%s%c%s%d",Demand,DELIM,prefix,index);
                   2255:        else
                   2256:                sprintf(path,"%s%c%d",Demand,DELIM,index);
1.1       noro     2257:        if ( !(fp = fopen(path,"wb") ) )
                   2258:                error("dp_save : cannot open a file");
1.7       noro     2259:        savevl(fp,VC); saveobj(fp,p); fclose(fp);
1.1       noro     2260: }
                   2261:
1.37      noro     2262: void dp_load(int index,DP *p)
1.1       noro     2263: {
                   2264:        FILE *fp;
                   2265:        char path[BUFSIZ];
                   2266:
                   2267:        if ( !Demand || ps[index] )
                   2268:                *p = ps[index];
                   2269:        else {
                   2270:                sprintf(path,"%s%c%d",Demand,DELIM,index);
                   2271:                if ( !(fp = fopen(path,"rb") ) )
                   2272:                        error("dp_load : cannot open a file");
1.45      noro     2273:                if ( PCoeffs )
                   2274:                        loadvl(fp);
                   2275:                else
                   2276:                        skipvl(fp);
                   2277:                loadobj(fp,(Obj *)p); fclose(fp);
1.1       noro     2278:        }
                   2279: }
                   2280:
1.37      noro     2281: int dp_load_t(int index,DP *p)
1.1       noro     2282: {
                   2283:        FILE *fp;
                   2284:        char path[BUFSIZ];
                   2285:
                   2286:        sprintf(path,"%s%c%d",Demand,DELIM,index);
                   2287:        if ( !(fp = fopen(path,"rb") ) )
                   2288:                return 0;
                   2289:        else {
1.45      noro     2290:                if ( PCoeffs )
                   2291:                        loadvl(fp);
                   2292:                else
                   2293:                        skipvl(fp);
                   2294:                loadobj(fp,(Obj *)p); fclose(fp); return 1;
1.1       noro     2295:        }
                   2296: }
                   2297:
                   2298: void init_stat() {
                   2299:        init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
1.14      noro     2300:        init_eg(&eg_pz); init_eg(&eg_np);
1.1       noro     2301:        init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
1.35      noro     2302:        ZR = NZR = TP = NMP = NBP = NFP = NDP = 0;
1.1       noro     2303: }
                   2304:
                   2305: void print_stat() {
1.12      noro     2306:        if ( !DP_Print && !Stat )
1.1       noro     2307:                return;
                   2308:        print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
1.14      noro     2309:        print_eg("PZ",&eg_pz); print_eg("NP",&eg_np);
1.1       noro     2310:        print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
                   2311:        fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,NBP,NMP,NFP,NDP,ZR,NZR);
                   2312: }
1.14      noro     2313:
                   2314: /*
                   2315:  * dp_nf used in gb()
                   2316:  *
                   2317:  */
                   2318:
                   2319: double pz_t_e, pz_t_d, pz_t_d1, pz_t_c, im_t_s, im_t_r;
                   2320:
                   2321: extern int GenTrace;
                   2322: extern NODE TraceList;
1.15      noro     2323: extern int mpi_mag;
                   2324:
1.37      noro     2325: void dp_mulc_d(DP p,P c,DP *r)
1.15      noro     2326: {
                   2327:        if ( Dist && BDY(Dist)
                   2328:                && HMAG(p) > mpi_mag
                   2329:                && p_mag((P)c) > mpi_mag ) {
                   2330:                if ( DP_NFStat ) fprintf(asir_out,"~");
                   2331:                dp_imul_d(p,(Q)c,r);
                   2332:        } else {
                   2333:                if ( DP_NFStat ) fprintf(asir_out,"_");
                   2334:                muldc(CO,p,c,r);
                   2335:        }
                   2336: }
1.14      noro     2337:
1.37      noro     2338: void _dp_nf(NODE b,DP g,DP *ps,int full,DP *rp)
1.16      noro     2339: {
                   2340:        DP u,p,d,s,t,mult;
                   2341:        P coef;
                   2342:        NODE l;
                   2343:        MP m,mr;
                   2344:        int sugar,psugar;
                   2345:
                   2346:        if ( !g ) {
                   2347:                *rp = 0; return;
                   2348:        }
                   2349:        sugar = g->sugar;
                   2350:        for ( d = 0; g; ) {
                   2351:                for ( u = 0, l = b; l; l = NEXT(l) ) {
                   2352:                        if ( dl_redble(BDY(g)->dl,psh[(int)BDY(l)]) ) {
                   2353:                                dp_load((int)BDY(l),&p);
                   2354:                                /* t+u = coef*(d+g) - mult*p (t = coef*d) */
                   2355:                                dp_red(d,g,p,&t,&u,&coef,&mult);
                   2356:                                psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
                   2357:                                sugar = MAX(sugar,psugar);
                   2358:                                if ( GenTrace ) {
                   2359:                                        LIST hist;
                   2360:                                        Q cq;
                   2361:                                        NODE node,node0;
                   2362:
                   2363:                                        STOQ((int)BDY(l),cq);
                   2364:                                        node0 = mknode(4,coef,cq,mult,ONE);
                   2365:                                        MKLIST(hist,node0);
                   2366:                                        MKNODE(node,hist,TraceList); TraceList = node;
                   2367:                                }
                   2368:                                if ( !u ) {
                   2369:                                        if ( d )
                   2370:                                                d->sugar = sugar;
                   2371:                                        *rp = d; return;
                   2372:                                }
                   2373:                                d = t;
                   2374:                                break;
                   2375:                        }
                   2376:                }
                   2377:                if ( u )
                   2378:                        g = u;
                   2379:                else if ( !full ) {
                   2380:                        if ( g ) {
                   2381:                                MKDP(g->nv,BDY(g),t); t->sugar = sugar; g = t;
                   2382:                        }
                   2383:                        *rp = g; return;
                   2384:                } else {
                   2385:                        m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
                   2386:                        NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
                   2387:                        addd(CO,d,t,&s); d = s;
                   2388:                        dp_rest(g,&t); g = t;
                   2389:                }
                   2390:        }
                   2391:        if ( d )
                   2392:                d->sugar = sugar;
                   2393:        *rp = d;
                   2394: }
                   2395:
1.37      noro     2396: void _dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *r)
1.14      noro     2397: {
1.37      noro     2398:        DP u,dp,rp,t,t1,red,shift;
1.14      noro     2399:        Q dc,rc,dcq,rcq,cont,hr,hred,cr,cred,mcred,c,gcd,cq;
                   2400:        NODE l;
                   2401:        int hmag,denom;
                   2402:        int sugar,psugar;
                   2403:        STRING imul;
1.37      noro     2404:        double t_0,tt,t_p,t_m,t_g,t_a;
1.14      noro     2405:        LIST hist;
                   2406:        NODE node;
                   2407:        Q rcred,mrcred;
                   2408:
                   2409:        if ( !g ) {
                   2410:                *r = 0; return;
                   2411:        }
                   2412:        pz_t_e = pz_t_d = pz_t_d1 = pz_t_c = 0;
1.15      noro     2413:        t_p = t_m = t_g = t_a = 0;
1.14      noro     2414:
                   2415:        denom = Denominator?Denominator:1;
                   2416:        hmag = multiple*HMAG(g)/denom;
                   2417:        sugar = g->sugar;
                   2418:
                   2419:        dc = 0; dp = 0; rc = ONE; rp = g;
                   2420:        MKSTR(imul,"dp_imul_index");
                   2421:
                   2422:        /* g = dc*dp+rc*rp */
                   2423:        for ( ; rp; ) {
                   2424:                for ( u = 0, l = b; l; l = NEXT(l) ) {
                   2425:                        if ( dl_redble(BDY(rp)->dl,psh[(int)BDY(l)]) ) {
                   2426:                                t_0 = get_rtime();
                   2427:                                dp_load((int)BDY(l),&red);
                   2428:                                hr = (Q)BDY(rp)->c; hred = (Q)BDY(red)->c;
                   2429:                                igcd_cofactor((Q)BDY(rp)->c,(Q)BDY(red)->c,&gcd,&cred,&cr);
1.15      noro     2430:                                tt = get_rtime(); t_p += tt-t_0;
1.14      noro     2431:
                   2432:                                dp_subd(rp,red,&shift);
1.37      noro     2433:                                dp_mulc_d(rp,(P)cr,&t);
1.14      noro     2434:                                chsgnp((P)cred,(P *)&mcred);
1.37      noro     2435:                                dp_mulc_d(red,(P)mcred,&t1);
1.15      noro     2436:                                muld(CO,shift,t1,&t1);
                   2437:                                addd(CO,t,t1,&u);
                   2438:                                t_m += get_rtime()-tt;
1.14      noro     2439:
                   2440:                                psugar = (BDY(rp)->dl->td - BDY(red)->dl->td) + red->sugar;
                   2441:                                sugar = MAX(sugar,psugar);
1.15      noro     2442:
1.14      noro     2443:                                if ( GenTrace ) {
                   2444:                                        /* u = cr*rp + (-cred)*shift*red */
                   2445:                                        STOQ((int)BDY(l),cq);
                   2446:                                        node = mknode(4,cr,cq,0,0);
                   2447:                                        mulq(cred,rc,&rcred);
                   2448:                                        chsgnnum((Num)rcred,(Num *)&mrcred);
                   2449:                                        muldc(CO,shift,(P)mrcred,(DP *)&ARG2(node));
                   2450:                                        MKLIST(hist,node);
                   2451:                                }
1.15      noro     2452:
1.14      noro     2453:                                if ( !u ) {
                   2454:                                        if ( dp )
                   2455:                                                dp->sugar = sugar;
                   2456:                                        *r = dp;
                   2457:                                        if ( GenTrace ) {
                   2458:                                                ARG3(BDY(hist)) = ONE;
                   2459:                                                MKNODE(node,hist,TraceList); TraceList = node;
                   2460:                                        }
                   2461:                                        goto final;
                   2462:                                }
                   2463:                                break;
                   2464:                        }
                   2465:                }
                   2466:                if ( u ) {
                   2467:                        if ( multiple && HMAG(u) > hmag ) {
                   2468:                                t_0 = get_rtime();
1.16      noro     2469:                                dp_ptozp_d(u,&rp);
1.15      noro     2470:                                tt = get_rtime(); t_g += tt-t_0;
                   2471:
                   2472:                                divsq((Q)BDY(u)->c,(Q)BDY(rp)->c,&cont);
1.14      noro     2473:                                if ( !dp_fcoeffs && DP_NFStat ) {
1.15      noro     2474:                                        fprintf(asir_out,
                   2475:                                                "(%d)",p_mag((P)cont)*100/p_mag((P)BDY(u)->c));
1.14      noro     2476:                                        fflush(asir_out);
                   2477:                                }
1.15      noro     2478:                                mulq(cr,dc,&dcq); mulq(cont,rc,&rcq);
1.14      noro     2479:                                igcd_cofactor(dcq,rcq,&gcd,&dc,&rc);
1.15      noro     2480:                                t_a = get_rtime()-tt;
                   2481:
1.14      noro     2482:                                hmag = multiple*HMAG(rp)/denom;
                   2483:                                if ( GenTrace ) {
                   2484:                                        ARG3(BDY(hist)) = (pointer)gcd;
                   2485:                                        MKNODE(node,hist,TraceList); TraceList = node;
                   2486:                                }
                   2487:                        } else {
1.15      noro     2488:                                rp = u;
1.14      noro     2489:                                t_0 = get_rtime();
1.15      noro     2490:                                mulq(cr,dc,&dc);
                   2491:                                t_a += get_rtime()-t_0;
1.14      noro     2492:                                if ( GenTrace ) {
                   2493:                                        ARG3(BDY(hist)) = (pointer)ONE;
                   2494:                                        MKNODE(node,hist,TraceList); TraceList = node;
                   2495:                                }
                   2496:                        }
                   2497:                } else if ( !full ) {
                   2498:                        if ( rp ) {
                   2499:                                MKDP(rp->nv,BDY(rp),t); t->sugar = sugar; rp = t;
                   2500:                        }
                   2501:                        *r = rp;
                   2502:                        goto final;
                   2503:                } else {
                   2504:                        t_0 = get_rtime();
                   2505:                        mulq((Q)BDY(rp)->c,rc,&c);
1.15      noro     2506:                        igcd_cofactor(dc,c,&dc,&dcq,&cq);
                   2507:                        muldc(CO,dp,(P)dcq,&t);
                   2508:                        dp_hm(rp,&t1); BDY(t1)->c = (P)cq;  addd(CO,t,t1,&dp);
                   2509:                        dp_rest(rp,&rp);
                   2510:                        t_a += get_rtime()-t_0;
1.14      noro     2511:                }
                   2512:        }
                   2513:        if ( GenTrace ) {
                   2514:                mulq(ARG3(BDY((LIST)BDY(TraceList))),dc,&cq);
                   2515:                ARG3(BDY((LIST)BDY(TraceList))) = (pointer)cq;
                   2516:        }
                   2517:        if ( dp )
                   2518:                dp->sugar = sugar;
                   2519:        *r = dp;
                   2520: final:
                   2521:        if ( DP_NFStat )
1.15      noro     2522:                fprintf(asir_out,
                   2523:                        "(%.3g %.3g %.3g %.3g %.3g %.3g %.3g %.3g)",
                   2524:                        t_p,t_m,t_g,t_a,
1.14      noro     2525:                        pz_t_e, pz_t_d, pz_t_d1, pz_t_c);
                   2526: }
                   2527:
                   2528: void imulv();
                   2529:
1.37      noro     2530: void dp_imul_d(DP p,Q q,DP *rp)
1.14      noro     2531: {
                   2532:        int nsep,ndist,i,j,k,l,n;
                   2533:        double t0,t1,t2;
                   2534:        Q *s;
                   2535:        pointer *b;
                   2536:        VECT c,cs,ri;
                   2537:        VECT *r;
                   2538:        MP m;
1.37      noro     2539:        NODE tn,dist,n0;
                   2540:        Obj obj;
1.14      noro     2541:        STRING imul;
                   2542:        extern LIST Dist;
                   2543:
                   2544:        if ( !p || !q ) {
                   2545:                *rp = 0; return;
                   2546:        }
                   2547:        dist = BDY(Dist);
                   2548:        for ( tn = dist, ndist = 0; tn; tn = NEXT(tn), ndist++ );
                   2549:        nsep = ndist + 1;
                   2550:        for ( m = BDY(p), n = 0; m; m = NEXT(m), n++ );
                   2551:        if ( n <= nsep ) {
                   2552:                muldc(CO,p,(P)q,rp); return;
                   2553:        }
                   2554:        MKSTR(imul,"imulv");
                   2555:        t0 = get_rtime();
                   2556:        dp_dtov(p,&c);
                   2557:        sepvect(c,nsep,&cs);
                   2558:        r = (VECT *)CALLOC(nsep,sizeof(VECT *));
                   2559:        for ( i = 0, tn = dist, b = BDY(cs); i < ndist; i++, tn = NEXT(tn) ) {
                   2560:                n0 = mknode(4,BDY(tn),imul,b[i],q);
1.37      noro     2561:                Pox_rpc(n0,&obj);
1.14      noro     2562:        }
                   2563:        t1 = get_rtime();
                   2564:        im_t_s += t1 - t0;
                   2565:        imulv(b[i],q,&r[i]);
                   2566:        t1 = get_rtime();
                   2567:        for ( i = 0, tn = dist; i < ndist; i++, tn = NEXT(tn) ) {
                   2568:                MKNODE(n0,BDY(tn),0);
1.37      noro     2569:                Pox_pop_local(n0,&obj); r[i] = (VECT)obj;
1.14      noro     2570:                if ( OID(r[i]) == O_ERR ) {
                   2571:                        printexpr(CO,(Obj)r[i]);
                   2572:                        error("dp_imul_d : aborted");
                   2573:                }
                   2574:        }
                   2575:        t2 = get_rtime();
                   2576:        im_t_r += t2 - t1;
                   2577:        s = (Q *)CALLOC(n,sizeof(Q));
                   2578:        for ( i = j = 0; i < nsep; i++ ) {
                   2579:                for ( k = 0, ri = r[i], l = ri->len; k < l; k++, j++ ) {
                   2580:                        s[j] = (Q)BDY(ri)[k];
                   2581:                }
                   2582:        }
                   2583:        dp_vtod(s,p,rp);
                   2584: }
                   2585:
1.37      noro     2586: void imulv(VECT w,Q c,VECT *rp)
1.14      noro     2587: {
                   2588:        int n,i;
                   2589:        VECT r;
                   2590:
                   2591:        n = w->len;
                   2592:        MKVECT(r,n); *rp = r;
                   2593:        for ( i = 0; i < n; i++ )
                   2594:                mulq((Q)BDY(w)[i],(Q)c,(Q *)&BDY(r)[i]);
1.32      noro     2595: }
                   2596:
1.37      noro     2597: void dptoca(DP p,unsigned int **rp)
1.32      noro     2598: {
                   2599:        int i;
                   2600:        MP m;
                   2601:        unsigned int *r;
                   2602:
                   2603:        if ( !p )
                   2604:                *rp = 0;
                   2605:        else {
                   2606:                for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ );
                   2607:                *rp = r = (unsigned int *)MALLOC_ATOMIC(i*sizeof(unsigned int));
                   2608:                for ( m = BDY(p), i = 0; m; m = NEXT(m), i++ )
                   2609:                        r[i] = ITOS(C(m));
                   2610:        }
1.14      noro     2611: }
                   2612:

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