Annotation of OpenXM_contrib2/asir2018/engine/F.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 <math.h>
! 52:
! 53: int use_new_hensel;
! 54:
! 55: void fctrp(VL vl,P f,DCP *dcp)
! 56: {
! 57: VL nvl;
! 58: DCP dc;
! 59:
! 60: if ( !f || NUM(f) ) {
! 61: NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
! 62: NEXT(dc) = 0; *dcp = dc;
! 63: return;
! 64: } else if ( !qpcheck((Obj)f) )
! 65: error("fctrp : invalid argument");
! 66: else {
! 67: clctv(vl,f,&nvl);
! 68: if ( !NEXT(nvl) )
! 69: ufctr(f,1,dcp);
! 70: else
! 71: mfctr(nvl,f,dcp);
! 72: }
! 73: }
! 74:
! 75: void fctr_wrt_v_p(VL vl,P f,V v,DCP *dcp)
! 76: {
! 77: VL nvl;
! 78: DCP dc;
! 79:
! 80: if ( !f || NUM(f) ) {
! 81: NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
! 82: NEXT(dc) = 0; *dcp = dc;
! 83: return;
! 84: } else if ( !qpcheck((Obj)f) )
! 85: error("fctrp : invalid argument");
! 86: else {
! 87: clctv(vl,f,&nvl);
! 88: if ( !NEXT(nvl) )
! 89: ufctr(f,1,dcp);
! 90: else
! 91: mfctr_wrt_v(nvl,f,v,dcp);
! 92: }
! 93: }
! 94:
! 95: void homfctr(VL vl,P g,DCP *dcp)
! 96: {
! 97: P s,t,u,f,h,x;
! 98: Z e;
! 99: int d,d0;
! 100: DCP dc,dct;
! 101:
! 102: substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
! 103: for ( dct = dc; dct; dct = NEXT(dct) )
! 104: if ( !NUM(dc->c) ) {
! 105: for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
! 106: exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
! 107: mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
! 108: subp(vl,f,h,&u); f = u;
! 109: }
! 110: dc->c = s;
! 111: }
! 112: *dcp = dc;
! 113: }
! 114:
! 115: void mfctr(VL vl,P f,DCP *dcp)
! 116: {
! 117: DCP dc,dc0,dct,dcs,dcr;
! 118: P p,pmin,ppmin,cmin,t;
! 119: VL mvl;
! 120: Q c;
! 121:
! 122: ptozp(f,1,&c,&p);
! 123: NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
! 124: msqfr(vl,p,&dct);
! 125: for ( ; dct; dct = NEXT(dct) ) {
! 126: mindegp(vl,COEF(dct),&mvl,&pmin);
! 127: #if 0
! 128: minlcdegp(vl,COEF(dct),&mvl,&pmin);
! 129: min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin);
! 130: #endif
! 131: pcp(mvl,pmin,&ppmin,&cmin);
! 132: if ( !NUM(cmin) ) {
! 133: mfctr(mvl,cmin,&dcs);
! 134: for ( dcr = NEXT(dcs); dcr; dcr = NEXT(dcr) ) {
! 135: DEG(dcr) = DEG(dct);
! 136: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 137: }
! 138: for ( ; NEXT(dc); dc = NEXT(dc) );
! 139: NEXT(dc) = NEXT(dcs);
! 140: }
! 141: mfctrmain(mvl,ppmin,&dcs);
! 142: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 143: DEG(dcr) = DEG(dct);
! 144: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 145: }
! 146: for ( ; NEXT(dc); dc = NEXT(dc) );
! 147: NEXT(dc) = dcs;
! 148: }
! 149: adjsgn(f,dc0); *dcp = dc0;
! 150: }
! 151:
! 152: void mfctr_wrt_v(VL vl,P f,V v,DCP *dcp)
! 153: {
! 154: DCP dc,dc0,dct,dcs,dcr;
! 155: P p,pmin,ppmin,cmin,t;
! 156: VL nvl,mvl;
! 157: Q c;
! 158:
! 159: ptozp(f,1,&c,&p);
! 160: NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
! 161: msqfr(vl,p,&dct);
! 162: for ( ; dct; dct = NEXT(dct) ) {
! 163: clctv(vl,COEF(dct),&nvl);
! 164: reordvar(nvl,v,&mvl);
! 165: reorderp(mvl,vl,COEF(dct),&pmin);
! 166: pcp(mvl,pmin,&ppmin,&cmin);
! 167: if ( !NUM(cmin) ) {
! 168: mfctrmain(mvl,cmin,&dcs);
! 169: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 170: DEG(dcr) = DEG(dct);
! 171: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 172: }
! 173: for ( ; NEXT(dc); dc = NEXT(dc) );
! 174: NEXT(dc) = dcs;
! 175: }
! 176: mfctrmain(mvl,ppmin,&dcs);
! 177: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 178: DEG(dcr) = DEG(dct);
! 179: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 180: }
! 181: for ( ; NEXT(dc); dc = NEXT(dc) );
! 182: NEXT(dc) = dcs;
! 183: }
! 184: adjsgn(f,dc0); *dcp = dc0;
! 185: }
! 186:
! 187: #if 0
! 188: void adjsgn(P p,DCP dc)
! 189: {
! 190: int sgn;
! 191: DCP dct;
! 192: P c;
! 193:
! 194: for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
! 195: if ( !EVENN(NM(dct->d)) )
! 196: sgn *= headsgn(COEF(dct));
! 197: if ( sgn < 0 ) {
! 198: chsgnp(COEF(dc),&c); COEF(dc) = c;
! 199: }
! 200: }
! 201: #else
! 202: void adjsgn(P p,DCP dc)
! 203: {
! 204: DCP dct;
! 205: P c;
! 206:
! 207: if ( headsgn(COEF(dc)) != headsgn(p) ) {
! 208: chsgnp(COEF(dc),&c); COEF(dc) = c;
! 209: }
! 210: for ( dct = NEXT(dc); dct; dct = NEXT(dct) )
! 211: if ( headsgn(COEF(dct)) < 0 ) {
! 212: chsgnp(COEF(dct),&c); COEF(dct) = c;
! 213: }
! 214: }
! 215: #endif
! 216:
! 217: int headsgn(P p)
! 218: {
! 219: if ( !p )
! 220: return 0;
! 221: else {
! 222: while ( !NUM(p) )
! 223: p = COEF(DC(p));
! 224: return sgnq((Q)p);
! 225: }
! 226: }
! 227:
! 228: void fctrwithmvp(VL vl,P f,V v,DCP *dcp)
! 229: {
! 230: VL nvl;
! 231: DCP dc;
! 232:
! 233: if ( NUM(f) ) {
! 234: NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
! 235: NEXT(dc) = 0; *dcp = dc;
! 236: return;
! 237: }
! 238:
! 239: clctv(vl,f,&nvl);
! 240: if ( !NEXT(nvl) )
! 241: ufctr(f,1,dcp);
! 242: else
! 243: mfctrwithmv(nvl,f,v,dcp);
! 244: }
! 245:
! 246: void mfctrwithmv(VL vl,P f,V v,DCP *dcp)
! 247: {
! 248: DCP dc,dc0,dct,dcs,dcr;
! 249: P p,pmin,ppmin,cmin,t;
! 250: VL mvl;
! 251: Q c;
! 252:
! 253: ptozp(f,1,&c,&p);
! 254: NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
! 255: msqfr(vl,p,&dct);
! 256: for ( ; dct; dct = NEXT(dct) ) {
! 257: if ( !v )
! 258: mindegp(vl,COEF(dct),&mvl,&pmin);
! 259: else {
! 260: reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
! 261: }
! 262: pcp(mvl,pmin,&ppmin,&cmin);
! 263: if ( !NUM(cmin) ) {
! 264: mfctrmain(mvl,cmin,&dcs);
! 265: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 266: DEG(dcr) = DEG(dct);
! 267: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 268: }
! 269: for ( ; NEXT(dc); dc = NEXT(dc) );
! 270: NEXT(dc) = dcs;
! 271: }
! 272: mfctrmain(mvl,ppmin,&dcs);
! 273: for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
! 274: DEG(dcr) = DEG(dct);
! 275: reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
! 276: }
! 277: for ( ; NEXT(dc); dc = NEXT(dc) );
! 278: NEXT(dc) = dcs;
! 279: }
! 280: *dcp = dc0;
! 281: }
! 282:
! 283: void ufctr(P f,int hint,DCP *dcp)
! 284: {
! 285: P p,c;
! 286: DCP dc,dct,dcs,dcr;
! 287:
! 288: ptozp(f,sgnq((Q)UCOEF(f)),(Q *)&c,&p);
! 289: NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
! 290: usqp(p,&dct);
! 291: for ( *dcp = dc; dct; dct = NEXT(dct) ) {
! 292: ufctrmain(COEF(dct),hint,&dcs);
! 293: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
! 294: DEG(dcr) = DEG(dct);
! 295: for ( ; NEXT(dc); dc = NEXT(dc) );
! 296: NEXT(dc) = dcs;
! 297: }
! 298: }
! 299:
! 300: void mfctrmain(VL vl,P p,DCP *dcp)
! 301: {
! 302: int i,j,k,*win,np,x;
! 303: VL nvl,tvl;
! 304: VN vn,vnt,vn1;
! 305: P p0,f,g,f0,g0,s,t,lp,m;
! 306: P *fp0,*fpt,*l,*tl;
! 307: DCP dc,dc0,dcl;
! 308: int count,nv;
! 309: int *nonzero;
! 310:
! 311: if ( !cmpz(DEG(DC(p)),ONE) ) {
! 312: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 313: *dcp = dc; return;
! 314: }
! 315: lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
! 316: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
! 317: W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
! 318: W_CALLOC(nv,struct oVN,vn1);
! 319: W_CALLOC(nv,int,nonzero);
! 320:
! 321: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 322: vn1[i].v = vn[i].v = tvl->v;
! 323: vn1[i].v = vn[i].v = 0;
! 324:
! 325: /* determine a heuristic bound of deg(GCD(p,p')) */
! 326: while ( 1 ) {
! 327: for ( i = 0; vn1[i].v; i++ )
! 328: vn1[i].n = ((unsigned int)random())%256+1;
! 329: substvp(nvl,LC(p),vn1,&p0);
! 330: if ( p0 ) {
! 331: substvp(nvl,p,vn1,&p0);
! 332: if ( sqfrchk(p0) ) {
! 333: ufctr(p0,1,&dc0);
! 334: if ( NEXT(NEXT(dc0)) == 0 ) {
! 335: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 336: *dcp = dc;
! 337: return;
! 338: } else {
! 339: for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );
! 340: break;
! 341: }
! 342: }
! 343: }
! 344: }
! 345:
! 346: /* determine the position of variables which is not allowed to
! 347: be set to 0 */
! 348: for ( i = 0; vn1[i].v; i++ ) {
! 349: x = vn1[i].n; vn1[i].n = 0;
! 350: substvp(nvl,LC(p),vn1,&p0);
! 351: if ( !p0 )
! 352: vn1[i].n = x;
! 353: else {
! 354: substvp(nvl,p,vn1,&p0);
! 355: if ( !sqfrchk(p0) )
! 356: vn1[i].n = x;
! 357: else {
! 358: ufctr(p0,1,&dc0);
! 359: for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );
! 360: if ( j > np )
! 361: vn1[i].n = x;
! 362: }
! 363: }
! 364: }
! 365: for ( i = 0; vn1[i].v; i++ )
! 366: if (vn1[i].n )
! 367: nonzero[i] = 1;
! 368:
! 369: count = 0;
! 370: while ( 1 ) {
! 371: while ( 1 ) {
! 372: count++;
! 373: for ( i = 0, j = 0; vn[i].v; i++ )
! 374: if ( vn[i].n )
! 375: vnt[j++].v = (V)((long)i);
! 376: vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
! 377: for ( i = 0; vn1[i].v; i++ )
! 378: if ( vn1[i].n ) {
! 379: if ( vn1[i].n > 0 )
! 380: vn1[i].n = sprime[vn1[i].n];
! 381: else
! 382: vn1[i].n = -sprime[-vn1[i].n];
! 383: }
! 384: if ( valideval(nvl,dcl,vn1) ) {
! 385: substvp(nvl,p,vn1,&p0);
! 386: if ( sqfrchk(p0) ) {
! 387: ufctr(p0,1,&dc0);
! 388: if ( NEXT(NEXT(dc0)) == 0 ) {
! 389: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 390: *dcp = dc;
! 391: return;
! 392: } else {
! 393: for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );
! 394: if ( i <= np )
! 395: goto MAIN;
! 396: if ( i < np )
! 397: np = i;
! 398: }
! 399: }
! 400: }
! 401: if ( nextbin(vnt,j) )
! 402: break;
! 403: }
! 404: while ( 1 ) {
! 405: next(vn);
! 406: for ( i = 0; vn[i].v; i++ )
! 407: if ( nonzero[i] && !vn[i].n )
! 408: break;
! 409: if ( !vn[i].v )
! 410: break;
! 411: }
! 412: }
! 413: MAIN :
! 414: #if 0
! 415: for ( i = 0; vn1[i].v; i++ )
! 416: fprintf(stderr,"%d ",vn1[i].n);
! 417: fprintf(stderr,"\n");
! 418: #endif
! 419: for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
! 420: fp0 = (P *)ALLOCA((np + 1)*sizeof(P));
! 421: fpt = (P *)ALLOCA((np + 1)*sizeof(P));
! 422: l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
! 423: for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) )
! 424: fp0[i-1] = COEF(dc);
! 425: #if 0
! 426: sort_by_deg(np,fp0,fpt);
! 427: sort_by_deg_rev(np-1,fpt+1,fp0+1);
! 428: #endif
! 429: #if 0
! 430: sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
! 431: sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
! 432: #endif
! 433: fp0[np] = 0;
! 434: win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
! 435: for ( k = 1, win[0] = 1, --np; ; ) {
! 436: P h0,lcg,lch;
! 437: Q c;
! 438:
! 439: g0 = fp0[win[0]];
! 440: for ( i = 1; i < k; i++ ) {
! 441: mulp(vl,g0,fp0[win[i]],&m); g0 = m;
! 442: }
! 443: mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,(Z)c,dcl,vn1,&lcg);
! 444: divsp(nvl,f0,g0,&h0);
! 445: mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,(Z)c,dcl,vn1,&lch);
! 446: mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);
! 447: if ( g ) {
! 448: *tl++ = g; divsp(vl,f,g,&t);
! 449: f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0);
! 450: for ( i = 0; i < k - 1; i++ )
! 451: for ( j = win[i] + 1; j < win[i + 1]; j++ )
! 452: fp0[j - i - 1] = fp0[j];
! 453: for ( j = win[k - 1] + 1; j <= np; j++ )
! 454: fp0[j - k] = fp0[j];
! 455: if ( ( np -= k ) < k ) break;
! 456: if ( np - win[0] + 1 < k )
! 457: if ( ++k <= np ) {
! 458: for ( i = 0; i < k; i++ ) win[i] = i + 1;
! 459: continue;
! 460: } else
! 461: break;
! 462: else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
! 463: } else {
! 464: if ( ncombi(1,np,k,win) == 0 ) {
! 465: if ( k == np ) break;
! 466: else {
! 467: for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
! 468: }
! 469: }
! 470: }
! 471: }
! 472: *tl++ = f; *tl = 0;
! 473: for ( dc0 = 0; *l; l++ ) {
! 474: NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
! 475: }
! 476: NEXT(dc) = 0; *dcp = dc0;
! 477: }
! 478:
! 479: void ufctrmain(P p,int hint,DCP *dcp)
! 480: {
! 481: ML list;
! 482: DCP dc;
! 483:
! 484: if ( NUM(p) || (UDEG(p) == 1) ) {
! 485: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 486: *dcp = dc;
! 487: } else if ( iscycm(p) )
! 488: cycm(VR(p),UDEG(p),dcp);
! 489: else if ( iscycp(p) )
! 490: cycp(VR(p),UDEG(p),dcp);
! 491: else {
! 492: if ( use_new_hensel )
! 493: hensel2(5,5,p,&list);
! 494: else
! 495: hensel(5,5,p,&list);
! 496: if ( list->n == 1 ) {
! 497: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 498: *dcp = dc;
! 499: } else
! 500: dtest(p,list,hint,dcp);
! 501: }
! 502: }
! 503:
! 504: void cycm(V v,int n,DCP *dcp)
! 505: {
! 506: register int i,j;
! 507: struct oMF *mfp;
! 508: DCP dc,dc0;
! 509:
! 510: if ( n == 1 ) {
! 511: NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
! 512: } else {
! 513: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
! 514: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
! 515: for ( i = 1, j = 0; i <= n; i++ )
! 516: if ( !(n%i) ) mfp[j++].m = i;
! 517: mkssum(v,1,1,-1,&mfp[0].f);
! 518: for ( i = 1; i < j; i++ )
! 519: calcphi(v,i,mfp);
! 520: for ( dc0 = 0, i = 0; i < j; i++ ) {
! 521: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
! 522: }
! 523: }
! 524: NEXT(dc) = 0; *dcp = dc0;
! 525: }
! 526:
! 527: void cycp(V v,int n,DCP *dcp)
! 528: {
! 529: register int i,j;
! 530: int n0;
! 531: struct oMF *mfp;
! 532: DCP dc,dc0;
! 533:
! 534: if ( n == 1 ) {
! 535: NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
! 536: } else {
! 537: n0 = n; n *= 2;
! 538: mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
! 539: bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
! 540: for ( i = 1, j = 0; i <= n; i++ )
! 541: if ( !(n%i) ) mfp[j++].m = i;
! 542: mkssum(v,1,1,-1,&mfp[0].f);
! 543: for ( i = 1; i < j; i++ )
! 544: calcphi(v,i,mfp);
! 545: for ( dc0 = 0, i = 0; i < j; i++ )
! 546: if ( n0 % mfp[i].m ) {
! 547: NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
! 548: }
! 549: }
! 550: NEXT(dc) = 0; *dcp = dc0;
! 551: }
! 552:
! 553: void calcphi(V v,int n,struct oMF *mfp)
! 554: {
! 555: register int i,m;
! 556: P t,s,tmp;
! 557:
! 558: for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
! 559: if ( !(m%mfp[i].m) ) {
! 560: mulp(CO,t,mfp[i].f,&s); t = s;
! 561: }
! 562: mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
! 563: if ( tmp )
! 564: error("calcphi: cannot happen");
! 565: }
! 566:
! 567: void mkssum(V v,int e,int s,int sgn,P *r)
! 568: {
! 569: register int i,sgnt;
! 570: DCP dc,dc0;
! 571: Z q;
! 572:
! 573: for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
! 574: if ( !dc0 ) {
! 575: NEWDC(dc0); dc = dc0;
! 576: } else {
! 577: NEWDC(NEXT(dc)); dc = NEXT(dc);
! 578: }
! 579: STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
! 580: }
! 581: NEXT(dc) = 0; MKP(v,dc0,*r);
! 582: }
! 583:
! 584: int iscycp(P f)
! 585: {
! 586: DCP dc;
! 587: dc = DC(f);
! 588:
! 589: if ( !UNIQ((Q)COEF(dc)) )
! 590: return ( 0 );
! 591: dc = NEXT(dc);
! 592: if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) )
! 593: return ( 0 );
! 594: return ( 1 );
! 595: }
! 596:
! 597: int iscycm(P f)
! 598: {
! 599: DCP dc;
! 600:
! 601: dc = DC(f);
! 602: if ( !UNIQ((Q)COEF(dc)) )
! 603: return ( 0 );
! 604: dc = NEXT(dc);
! 605: if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) )
! 606: return ( 0 );
! 607: return ( 1 );
! 608: }
! 609:
! 610: void sortfs(DCP *dcp)
! 611: {
! 612: int i,k,n,k0,d;
! 613: DCP dc,dct,t;
! 614: DCP *a;
! 615:
! 616: dc = *dcp;
! 617: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
! 618: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
! 619: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
! 620: a[i] = dct;
! 621: a[n] = 0;
! 622:
! 623: for ( i = 0; i < n; i++ ) {
! 624: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
! 625: if ( (int)UDEG(COEF(a[k])) < d ) {
! 626: k0 = k;
! 627: d = UDEG(COEF(a[k]));
! 628: }
! 629: if ( i != k0 ) {
! 630: t = a[i]; a[i] = a[k0]; a[k0] = t;
! 631: }
! 632: }
! 633: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
! 634: dct = NEXT(dct) = a[i];
! 635: NEXT(dct) = 0;
! 636: }
! 637:
! 638: void sortfsrev(DCP *dcp)
! 639: {
! 640: int i,k,n,k0,d;
! 641: DCP dc,dct,t;
! 642: DCP *a;
! 643:
! 644: dc = *dcp;
! 645: for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
! 646: a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
! 647: for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ )
! 648: a[i] = dct;
! 649: a[n] = 0;
! 650:
! 651: for ( i = 0; i < n; i++ ) {
! 652: for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ )
! 653: if ( (int)UDEG(COEF(a[k])) > d ) {
! 654: k0 = k;
! 655: d = UDEG(COEF(a[k]));
! 656: }
! 657: if ( i != k0 ) {
! 658: t = a[i]; a[i] = a[k0]; a[k0] = t;
! 659: }
! 660: }
! 661: for ( *dcp = dct = a[0], i = 1; i < n; i++ )
! 662: dct = NEXT(dct) = a[i];
! 663: NEXT(dct) = 0;
! 664: }
! 665:
! 666: void nthrootchk(P f,struct oDUM *dc,ML fp,DCP *dcp)
! 667: {
! 668: register int i,k;
! 669: int e,n,dr,tmp,t;
! 670: int *tmpp,**tmppp;
! 671: int **pp,**wpp;
! 672: LUM fpa,tpa,f0l;
! 673: UM wt,wq,ws,dvr,f0,fe;
! 674: Z lc;
! 675: int lcm;
! 676: int m,b;
! 677:
! 678: m = fp->mod; b = fp->bound; fe = *((UM *)fp->c);
! 679: e = dc->n; f0 = dc->f; nthrootz((Z)COEF(DC(f)),e,&lc);
! 680: if ( !lc ) {
! 681: *dcp = 0;
! 682: return;
! 683: }
! 684: lcm = remqi((Q)lc,m); W_LUMALLOC(DEG(f0),b,f0l);
! 685: for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
! 686: i >= 0; i-- )
! 687: *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
! 688: dtestroot(m,1,f,f0l,dc,dcp);
! 689: if ( *dcp )
! 690: return;
! 691: n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
! 692: ptolum(m,b,f,fpa);
! 693: dvr = W_UMALLOC(n); wq = W_UMALLOC(n);
! 694: wt = W_UMALLOC(n); ws = W_UMALLOC(n);
! 695: cpyum(fe,dvr); divum(m,dvr,f0,wq);
! 696: t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
! 697: for ( k = 0; k <= DEG(wq); k++ )
! 698: COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
! 699: for ( i = 1; i < b; i++ ) {
! 700: pwrlum(m,i+1,f0l,e,tpa);
! 701: for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa);
! 702: k <= n; k++ )
! 703: COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
! 704: degum(wt,n); dr = divum(m,wt,dvr,ws);
! 705: if ( dr >= 0 ) {
! 706: *dcp = 0;
! 707: return;
! 708: }
! 709: for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
! 710: pp[k][i] = COEF(ws)[k];
! 711: dtestroot(m,i+1,f,f0l,dc,dcp);
! 712: if ( *dcp )
! 713: return;
! 714: }
! 715: }
! 716:
! 717: void sqfrp(VL vl,P f,DCP *dcp)
! 718: {
! 719: P c,p;
! 720: DCP dc,dc0;
! 721:
! 722: if ( !f || NUM(f) ) {
! 723: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
! 724: COEF(dc0) = f; NEXT(dc0) = 0;
! 725: } else if ( !qpcheck((Obj)f) )
! 726: error("sqfrp : invalid argument");
! 727: else {
! 728: NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
! 729: ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
! 730: COEF(dc0) = c; NEXT(dc0) = dc;
! 731: adjsgn(f,dc0);
! 732: }
! 733: }
! 734:
! 735: /*
! 736: * f : must be a poly with int coef, ignore int content
! 737: */
! 738: void msqfr(VL vl,P f,DCP *dcp)
! 739: {
! 740: DCP dc,dct,dcs;
! 741: P c,p,t,s,pc;
! 742: VL mvl;
! 743:
! 744: ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
! 745: if ( NUM(p) ) {
! 746: *dcp = dc;
! 747: return;
! 748: }
! 749: mindegp(vl,p,&mvl,&s);
! 750: #if 0
! 751: minlcdegp(vl,p,&mvl,&s);
! 752: min_common_vars_in_coefp(vl,p,&mvl,&s);
! 753: #endif
! 754: pcp(mvl,s,&p,&pc);
! 755: if ( !NUM(pc) ) {
! 756: msqfr(mvl,pc,&dcs);
! 757: if ( !dc )
! 758: dc = dcs;
! 759: else {
! 760: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
! 761: NEXT(dct) = dcs;
! 762: }
! 763: }
! 764: msqfrmain(mvl,p,&dcs);
! 765: for ( dct = dcs; dct; dct = NEXT(dct) ) {
! 766: reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
! 767: }
! 768: if ( !dc )
! 769: *dcp = dcs;
! 770: else {
! 771: for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
! 772: NEXT(dct) = dcs; *dcp = dc;
! 773: }
! 774: }
! 775:
! 776: void usqp(P f,DCP *dcp)
! 777: {
! 778: int index,nindex;
! 779: P g,c,h;
! 780: DCP dc;
! 781:
! 782: ptozp(f,1,(Q *)&c,&h);
! 783: if ( sgnq((Q)LC(h)) < 0 )
! 784: chsgnp(h,&g);
! 785: else
! 786: g = h;
! 787: for ( index = 0, dc = 0; !dc; index = nindex )
! 788: hsq(index,5,g,&nindex,&dc);
! 789: *dcp = dc;
! 790: }
! 791:
! 792: void msqfrmain(VL vl,P p,DCP *dcp)
! 793: {
! 794: int i,j;
! 795: VL nvl,tvl;
! 796: VN vn,vnt,vn1;
! 797: P gg,tmp,p0,g;
! 798: DCP dc,dc0,dcr,dcr0;
! 799: int nv,d,d1;
! 800: int found;
! 801: VN svn1;
! 802: P sp0;
! 803: DCP sdc0;
! 804:
! 805: if ( deg(VR(p),p) == 1 ) {
! 806: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 807: *dcp = dc;
! 808: return;
! 809: }
! 810: clctv(vl,p,&nvl);
! 811: for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
! 812: if ( nv == 1 ) {
! 813: usqp(p,dcp);
! 814: return;
! 815: }
! 816: #if 0
! 817: if ( heusqfrchk(nvl,p) ) {
! 818: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 819: *dcp = dc;
! 820: return;
! 821: }
! 822: #endif
! 823: W_CALLOC(nv,struct oVN,vn);
! 824: W_CALLOC(nv,struct oVN,vnt);
! 825: W_CALLOC(nv,struct oVN,vn1);
! 826: W_CALLOC(nv,struct oVN,svn1);
! 827: for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ )
! 828: vn1[i].v = vn[i].v = tvl->v;
! 829: vn1[i].v = vn[i].v = 0;
! 830:
! 831: /* determine a heuristic bound of deg(GCD(p,p')) */
! 832: while ( 1 ) {
! 833: for ( i = 0; vn1[i].v; i++ )
! 834: vn1[i].n = ((unsigned int)random())%256+1;
! 835: substvp(nvl,LC(p),vn1,&tmp);
! 836: if ( tmp ) {
! 837: substvp(nvl,p,vn1,&p0);
! 838: usqp(p0,&dc0);
! 839: for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
! 840: if ( DEG(dc) )
! 841: d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
! 842: if ( d1 == 0 ) {
! 843: /* p is squarefree */
! 844: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
! 845: *dcp = dc;
! 846: return;
! 847: } else {
! 848: d = d1+1; /* XXX : try searching better evaluation */
! 849: found = 0;
! 850: break;
! 851: }
! 852: }
! 853: }
! 854:
! 855: for ( dcr0 = 0, g = p; ; ) {
! 856: while ( 1 ) {
! 857: for ( i = 0, j = 0; vn[i].v; i++ )
! 858: if ( vn[i].n ) vnt[j++].v = (V)((long)i);
! 859: vnt[j].n = 0;
! 860:
! 861: mulsgn(vn,vnt,j,vn1);
! 862: substvp(nvl,LC(g),vn1,&tmp);
! 863: if ( tmp ) {
! 864: substvp(nvl,g,vn1,&p0);
! 865: usqp(p0,&dc0);
! 866: for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
! 867: if ( DEG(dc) )
! 868: d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
! 869:
! 870: if ( d1 == 0 ) {
! 871: NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
! 872: if ( !dcr0 )
! 873: dcr0 = dc;
! 874: else
! 875: NEXT(dcr) = dc;
! 876: *dcp = dcr0;
! 877: return;
! 878: }
! 879:
! 880: if ( d < d1 )
! 881: goto END;
! 882: if ( d > d1 ) {
! 883: d = d1;
! 884: found = 1;
! 885: bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
! 886: sp0 = p0; sdc0 = dc0;
! 887: goto END;
! 888: }
! 889: /* d1 == d */
! 890: if ( found ) {
! 891: found = 0;
! 892: #if 0
! 893: if ( d > d1 ) {
! 894: d = d1;
! 895: /*} */
! 896: #endif
! 897: msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
! 898: if ( dc ) {
! 899: if ( !dcr0 )
! 900: dcr0 = dc;
! 901: else
! 902: NEXT(dcr) = dc;
! 903: for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
! 904: if ( NUM(g) ) {
! 905: NEXT(dcr) = 0; *dcp = dcr0;
! 906: return;
! 907: }
! 908: d = deg(VR(g),g);
! 909: }
! 910: }
! 911: }
! 912: END:
! 913: if ( nextbin(vnt,j) )
! 914: break;
! 915: }
! 916: next(vn);
! 917: }
! 918: }
! 919:
! 920: void msqfrmainmain(VL vl,P p,VN vn,P p0,DCP dc0,DCP *dcp,P *pp)
! 921: {
! 922: int i,j,k,np;
! 923: DCP *a;
! 924: DCP dc,dcr,dcr0,dct;
! 925: P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
! 926: Z q;
! 927: V v;
! 928:
! 929: for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
! 930: a = (DCP *)ALLOCA((np + 1)*sizeof(DCP));
! 931: for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) )
! 932: a[np-i-1] = dc;
! 933:
! 934: for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p);
! 935: i < np; i++ ) {
! 936: if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
! 937: NEXTDC(dcr0,dcr);
! 938: DEG(dcr) = DEG(a[i]);
! 939: COEF(dcr) = f;
! 940: f = (P)ONE;
! 941: } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
! 942: diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
! 943: if ( divtpz(vl,f,t,&s) ) {
! 944: NEXTDC(dcr0,dcr);
! 945: DEG(dcr) = DEG(a[i]);
! 946: COEF(dcr) = s;
! 947: f = (P)ONE;
! 948: } else
! 949: break;
! 950: } else {
! 951: for ( t = f, t0 = f0,
! 952: j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
! 953: diffp(vl,t,v,&s); t = s;
! 954: diffp(vl,t0,v,&s); t0 = s;
! 955: }
! 956: factorialz(k,&q);
! 957: divsp(vl,t,(P)q,&s);
! 958: for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
! 959: if ( DEG(dct) ) {
! 960: MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx);
! 961: divsp(vl,s,xx,&d);
! 962: } else {
! 963: xx = (P)ONE; d = s;
! 964: }
! 965: divsp(vl,t0,xx,&t);
! 966: divsp(vl,t,(P)q,&s);
! 967: ptozp(s,1,(Q *)&t,&d0);
! 968:
! 969: for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
! 970: if ( DEG(dct) )
! 971: divsp(vl,COEF(a[i]),xx,&g0);
! 972: else {
! 973: xx = (P)ONE; g0 = COEF(a[i]);
! 974: }
! 975:
! 976: pcp(vl,d,&t,&u); d = t;
! 977: ptozp(g0,1,(Q *)&u,&t); g0 = t;
! 978:
! 979: {
! 980: DCP dca,dcb;
! 981:
! 982: fctrp(vl,LC(d),&dca);
! 983: for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
! 984: if ( NUM(COEF(dcb)) ) {
! 985: mulp(vl,u,COEF(dcb),&t); u = t;
! 986: } else {
! 987: Z qq;
! 988: P tt;
! 989:
! 990: pwrp(vl,COEF(dcb),DEG(a[i]),&s);
! 991: for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
! 992: STOQ(j,qq);
! 993: if ( cmpz(qq,DEG(dcb)) > 0 )
! 994: qq = DEG(dcb);
! 995: pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
! 996: }
! 997: }
! 998: divsp(vl,d0,g0,&h0);
! 999: }
! 1000: mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s);
! 1001: if ( s ) {
! 1002: mulp(vl,s,xx,&g);
! 1003: pwrp(vl,g,DEG(a[i]),&t);
! 1004: if ( divtpz(vl,f,t,&s) ) {
! 1005: NEXTDC(dcr0,dcr);
! 1006: DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
! 1007: f = s; substvp(vl,f,vn,&f0);
! 1008: } else
! 1009: break;
! 1010: } else
! 1011: break;
! 1012: }
! 1013: }
! 1014: *pp = f;
! 1015: if ( dcr0 )
! 1016: NEXT(dcr) = 0;
! 1017: *dcp = dcr0;
! 1018: }
! 1019:
! 1020: void mfctrhen2(VL vl,VN vn,P f,P f0,P g0,P h0,P lcg,P lch,P *gp)
! 1021: {
! 1022: V v;
! 1023: P f1,lc,lc0,lcg0,lch0;
! 1024: P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
! 1025: Q bb,s;
! 1026: Z qk;
! 1027: Q cbd;
! 1028: int dbd;
! 1029: int d;
! 1030:
! 1031: if ( NUM(g0) ) {
! 1032: *gp = (P)ONE;
! 1033: return;
! 1034: }
! 1035:
! 1036: v = VR(f); d = deg(v,f);
! 1037: if ( d == deg(v,g0) ) {
! 1038: pcp(vl,f,gp,&tmp);
! 1039: return;
! 1040: }
! 1041:
! 1042: mulp(vl,lcg,lch,&lc);
! 1043: if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
! 1044: *gp = 0;
! 1045: return;
! 1046: }
! 1047: mulp(vl,(P)s,f,&f1);
! 1048: dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
! 1049:
! 1050: substvp(vl,lc,vn,&lc0);
! 1051: divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0);
! 1052: substvp(vl,lcg,vn,&lcg0);
! 1053: divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0);
! 1054: substvp(vl,lch,vn,&lch0);
! 1055: divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0);
! 1056: addq(cbd,cbd,&bb);
! 1057: henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
! 1058: henmv(vl,vn,f1,gk,hk,ak,bk,
! 1059: lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh);
! 1060:
! 1061: if ( divtpz(vl,f1,ggg,&gggr) )
! 1062: pcp(vl,ggg,gp,&tmp);
! 1063: else
! 1064: *gp = 0;
! 1065: }
! 1066:
! 1067: int sqfrchk(P p)
! 1068: {
! 1069: Q c;
! 1070: P f;
! 1071: DCP dc;
! 1072:
! 1073: ptozp(p,sgnq((Q)UCOEF(p)),&c,&f); usqp(f,&dc);
! 1074: if ( NEXT(dc) || !UNIQ(DEG(dc)) )
! 1075: return ( 0 );
! 1076: else
! 1077: return ( 1 );
! 1078: }
! 1079:
! 1080: int cycchk(P p)
! 1081: {
! 1082: Q c;
! 1083: P f;
! 1084:
! 1085: ptozp(p,sgnq((Q)UCOEF(p)),&c,&f);
! 1086: if ( iscycp(f) || iscycm(f) )
! 1087: return 0;
! 1088: else
! 1089: return 1;
! 1090: }
! 1091:
! 1092: int zerovpchk(VL vl,P p,VN vn)
! 1093: {
! 1094: P t;
! 1095:
! 1096: substvp(vl,p,vn,&t);
! 1097: if ( t )
! 1098: return ( 0 );
! 1099: else
! 1100: return ( 1 );
! 1101: }
! 1102:
! 1103: int valideval(VL vl,DCP dc,VN vn)
! 1104: {
! 1105: DCP dct;
! 1106: Z *a;
! 1107: int i,j,n;
! 1108: Z q,r;
! 1109:
! 1110: for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
! 1111: W_CALLOC(n,Z,a);
! 1112: for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
! 1113: substvp(vl,COEF(dct),vn,(P *)&a[i]);
! 1114: if ( !a[i] )
! 1115: return ( 0 );
! 1116:
! 1117: for ( j = 0; j < i; j++ ) {
! 1118: divqrz(a[j],a[i],&q,&r);
! 1119: if ( !r )
! 1120: return ( 0 );
! 1121: divqrz(a[i],a[j],&q,&r);
! 1122: if ( !r )
! 1123: return ( 0 );
! 1124: }
! 1125: }
! 1126: return ( 1 );
! 1127: }
! 1128:
! 1129: void estimatelc(VL vl,Z c,DCP dc,VN vn,P *lcp)
! 1130: {
! 1131: int i;
! 1132: DCP dct;
! 1133: P r,s,t;
! 1134: Z c0,c1,c2;
! 1135:
! 1136: for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
! 1137: if ( NUM(COEF(dct)) ) {
! 1138: mulp(vl,r,COEF(dct),&s); r = s;
! 1139: } else {
! 1140: substvp(vl,COEF(dct),vn,(P *)&c0);
! 1141: for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
! 1142: divz(c1,c0,&c2);
! 1143: if ( !INT(c2) )
! 1144: break;
! 1145: else
! 1146: c1 = c2;
! 1147: }
! 1148: if ( i ) {
! 1149: STOQ(i,c1);
! 1150: pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
! 1151: }
! 1152: }
! 1153: }
! 1154: *lcp = r;
! 1155: }
! 1156:
! 1157: void monomialfctr(VL vl,P p,P *pr,DCP *dcp)
! 1158: {
! 1159: VL nvl,avl;
! 1160: Z d;
! 1161: P f,t,s;
! 1162: DCP dc0,dc;
! 1163:
! 1164: clctv(vl,p,&nvl);
! 1165: for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
! 1166: getmindeg(avl->v,f,&d);
! 1167: if ( d ) {
! 1168: MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
! 1169: pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
! 1170: }
! 1171: }
! 1172: if ( dc0 )
! 1173: NEXT(dc) = 0;
! 1174: *pr = f; *dcp = dc0;
! 1175: }
! 1176:
! 1177: void afctr(VL vl,P p0,P p,DCP *dcp)
! 1178: {
! 1179: DCP dc,dc0,dcr,dct,dcs;
! 1180: P t;
! 1181: VL nvl;
! 1182:
! 1183: if ( VR(p) == VR(p0) ) {
! 1184: NEWDC(dc);
! 1185: DEG(dc) = ONE;
! 1186: COEF(dc) = p;
! 1187: NEXT(dc) = 0;
! 1188: *dcp = dc;
! 1189: return;
! 1190: }
! 1191:
! 1192: clctv(vl,p,&nvl);
! 1193: if ( !NEXT(nvl) )
! 1194: ufctr(p,1,&dc);
! 1195: else {
! 1196: sqa(vl,p0,p,&dc);
! 1197: for ( dct = dc; dct; dct = NEXT(dct) ) {
! 1198: pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
! 1199: }
! 1200: }
! 1201: if ( NUM(COEF(dc)) )
! 1202: dcr = NEXT(dc);
! 1203: else
! 1204: dcr = dc;
! 1205: for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
! 1206: afctrmain(vl,p0,COEF(dcr),1,&dcs);
! 1207:
! 1208: for ( dct = dcs; dct; dct = NEXT(dct) )
! 1209: DEG(dct) = DEG(dcr);
! 1210: if ( !dc0 )
! 1211: dc0 = dcs;
! 1212: else {
! 1213: for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
! 1214: NEXT(dct) = dcs;
! 1215: }
! 1216: }
! 1217: *dcp = dc0;
! 1218: }
! 1219:
! 1220: void afctrmain(VL vl,P p0,P p,int init,DCP *dcp)
! 1221: {
! 1222: P x,y,s,m,a,t,u,pt,pt1,res,g;
! 1223: Z q;
! 1224: DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
! 1225: V v,v0;
! 1226:
! 1227: if ( !cmpz(DEG(DC(p)),ONE) ) {
! 1228: NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
! 1229: pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
! 1230: return;
! 1231: }
! 1232:
! 1233: v = VR(p); MKV(v,x);
! 1234: v0 = VR(p0); MKV(v0,y);
! 1235: STOQ(init,q),s = (P)q;
! 1236: mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a);
! 1237: substp(vl,p,v,t,&pt);
! 1238: remsdcp(vl,pt,p0,&pt1);
! 1239:
! 1240: /*
! 1241: if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) )
! 1242: resultp(vl,v0,p0,pt1,&res);
! 1243: else
! 1244: srcrnorm(vl,v0,pt1,p0,&res);
! 1245: */
! 1246: #if 0
! 1247: srcr(vl,v0,pt1,p0,&res);
! 1248: #endif
! 1249: resultp(vl,v0,p0,pt1,&res);
! 1250: usqp(res,&dcsq0);
! 1251: for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
! 1252: if ( UNIQ(DEG(dct)) )
! 1253: ufctr(COEF(dct),deg(v0,p0),&dcs);
! 1254: else
! 1255: ufctr(COEF(dct),1,&dcs);
! 1256: for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
! 1257: DEG(dcr) = DEG(dct);
! 1258: if ( !dc0 ) {
! 1259: dc0 = NEXT(dcs);
! 1260: dc = dc0;
! 1261: } else {
! 1262: for ( ; NEXT(dc); dc = NEXT(dc) );
! 1263: NEXT(dc) = NEXT(dcs);
! 1264: }
! 1265: }
! 1266: sortfs(&dc0);
! 1267:
! 1268: for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
! 1269: if ( !UNIQ(DEG(dc)) ) {
! 1270: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
! 1271: gcda(vl,p0,s,g,&u);
! 1272: if ( !NUM(u) && (VR(u) == v)) {
! 1273: afctrmain(vl,p0,u,init+1,&dct);
! 1274: for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
! 1275: mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
! 1276: }
! 1277: pdiva(vl,p0,g,t,&s); g = s;
! 1278: if ( !dcr0 )
! 1279: dcr0 = dct;
! 1280: else
! 1281: NEXT(dcr) = dct;
! 1282: for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
! 1283: }
! 1284: } else {
! 1285: substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s);
! 1286: gcda(vl,p0,s,g,&u);
! 1287: if ( !NUM(u) && (VR(u) == v)) {
! 1288: NEXTDC(dcr0,dcr);
! 1289: DEG(dcr) = ONE;
! 1290: COEF(dcr) = u;
! 1291: pdiva(vl,p0,g,u,&t); g = t;
! 1292: }
! 1293: }
! 1294: }
! 1295: if ( dcr0 )
! 1296: NEXT(dcr) = 0;
! 1297: *dcp = dcr0;
! 1298: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to F.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine