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

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

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