Annotation of OpenXM_contrib2/asir2018/engine/PU.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:
! 52: void reorderp(VL nvl,VL ovl,P p,P *pr)
! 53: {
! 54: DCP dc;
! 55: P x,m,s,t,c;
! 56:
! 57: if ( !p )
! 58: *pr = 0;
! 59: else if ( NUM(p) )
! 60: *pr = p;
! 61: else {
! 62: MKV(VR(p),x);
! 63: for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 64: reorderp(nvl,ovl,COEF(dc),&c);
! 65: if ( DEG(dc) ) {
! 66: pwrp(nvl,x,DEG(dc),&t); mulp(nvl,c,t,&m);
! 67: addp(nvl,m,s,&t);
! 68: } else
! 69: addp(nvl,s,c,&t);
! 70: s = t;
! 71: }
! 72: *pr = s;
! 73: }
! 74: }
! 75:
! 76: void substp(VL vl,P p,V v0,P p0,P *pr)
! 77: {
! 78: P x,t,m,c,s,a;
! 79: DCP dc;
! 80: Z d;
! 81:
! 82: if ( !p )
! 83: *pr = 0;
! 84: else if ( NUM(p) )
! 85: *pr = p;
! 86: else if ( VR(p) != v0 ) {
! 87: MKV(VR(p),x);
! 88: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 89: substp(vl,COEF(dc),v0,p0,&t);
! 90: if ( DEG(dc) ) {
! 91: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 92: addp(vl,m,c,&a);
! 93: c = a;
! 94: } else {
! 95: addp(vl,t,c,&a);
! 96: c = a;
! 97: }
! 98: }
! 99: *pr = c;
! 100: } else {
! 101: dc = DC(p);
! 102: c = COEF(dc);
! 103: for ( d = DEG(dc), dc = NEXT(dc);
! 104: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 105: subz(d,DEG(dc),(Z *)&t); pwrp(vl,p0,(Z)t,&s);
! 106: mulp(vl,s,c,&m);
! 107: addp(vl,m,COEF(dc),&c);
! 108: }
! 109: if ( d ) {
! 110: pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
! 111: c = m;
! 112: }
! 113: *pr = c;
! 114: }
! 115: }
! 116:
! 117: void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr);
! 118:
! 119: void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr)
! 120: {
! 121: P x,t,m,c,s,a,p0,c1;
! 122: DCP dc;
! 123: Z d;
! 124: V v;
! 125: int i;
! 126:
! 127: if ( !p )
! 128: *pr = 0;
! 129: else if ( NUM(p) )
! 130: *pr = p;
! 131: else {
! 132: v = VR(p);
! 133: for ( i = 0; i < nv; i++ ) if ( vvect[i] == v ) break;
! 134: if ( svect[i] && OID(svect[i]) < 0 ) {
! 135: MKV(VR(p),x);
! 136: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 137: substpp(vl,COEF(dc),vvect,svect,nv,&t);
! 138: if ( DEG(dc) ) {
! 139: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 140: addp(vl,m,c,&a);
! 141: c = a;
! 142: } else {
! 143: addp(vl,t,c,&a);
! 144: c = a;
! 145: }
! 146: }
! 147: *pr = c;
! 148: } else {
! 149: p0 = svect[i];
! 150: dc = DC(p);
! 151: substpp(vl,COEF(dc),vvect,svect,nv,&c);
! 152: for ( d = DEG(dc), dc = NEXT(dc);
! 153: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 154: subz(d,DEG(dc),(Z *)&t); pwrp(vl,p0,(Z)t,&s);
! 155: mulp(vl,s,c,&m);
! 156: substpp(vl,COEF(dc),vvect,svect,nv,&c1);
! 157: addp(vl,m,c1,&c);
! 158: }
! 159: if ( d ) {
! 160: pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
! 161: c = m;
! 162: }
! 163: *pr = c;
! 164: }
! 165: }
! 166: }
! 167:
! 168: void detp(VL vl,P **rmat,int n,P *dp)
! 169: {
! 170: int i,j,k,l,sgn,nmin,kmin,lmin,ntmp;
! 171: P mjj,mij,t,s,u,d;
! 172: P **mat;
! 173: P *mi,*mj;
! 174:
! 175: mat = (P **)almat_pointer(n,n);
! 176: for ( i = 0; i < n; i++ )
! 177: for ( j = 0; j < n; j++ )
! 178: mat[i][j] = rmat[i][j];
! 179: for ( j = 0, d = (P)ONE, sgn = 1; j < n; j++ ) {
! 180: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 181: if ( i == n ) {
! 182: *dp = 0; return;
! 183: }
! 184: nmin = nmonop(mat[i][j]);
! 185: kmin=i; lmin=j;
! 186: for ( k = j; k < n; k++ )
! 187: for ( l = j; l < n; l++ )
! 188: if ( mat[k][l] && ((ntmp=nmonop(mat[k][l])) < nmin) ) {
! 189: kmin = k; lmin = l; nmin = ntmp;
! 190: }
! 191: if ( kmin != j ) {
! 192: mj = mat[j]; mat[j] = mat[kmin]; mat[kmin] = mj; sgn = -sgn;
! 193: }
! 194: if ( lmin != j ) {
! 195: for ( k = j; k < n; k++ ) {
! 196: t = mat[k][j]; mat[k][j] = mat[k][lmin]; mat[k][lmin] = t;
! 197: }
! 198: sgn = -sgn;
! 199: }
! 200: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 201: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 202: mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
! 203: subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
! 204: }
! 205: d = mjj;
! 206: }
! 207: if ( sgn < 0 )
! 208: chsgnp(d,dp);
! 209: else
! 210: *dp = d;
! 211: }
! 212:
! 213: void invmatp(VL vl,P **rmat,int n,P ***imatp,P *dnp)
! 214: {
! 215: int i,j,k,l,n2;
! 216: P mjj,mij,t,s,u,d;
! 217: P **mat,**imat;
! 218: P *mi,*mj,*w;
! 219:
! 220: n2 = n<<1;
! 221: mat = (P **)almat_pointer(n,n2);
! 222: for ( i = 0; i < n; i++ ) {
! 223: for ( j = 0; j < n; j++ )
! 224: mat[i][j] = rmat[i][j];
! 225: mat[i][i+n] = (P)ONE;
! 226: }
! 227: for ( j = 0, d = (P)ONE; j < n; j++ ) {
! 228: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 229: if ( i == n ) {
! 230: error("invmatp : input is singular");
! 231: }
! 232: for ( k = i; k < n; k++ )
! 233: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 234: i = k;
! 235: if ( j != i ) {
! 236: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj;
! 237: }
! 238: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 239: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n2; k++ ) {
! 240: mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
! 241: subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
! 242: }
! 243: d = mjj;
! 244: }
! 245: /* backward substitution */
! 246: w = (P *)ALLOCA(n2*sizeof(P));
! 247: for ( i = n-2; i >= 0; i-- ) {
! 248: bzero(w,n2*sizeof(P));
! 249: for ( k = i+1; k < n; k++ )
! 250: for ( l = k, u = mat[i][l]; l < n2; l++ ) {
! 251: mulp(vl,mat[k][l],u,&t); addp(vl,w[l],t,&s); w[l] = s;
! 252: }
! 253: for ( j = i, u = mat[i][j]; j < n2; j++ ) {
! 254: mulp(vl,mat[i][j],d,&t); subp(vl,t,w[j],&s);
! 255: divsp(vl,s,u,&mat[i][j]);
! 256: }
! 257: }
! 258: imat = (P **)almat_pointer(n,n);
! 259: for ( i = 0; i < n; i++ )
! 260: for ( j = 0; j < n; j++ )
! 261: imat[i][j] = mat[i][j+n];
! 262: *imatp = imat;
! 263: *dnp = d;
! 264: }
! 265:
! 266: void reordvar(VL vl,V v,VL *nvlp)
! 267: {
! 268: VL nvl,nvl0;
! 269:
! 270: for ( NEWVL(nvl0), nvl0->v = v, nvl = nvl0;
! 271: vl; vl = NEXT(vl) )
! 272: if ( vl->v == v )
! 273: continue;
! 274: else {
! 275: NEWVL(NEXT(nvl));
! 276: nvl = NEXT(nvl);
! 277: nvl->v = vl->v;
! 278: }
! 279: NEXT(nvl) = 0;
! 280: *nvlp = nvl0;
! 281: }
! 282:
! 283: void gcdprsp(VL vl,P p1,P p2,P *pr)
! 284: {
! 285: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
! 286: V v1,v2;
! 287:
! 288: if ( !p1 )
! 289: ptozp0(p2,pr);
! 290: else if ( !p2 )
! 291: ptozp0(p1,pr);
! 292: else if ( NUM(p1) || NUM(p2) )
! 293: *pr = (P)ONE;
! 294: else {
! 295: ptozp0(p1,&g1); ptozp0(p2,&g2);
! 296: if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
! 297: gcdcp(vl,g1,&gc1); divsp(vl,g1,gc1,&gp1);
! 298: gcdcp(vl,g2,&gc2); divsp(vl,g2,gc2,&gp2);
! 299: gcdprsp(vl,gc1,gc2,&gcr);
! 300: sprs(vl,v1,gp1,gp2,&g);
! 301:
! 302: if ( VR(g) == v1 ) {
! 303: ptozp0(g,&gp);
! 304: gcdcp(vl,gp,&gc); divsp(vl,gp,gc,&gp1);
! 305: mulp(vl,gp1,gcr,pr);
! 306:
! 307: } else
! 308: *pr = gcr;
! 309: } else {
! 310: while ( v1 != vl->v && v2 != vl->v )
! 311: vl = NEXT(vl);
! 312: if ( v1 == vl->v ) {
! 313: gcdcp(vl,g1,&gc1); gcdprsp(vl,gc1,g2,pr);
! 314: } else {
! 315: gcdcp(vl,g2,&gc2); gcdprsp(vl,gc2,g1,pr);
! 316: }
! 317: }
! 318: }
! 319: }
! 320:
! 321: void gcdcp(VL vl,P p,P *pr)
! 322: {
! 323: P g,g1;
! 324: DCP dc;
! 325:
! 326: if ( NUM(p) )
! 327: *pr = (P)ONE;
! 328: else {
! 329: dc = DC(p);
! 330: g = COEF(dc);
! 331: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 332: gcdprsp(vl,g,COEF(dc),&g1);
! 333: g = g1;
! 334: }
! 335: *pr = g;
! 336: }
! 337: }
! 338:
! 339: void sprs(VL vl,V v,P p1,P p2,P *pr)
! 340: {
! 341: P q1,q2,m,m1,m2,x,h,r,g1,g2;
! 342: int d;
! 343: Z dq;
! 344: VL nvl;
! 345:
! 346: reordvar(vl,v,&nvl);
! 347: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 348:
! 349: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 350: *pr = 0;
! 351: return;
! 352: }
! 353:
! 354: if ( deg(v,q1) >= deg(v,q2) ) {
! 355: g1 = q1; g2 = q2;
! 356: } else {
! 357: g2 = q1; g1 = q2;
! 358: }
! 359:
! 360: for ( h = (P)ONE, x = (P)ONE; ; ) {
! 361: if ( !deg(v,g2) )
! 362: break;
! 363:
! 364: premp(nvl,g1,g2,&r);
! 365: if ( !r )
! 366: break;
! 367:
! 368: d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
! 369: pwrp(nvl,h,dq,&m); mulp(nvl,m,x,&m1); g1 = g2;
! 370: divsp(nvl,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
! 371: pwrp(nvl,x,dq,&m1); mulp(nvl,m1,h,&m2);
! 372: divsp(nvl,m2,m,&h);
! 373: }
! 374: *pr = g2;
! 375: }
! 376:
! 377: void resultp(VL vl,V v,P p1,P p2,P *pr)
! 378: {
! 379: P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
! 380: int d,d1,d2,j,k;
! 381: VL nvl;
! 382: Z dq;
! 383:
! 384: if ( !p1 || !p2 ) {
! 385: *pr = 0; return;
! 386: }
! 387: reordvar(vl,v,&nvl);
! 388: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 389:
! 390: if ( VR(q1) != v )
! 391: if ( VR(q2) != v ) {
! 392: *pr = 0;
! 393: return;
! 394: } else {
! 395: d = deg(v,q2); STOQ(d,dq);
! 396: pwrp(vl,q1,dq,pr);
! 397: return;
! 398: }
! 399: else if ( VR(q2) != v ) {
! 400: d = deg(v,q1); STOQ(d,dq);
! 401: pwrp(vl,q2,dq,pr);
! 402: return;
! 403: }
! 404:
! 405: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 406: *pr = 0;
! 407: return;
! 408: }
! 409:
! 410: d1 = deg(v,q1); d2 = deg(v,q2);
! 411: if ( d1 > d2 ) {
! 412: g1 = q1; g2 = q2; adj = (P)ONE;
! 413: } else if ( d1 < d2 ) {
! 414: g2 = q1; g1 = q2;
! 415: if ( (d1 % 2) && (d2 % 2) ) {
! 416: STOQ(-1,dq); adj = (P)dq;
! 417: } else
! 418: adj = (P)ONE;
! 419: } else {
! 420: premp(nvl,q1,q2,&t);
! 421: d = deg(v,t); STOQ(d,dq); pwrp(nvl,LC(q2),dq,&adj);
! 422: g1 = q2; g2 = t;
! 423: if ( d1 % 2 ) {
! 424: chsgnp(adj,&t); adj = t;
! 425: }
! 426: }
! 427: d1 = deg(v,g1); j = d1 - 1;
! 428:
! 429: for ( lc = (P)ONE; ; ) {
! 430: if ( ( k = deg(v,g2) ) < 0 ) {
! 431: *pr = 0;
! 432: return;
! 433: }
! 434:
! 435: if ( k == j )
! 436: if ( !k ) {
! 437: divsp(nvl,g2,adj,pr);
! 438: return;
! 439: } else {
! 440: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
! 441: divsp(nvl,r,m,&q);
! 442: g1 = g2; g2 = q;
! 443: lc = LC(g1); /* g1 is not const */
! 444: j = k - 1;
! 445: }
! 446: else {
! 447: d = j - k; STOQ(d,dq);
! 448: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
! 449: mulp(nvl,g2,m,&m1);
! 450: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
! 451: if ( k == 0 ) {
! 452: divsp(nvl,t,adj,pr);
! 453: return;
! 454: } else {
! 455: premp(nvl,g1,g2,&r);
! 456: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
! 457: divsp(nvl,r,m2,&q);
! 458: g1 = t; g2 = q;
! 459: if ( d % 2 ) {
! 460: chsgnp(g2,&t); g2 = t;
! 461: }
! 462: lc = LC(g1); /* g1 is not const */
! 463: j = k - 1;
! 464: }
! 465: }
! 466: }
! 467: }
! 468:
! 469: void srch2(VL vl,V v,P p1,P p2,P *pr)
! 470: {
! 471: P q1,q2,m,m1,m2,lc,q,r,t,s,g1,g2,adj;
! 472: int d,d1,d2,j,k;
! 473: VL nvl;
! 474: Z dq;
! 475:
! 476: reordvar(vl,v,&nvl);
! 477: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 478:
! 479: if ( VR(q1) != v )
! 480: if ( VR(q2) != v ) {
! 481: *pr = 0;
! 482: return;
! 483: } else {
! 484: d = deg(v,q2); STOQ(d,dq);
! 485: pwrp(vl,q1,dq,pr);
! 486: return;
! 487: }
! 488: else if ( VR(q2) != v ) {
! 489: d = deg(v,q1); STOQ(d,dq);
! 490: pwrp(vl,q2,dq,pr);
! 491: return;
! 492: }
! 493:
! 494: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 495: *pr = 0;
! 496: return;
! 497: }
! 498:
! 499: if ( deg(v,q1) >= deg(v,q2) ) {
! 500: g1 = q1; g2 = q2;
! 501: } else {
! 502: g2 = q1; g1 = q2;
! 503: }
! 504:
! 505:
! 506: if ( ( d1 = deg(v,g1) ) > ( d2 = deg(v,g2) ) ) {
! 507: j = d1 - 1;
! 508: adj = (P)ONE;
! 509: } else {
! 510: premp(nvl,g1,g2,&t);
! 511: d = deg(v,t); STOQ(d,dq);
! 512: pwrp(nvl,LC(g2),dq,&adj);
! 513: g1 = g2; g2 = t;
! 514: j = deg(v,g1) - 1;
! 515: }
! 516:
! 517: for ( lc = (P)ONE; ; ) {
! 518: if ( ( k = deg(v,g2) ) < 0 ) {
! 519: *pr = 0;
! 520: return;
! 521: }
! 522:
! 523: ptozp(g1,1,(Q *)&t,&s); g1 = s; ptozp(g2,1,(Q *)&t,&s); g2 = s;
! 524: if ( k == j )
! 525: if ( !k ) {
! 526: divsp(nvl,g2,adj,pr);
! 527: return;
! 528: } else {
! 529: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
! 530: divsp(nvl,r,m,&q);
! 531: g1 = g2; g2 = q;
! 532: lc = LC(g1); /* g1 is not const */
! 533: j = k - 1;
! 534: }
! 535: else {
! 536: d = j - k; STOQ(d,dq);
! 537: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
! 538: mulp(nvl,g2,m,&m1);
! 539: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
! 540: if ( k == 0 ) {
! 541: divsp(nvl,t,adj,pr);
! 542: return;
! 543: } else {
! 544: premp(nvl,g1,g2,&r);
! 545: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
! 546: divsp(nvl,r,m2,&q);
! 547: g1 = t; g2 = q;
! 548: if ( d % 2 ) {
! 549: chsgnp(g2,&t); g2 = t;
! 550: }
! 551: lc = LC(g1); /* g1 is not const */
! 552: j = k - 1;
! 553: }
! 554: }
! 555: }
! 556: }
! 557:
! 558: void srcr(VL vl,V v,P p1,P p2,P *pr)
! 559: {
! 560: P q1,q2,c,c1;
! 561: P tg,tg1,tg2,resg;
! 562: int index,mod;
! 563: Q a,b,f,q,s,u,w;
! 564: Z n,m,t;
! 565: VL nvl;
! 566:
! 567: reordvar(vl,v,&nvl);
! 568: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 569:
! 570: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 571: *pr = 0;
! 572: return;
! 573: }
! 574: if ( VR(q1) != v ) {
! 575: pwrp(vl,q1,DEG(DC(q2)),pr);
! 576: return;
! 577: }
! 578: if ( VR(q2) != v ) {
! 579: pwrp(vl,q2,DEG(DC(q1)),pr);
! 580: return;
! 581: }
! 582: norm1c(q1,&a); norm1c(q2,&b);
! 583: n = DEG(DC(q1)); m = DEG(DC(q2));
! 584: pwrq(a,(Q)m,&w); pwrq(b,(Q)n,&s); mulq(w,s,&u);
! 585: factorialz(QTOS(n)+QTOS(m),&t);
! 586: mulq(u,(Q)t,&s); addq(s,s,&f);
! 587: for ( index = 0, q = (Q)ONE, c = 0; cmpq(f,q) >= 0; ) {
! 588: mod = get_lprime(index++);
! 589: ptomp(mod,LC(q1),&tg);
! 590: if ( !tg )
! 591: continue;
! 592: ptomp(mod,LC(q2),&tg);
! 593: if ( !tg )
! 594: continue;
! 595: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 596: srchmp(nvl,mod,v,tg1,tg2,&resg);
! 597: chnremp(nvl,mod,c,(Z)q,resg,&c1); c = c1;
! 598: STOQ(mod,t); mulq(q,(Q)t,&s); q = s;
! 599: }
! 600: *pr = c;
! 601: }
! 602:
! 603: void res_ch_det(VL vl,V v,P p1,P p2,P *pr)
! 604: {
! 605: P q1,q2,c,c1;
! 606: P tg,tg1,tg2,resg;
! 607: int index,mod;
! 608: Q a,b,f,q,s,u,w;
! 609: Z n,m,t;
! 610: VL nvl;
! 611:
! 612: reordvar(vl,v,&nvl);
! 613: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 614:
! 615: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 616: *pr = 0;
! 617: return;
! 618: }
! 619: if ( VR(q1) != v ) {
! 620: pwrp(vl,q1,DEG(DC(q2)),pr);
! 621: return;
! 622: }
! 623: if ( VR(q2) != v ) {
! 624: pwrp(vl,q2,DEG(DC(q1)),pr);
! 625: return;
! 626: }
! 627: norm1c(q1,&a); norm1c(q2,&b);
! 628: n = DEG(DC(q1)); m = DEG(DC(q2));
! 629: pwrq(a,(Q)m,&w); pwrq(b,(Q)n,&s); mulq(w,s,&u);
! 630: factorialz(QTOS(n)+QTOS(m),&t);
! 631: mulq(u,(Q)t,&s); addq(s,s,&f);
! 632: for ( index = 0, q = (Q)ONE, c = 0; cmpq(f,q) >= 0; ) {
! 633: mod = get_lprime(index++);
! 634: ptomp(mod,LC(q1),&tg);
! 635: if ( !tg )
! 636: continue;
! 637: ptomp(mod,LC(q2),&tg);
! 638: if ( !tg )
! 639: continue;
! 640: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 641: res_detmp(nvl,mod,v,tg1,tg2,&resg);
! 642: chnremp(nvl,mod,c,(Z)q,resg,&c1); c = c1;
! 643: STOQ(mod,t); mulq(q,(Q)t,&s); q = s;
! 644: }
! 645: *pr = c;
! 646: }
! 647:
! 648: void res_detmp(VL vl,int mod,V v,P p1,P p2,P *dp)
! 649: {
! 650: int n1,n2,n,sgn;
! 651: int i,j,k;
! 652: P mjj,mij,t,s,u,d;
! 653: P **mat;
! 654: P *mi,*mj;
! 655: DCP dc;
! 656:
! 657: n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
! 658: mat = (P **)almat_pointer(n,n);
! 659: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 660: mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
! 661: for ( i = 1; i < n2; i++ )
! 662: for ( j = 0; j <= n1; j++ )
! 663: mat[i][i+j] = mat[0][j];
! 664: for ( dc = DC(p2); dc; dc = NEXT(dc) )
! 665: mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
! 666: for ( i = 1; i < n1; i++ )
! 667: for ( j = 0; j <= n2; j++ )
! 668: mat[i+n2][i+j] = mat[n2][j];
! 669: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
! 670: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 671: if ( i == n ) {
! 672: *dp = 0; return;
! 673: }
! 674: for ( k = i; k < n; k++ )
! 675: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 676: i = k;
! 677: if ( j != i ) {
! 678: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
! 679: }
! 680: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 681: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 682: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
! 683: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
! 684: }
! 685: d = mjj;
! 686: }
! 687: if ( sgn > 0 )
! 688: *dp = d;
! 689: else
! 690: chsgnmp(mod,d,dp);
! 691: }
! 692:
! 693: #if 0
! 694: showmat(VL vl,P **mat,int n)
! 695: {
! 696: int i,j;
! 697: P t;
! 698:
! 699: for ( i = 0; i < n; i++ ) {
! 700: for ( j = 0; j < n; j++ ) {
! 701: mptop(mat[i][j],&t); asir_printp(vl,t); fprintf(out," ");
! 702: }
! 703: fprintf(out,"\n");
! 704: }
! 705: fflush(out);
! 706: }
! 707:
! 708: showmp(VL vl,P p)
! 709: {
! 710: P t;
! 711:
! 712: mptop(p,&t); asir_printp(vl,t); fprintf(out,"\n");
! 713: }
! 714: #endif
! 715:
! 716: void premp(VL vl,P p1,P p2,P *pr)
! 717: {
! 718: P m,m1,m2;
! 719: P *pw;
! 720: DCP dc;
! 721: V v1,v2;
! 722: register int i,j;
! 723: int n1,n2,d;
! 724:
! 725: if ( NUM(p1) )
! 726: if ( NUM(p2) )
! 727: *pr = 0;
! 728: else
! 729: *pr = p1;
! 730: else if ( NUM(p2) )
! 731: *pr = 0;
! 732: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 733: n1 = deg(v1,p1); n2 = deg(v1,p2);
! 734: if ( n1 < n2 ) {
! 735: *pr = p1;
! 736: return;
! 737: }
! 738: pw = (P *)ALLOCA((n1+1)*sizeof(P));
! 739: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
! 740:
! 741: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 742: pw[QTOS(DEG(dc))] = COEF(dc);
! 743:
! 744: for ( i = n1; i >= n2; i-- ) {
! 745: if ( pw[i] ) {
! 746: m = pw[i];
! 747: for ( j = i; j >= 0; j-- ) {
! 748: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 749: }
! 750:
! 751: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
! 752: mulp(vl,COEF(dc),m,&m1);
! 753: subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
! 754: pw[QTOS(DEG(dc))+d] = m2;
! 755: }
! 756: } else
! 757: for ( j = i; j >= 0; j-- )
! 758: if ( pw[j] ) {
! 759: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 760: }
! 761: }
! 762: plisttop(pw,v1,n2-1,pr);
! 763: } else {
! 764: while ( v1 != vl->v && v2 != vl->v )
! 765: vl = NEXT(vl);
! 766: if ( v1 == vl->v )
! 767: *pr = 0;
! 768: else
! 769: *pr = p1;
! 770: }
! 771: }
! 772:
! 773: void ptozp0(P p,P *pr)
! 774: {
! 775: Q c;
! 776:
! 777: if ( qpcheck((Obj)p) )
! 778: ptozp(p,1,&c,pr);
! 779: else
! 780: *pr = p;
! 781: }
! 782:
! 783: void mindegp(VL vl,P p,VL *mvlp,P *pr)
! 784: {
! 785: P t;
! 786: VL nvl,tvl,avl;
! 787: V v;
! 788: int n,d;
! 789:
! 790: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 791: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 792: n = getdeg(avl->v,p);
! 793: if ( n < d ) {
! 794: v = avl->v; d = n;
! 795: }
! 796: }
! 797: if ( v != nvl->v ) {
! 798: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 799: *pr = t; *mvlp = tvl;
! 800: } else {
! 801: *pr = p; *mvlp = nvl;
! 802: }
! 803: }
! 804:
! 805: void maxdegp(VL vl,P p,VL *mvlp,P *pr)
! 806: {
! 807: P t;
! 808: VL nvl,tvl,avl;
! 809: V v;
! 810: int n,d;
! 811:
! 812: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 813: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 814: n = getdeg(avl->v,p);
! 815: if ( n > d ) {
! 816: v = avl->v; d = n;
! 817: }
! 818: }
! 819: if ( v != nvl->v ) {
! 820: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 821: *pr = t; *mvlp = tvl;
! 822: } else {
! 823: *pr = p; *mvlp = nvl;
! 824: }
! 825: }
! 826:
! 827: void min_common_vars_in_coefp(VL vl,P p,VL *mvlp,P *pr)
! 828: {
! 829: P u,p0;
! 830: VL tvl,cvl,svl,uvl,avl,vl0;
! 831: int n,n0,d,d0;
! 832: DCP dc;
! 833:
! 834: clctv(vl,p,&tvl); vl = tvl;
! 835: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
! 836: for ( avl = vl; avl; avl = NEXT(avl) ) {
! 837: if ( VR(p) != avl->v ) {
! 838: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 839: } else {
! 840: uvl = vl; u = p;
! 841: }
! 842: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
! 843: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
! 844: }
! 845: if ( !cvl ) {
! 846: *mvlp = uvl; *pr = u; return;
! 847: } else {
! 848: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
! 849: if ( n < n0 ) {
! 850: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
! 851: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
! 852: vl0 = uvl; p0 = u; n0 = n; d0 = d;
! 853: }
! 854: }
! 855: }
! 856: *pr = p0; *mvlp = vl0;
! 857: }
! 858:
! 859: void minlcdegp(VL vl,P p,VL *mvlp,P *pr)
! 860: {
! 861: P u,p0;
! 862: VL tvl,uvl,avl,vl0;
! 863: int d,d0;
! 864:
! 865: clctv(vl,p,&tvl); vl = tvl;
! 866: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
! 867: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
! 868: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 869: d = homdeg(COEF(DC(u)));
! 870: if ( d < d0 ) {
! 871: vl0 = uvl; p0 = u; d0 = d;
! 872: }
! 873: }
! 874: *pr = p0; *mvlp = vl0;
! 875: }
! 876:
! 877: void sort_by_deg(int n,P *p,P *pr)
! 878: {
! 879: int j,k,d,k0;
! 880: V v;
! 881:
! 882: for ( j = 0; j < n; j++ ) {
! 883: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 884: k < n; k++ )
! 885: if ( deg(v,p[k]) < d ) {
! 886: k0 = k;
! 887: d = deg(v,p[k]);
! 888: }
! 889: pr[j] = p[k0];
! 890: if ( j != k0 )
! 891: p[k0] = p[j];
! 892: }
! 893: }
! 894:
! 895: void sort_by_deg_rev(int n,P *p,P *pr)
! 896: {
! 897: int j,k,d,k0;
! 898: V v;
! 899:
! 900: for ( j = 0; j < n; j++ ) {
! 901: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 902: k < n; k++ )
! 903: if ( deg(v,p[k]) > d ) {
! 904: k0 = k;
! 905: d = deg(v,p[k]);
! 906: }
! 907: pr[j] = p[k0];
! 908: if ( j != k0 )
! 909: p[k0] = p[j];
! 910: }
! 911: }
! 912:
! 913:
! 914: void getmindeg(V v,P p,Z *dp)
! 915: {
! 916: Z dt,d;
! 917: DCP dc;
! 918:
! 919: if ( !p || NUM(p) )
! 920: *dp = 0;
! 921: else if ( v == VR(p) ) {
! 922: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
! 923: *dp = DEG(dc);
! 924: } else {
! 925: dc = DC(p);
! 926: getmindeg(v,COEF(dc),&d);
! 927: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 928: getmindeg(v,COEF(dc),&dt);
! 929: if ( cmpz(dt,d) < 0 )
! 930: d = dt;
! 931: }
! 932: *dp = d;
! 933: }
! 934: }
! 935:
! 936: void minchdegp(VL vl,P p,VL *mvlp,P *pr)
! 937: {
! 938: P t;
! 939: VL tvl,nvl,avl;
! 940: int n,d,z;
! 941: V v;
! 942:
! 943: if ( NUM(p) ) {
! 944: *mvlp = vl;
! 945: *pr = p;
! 946: return;
! 947: }
! 948: clctv(vl,p,&nvl);
! 949: v = nvl->v;
! 950: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
! 951: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 952: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
! 953: if ( n < d ) {
! 954: v = avl->v; d = n;
! 955: }
! 956: }
! 957: if ( v != nvl->v ) {
! 958: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 959: *pr = t; *mvlp = tvl;
! 960: } else {
! 961: *pr = p; *mvlp = nvl;
! 962: }
! 963: }
! 964:
! 965: int getchomdeg(V v,P p)
! 966: {
! 967: int m,m1;
! 968: DCP dc;
! 969:
! 970: if ( !p || NUM(p) )
! 971: return ( 0 );
! 972: else if ( VR(p) == v )
! 973: return ( dbound(v,p) );
! 974: else {
! 975: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 976: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
! 977: m = MAX(m,m1);
! 978: }
! 979: return ( m );
! 980: }
! 981: }
! 982:
! 983: int getlchomdeg(V v,P p,int *d)
! 984: {
! 985: int m0,m1,d0,d1;
! 986: DCP dc;
! 987:
! 988: if ( !p || NUM(p) ) {
! 989: *d = 0;
! 990: return ( 0 );
! 991: } else if ( VR(p) == v ) {
! 992: *d = QTOS(DEG(DC(p)));
! 993: return ( homdeg(LC(p)) );
! 994: } else {
! 995: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
! 996: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
! 997: if ( d1 > d0 ) {
! 998: m0 = m1;
! 999: d0 = d1;
! 1000: } else if ( d1 == d0 )
! 1001: m0 = MAX(m0,m1);
! 1002: }
! 1003: *d = d0;
! 1004: return ( m0 );
! 1005: }
! 1006: }
! 1007:
! 1008: int nmonop(P p)
! 1009: {
! 1010: int s;
! 1011: DCP dc;
! 1012:
! 1013: if ( !p )
! 1014: return 0;
! 1015: else
! 1016: switch ( OID(p) ) {
! 1017: case O_N:
! 1018: return 1; break;
! 1019: case O_P:
! 1020: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
! 1021: s += nmonop(COEF(dc));
! 1022: return s; break;
! 1023: default:
! 1024: return 0;
! 1025: }
! 1026: }
! 1027:
! 1028: int qpcheck(Obj p)
! 1029: {
! 1030: DCP dc;
! 1031:
! 1032: if ( !p )
! 1033: return 1;
! 1034: else
! 1035: switch ( OID(p) ) {
! 1036: case O_N:
! 1037: return RATN((Num)p)?1:0;
! 1038: case O_P:
! 1039: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
! 1040: if ( !qpcheck((Obj)COEF(dc)) )
! 1041: return 0;
! 1042: return 1;
! 1043: default:
! 1044: return 0;
! 1045: }
! 1046: }
! 1047:
! 1048: /* check if p is univariate and all coeffs are INT or LM */
! 1049:
! 1050: int uzpcheck(Obj p)
! 1051: {
! 1052: DCP dc;
! 1053: P c;
! 1054:
! 1055: if ( !p )
! 1056: return 1;
! 1057: else
! 1058: switch ( OID(p) ) {
! 1059: case O_N:
! 1060: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
! 1061: case O_P:
! 1062: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
! 1063: c = COEF(dc);
! 1064: if ( !NUM(c) || !uzpcheck((Obj)c) )
! 1065: return 0;
! 1066: }
! 1067: return 1;
! 1068: default:
! 1069: return 0;
! 1070: }
! 1071: }
! 1072:
! 1073: int p_mag(P p)
! 1074: {
! 1075: int s;
! 1076: DCP dc;
! 1077:
! 1078: if ( !p )
! 1079: return 0;
! 1080: else if ( OID(p) == O_N )
! 1081: return z_bits((Q)p);
! 1082: else {
! 1083: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
! 1084: s += p_mag(COEF(dc));
! 1085: return s;
! 1086: }
! 1087: }
! 1088:
! 1089: int maxblenp(P p)
! 1090: {
! 1091: int s,t;
! 1092: DCP dc;
! 1093:
! 1094: if ( !p )
! 1095: return 0;
! 1096: else if ( OID(p) == O_N )
! 1097: return z_bits((Q)p);
! 1098: else {
! 1099: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 1100: t = maxblenp(COEF(dc));
! 1101: s = MAX(t,s);
! 1102: }
! 1103: return s;
! 1104: }
! 1105: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PU.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine