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

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

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