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