Annotation of OpenXM_contrib2/asir2000/engine/PU.c, Revision 1.16
1.3 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.4 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.3 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.16 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.15 2018/02/28 03:55:38 noro Exp $
1.3 noro 49: */
1.1 noro 50: #include "ca.h"
51:
1.8 noro 52: void reorderp(VL nvl,VL ovl,P p,P *pr)
1.1 noro 53: {
1.16 ! noro 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: }
1.1 noro 74: }
1.16 ! noro 75:
1.8 noro 76: void substp(VL vl,P p,V v0,P p0,P *pr)
1.1 noro 77: {
1.16 ! noro 78: P x,t,m,c,s,a;
! 79: DCP dc;
! 80: Q 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: subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)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: }
1.1 noro 115: }
116:
1.13 noro 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: {
1.16 ! noro 121: P x,t,m,c,s,a,p0,c1;
! 122: DCP dc;
! 123: Q 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: subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)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: }
1.13 noro 166: }
167:
1.8 noro 168: void detp(VL vl,P **rmat,int n,P *dp)
1.1 noro 169: {
1.16 ! noro 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;
1.7 noro 211: }
212:
1.8 noro 213: void invmatp(VL vl,P **rmat,int n,P ***imatp,P *dnp)
1.7 noro 214: {
1.16 ! noro 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;
1.1 noro 264: }
265:
1.8 noro 266: void reordvar(VL vl,V v,VL *nvlp)
1.1 noro 267: {
1.16 ! noro 268: VL nvl,nvl0;
1.1 noro 269:
1.16 ! noro 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;
1.1 noro 281: }
282:
1.8 noro 283: void gcdprsp(VL vl,P p1,P p2,P *pr)
1.1 noro 284: {
1.16 ! noro 285: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
! 286: V v1,v2;
1.1 noro 287:
1.16 ! noro 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: }
1.1 noro 319: }
320:
1.8 noro 321: void gcdcp(VL vl,P p,P *pr)
1.1 noro 322: {
1.16 ! noro 323: P g,g1;
! 324: DCP dc;
1.1 noro 325:
1.16 ! noro 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: }
1.1 noro 337: }
338:
1.8 noro 339: void sprs(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 340: {
1.16 ! noro 341: P q1,q2,m,m1,m2,x,h,r,g1,g2;
! 342: int d;
! 343: Q 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;
1.1 noro 375: }
376:
1.8 noro 377: void resultp(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 378: {
1.16 ! noro 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: Q 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: }
1.1 noro 467: }
468:
1.8 noro 469: void srch2(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 470: {
1.16 ! noro 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: Q 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: }
1.1 noro 556: }
557:
1.8 noro 558: void srcr(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 559: {
1.16 ! noro 560: P q1,q2,c,c1;
! 561: P tg,tg1,tg2,resg;
! 562: int index,mod;
! 563: Q a,b,f,q,t,s,u,n,m;
! 564: VL nvl;
! 565:
! 566: reordvar(vl,v,&nvl);
! 567: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 568:
! 569: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 570: *pr = 0;
! 571: return;
! 572: }
! 573: if ( VR(q1) != v ) {
! 574: pwrp(vl,q1,DEG(DC(q2)),pr);
! 575: return;
! 576: }
! 577: if ( VR(q2) != v ) {
! 578: pwrp(vl,q2,DEG(DC(q1)),pr);
! 579: return;
! 580: }
! 581: norm1c(q1,&a); norm1c(q2,&b);
! 582: n = DEG(DC(q1)); m = DEG(DC(q2));
! 583: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
! 584: factorial(QTOS(n)+QTOS(m),&t);
! 585: mulq(u,t,&s); addq(s,s,&f);
! 586: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
! 587: mod = get_lprime(index++);
! 588: ptomp(mod,LC(q1),&tg);
! 589: if ( !tg )
! 590: continue;
! 591: ptomp(mod,LC(q2),&tg);
! 592: if ( !tg )
! 593: continue;
! 594: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 595: srchmp(nvl,mod,v,tg1,tg2,&resg);
! 596: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
! 597: STOQ(mod,t); mulq(q,t,&s); q = s;
! 598: }
! 599: *pr = c;
1.1 noro 600: }
601:
1.8 noro 602: void res_ch_det(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 603: {
1.16 ! noro 604: P q1,q2,c,c1;
! 605: P tg,tg1,tg2,resg;
! 606: int index,mod;
! 607: Q a,b,f,q,t,s,u,n,m;
! 608: VL nvl;
! 609:
! 610: reordvar(vl,v,&nvl);
! 611: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 612:
! 613: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 614: *pr = 0;
! 615: return;
! 616: }
! 617: if ( VR(q1) != v ) {
! 618: pwrp(vl,q1,DEG(DC(q2)),pr);
! 619: return;
! 620: }
! 621: if ( VR(q2) != v ) {
! 622: pwrp(vl,q2,DEG(DC(q1)),pr);
! 623: return;
! 624: }
! 625: norm1c(q1,&a); norm1c(q2,&b);
! 626: n = DEG(DC(q1)); m = DEG(DC(q2));
! 627: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
! 628: factorial(QTOS(n)+QTOS(m),&t);
! 629: mulq(u,t,&s); addq(s,s,&f);
! 630: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
! 631: mod = get_lprime(index++);
! 632: ptomp(mod,LC(q1),&tg);
! 633: if ( !tg )
! 634: continue;
! 635: ptomp(mod,LC(q2),&tg);
! 636: if ( !tg )
! 637: continue;
! 638: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 639: res_detmp(nvl,mod,v,tg1,tg2,&resg);
! 640: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
! 641: STOQ(mod,t); mulq(q,t,&s); q = s;
! 642: }
! 643: *pr = c;
1.1 noro 644: }
645:
1.8 noro 646: void res_detmp(VL vl,int mod,V v,P p1,P p2,P *dp)
1.1 noro 647: {
1.16 ! noro 648: int n1,n2,n,sgn;
! 649: int i,j,k;
! 650: P mjj,mij,t,s,u,d;
! 651: P **mat;
! 652: P *mi,*mj;
! 653: DCP dc;
! 654:
! 655: n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
! 656: mat = (P **)almat_pointer(n,n);
! 657: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 658: mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
! 659: for ( i = 1; i < n2; i++ )
! 660: for ( j = 0; j <= n1; j++ )
! 661: mat[i][i+j] = mat[0][j];
! 662: for ( dc = DC(p2); dc; dc = NEXT(dc) )
! 663: mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
! 664: for ( i = 1; i < n1; i++ )
! 665: for ( j = 0; j <= n2; j++ )
! 666: mat[i+n2][i+j] = mat[n2][j];
! 667: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
! 668: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 669: if ( i == n ) {
! 670: *dp = 0; return;
! 671: }
! 672: for ( k = i; k < n; k++ )
! 673: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 674: i = k;
! 675: if ( j != i ) {
! 676: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
! 677: }
! 678: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 679: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 680: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
! 681: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
! 682: }
! 683: d = mjj;
! 684: }
! 685: if ( sgn > 0 )
! 686: *dp = d;
! 687: else
! 688: chsgnmp(mod,d,dp);
1.1 noro 689: }
690:
691: #if 0
1.8 noro 692: showmat(VL vl,P **mat,int n)
1.1 noro 693: {
1.16 ! noro 694: int i,j;
! 695: P t;
1.1 noro 696:
1.16 ! noro 697: for ( i = 0; i < n; i++ ) {
! 698: for ( j = 0; j < n; j++ ) {
! 699: mptop(mat[i][j],&t); asir_printp(vl,t); fprintf(out," ");
! 700: }
! 701: fprintf(out,"\n");
! 702: }
! 703: fflush(out);
1.1 noro 704: }
705:
1.8 noro 706: showmp(VL vl,P p)
1.1 noro 707: {
1.16 ! noro 708: P t;
1.1 noro 709:
1.16 ! noro 710: mptop(p,&t); asir_printp(vl,t); fprintf(out,"\n");
1.1 noro 711: }
712: #endif
713:
1.8 noro 714: void premp(VL vl,P p1,P p2,P *pr)
1.1 noro 715: {
1.16 ! noro 716: P m,m1,m2;
! 717: P *pw;
! 718: DCP dc;
! 719: V v1,v2;
! 720: register int i,j;
! 721: int n1,n2,d;
! 722:
! 723: if ( NUM(p1) )
! 724: if ( NUM(p2) )
! 725: *pr = 0;
! 726: else
! 727: *pr = p1;
! 728: else if ( NUM(p2) )
! 729: *pr = 0;
! 730: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 731: n1 = deg(v1,p1); n2 = deg(v1,p2);
1.15 noro 732: if ( n1 < n2 ) {
733: *pr = p1;
734: return;
735: }
1.16 ! noro 736: pw = (P *)ALLOCA((n1+1)*sizeof(P));
! 737: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
1.1 noro 738:
1.16 ! noro 739: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 740: pw[QTOS(DEG(dc))] = COEF(dc);
1.1 noro 741:
1.16 ! noro 742: for ( i = n1; i >= n2; i-- ) {
! 743: if ( pw[i] ) {
! 744: m = pw[i];
! 745: for ( j = i; j >= 0; j-- ) {
! 746: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 747: }
! 748:
! 749: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
! 750: mulp(vl,COEF(dc),m,&m1);
! 751: subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
! 752: pw[QTOS(DEG(dc))+d] = m2;
! 753: }
! 754: } else
! 755: for ( j = i; j >= 0; j-- )
! 756: if ( pw[j] ) {
! 757: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 758: }
! 759: }
! 760: plisttop(pw,v1,n2-1,pr);
! 761: } else {
! 762: while ( v1 != vl->v && v2 != vl->v )
! 763: vl = NEXT(vl);
! 764: if ( v1 == vl->v )
! 765: *pr = 0;
! 766: else
! 767: *pr = p1;
! 768: }
1.1 noro 769: }
770:
1.8 noro 771: void ptozp0(P p,P *pr)
1.1 noro 772: {
1.16 ! noro 773: Q c;
1.1 noro 774:
1.16 ! noro 775: if ( qpcheck((Obj)p) )
! 776: ptozp(p,1,&c,pr);
! 777: else
! 778: *pr = p;
1.1 noro 779: }
780:
1.8 noro 781: void mindegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 782: {
1.16 ! noro 783: P t;
! 784: VL nvl,tvl,avl;
! 785: V v;
! 786: int n,d;
! 787:
! 788: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 789: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 790: n = getdeg(avl->v,p);
! 791: if ( n < d ) {
! 792: v = avl->v; d = n;
! 793: }
! 794: }
! 795: if ( v != nvl->v ) {
! 796: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 797: *pr = t; *mvlp = tvl;
! 798: } else {
! 799: *pr = p; *mvlp = nvl;
! 800: }
1.1 noro 801: }
802:
1.8 noro 803: void maxdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 804: {
1.16 ! noro 805: P t;
! 806: VL nvl,tvl,avl;
! 807: V v;
! 808: int n,d;
! 809:
! 810: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 811: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 812: n = getdeg(avl->v,p);
! 813: if ( n > d ) {
! 814: v = avl->v; d = n;
! 815: }
! 816: }
! 817: if ( v != nvl->v ) {
! 818: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 819: *pr = t; *mvlp = tvl;
! 820: } else {
! 821: *pr = p; *mvlp = nvl;
! 822: }
1.1 noro 823: }
824:
1.8 noro 825: void min_common_vars_in_coefp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 826: {
1.16 ! noro 827: P u,p0;
! 828: VL tvl,cvl,svl,uvl,avl,vl0;
! 829: int n,n0,d,d0;
! 830: DCP dc;
! 831:
! 832: clctv(vl,p,&tvl); vl = tvl;
! 833: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
! 834: for ( avl = vl; avl; avl = NEXT(avl) ) {
! 835: if ( VR(p) != avl->v ) {
! 836: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 837: } else {
! 838: uvl = vl; u = p;
! 839: }
! 840: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
! 841: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
! 842: }
! 843: if ( !cvl ) {
! 844: *mvlp = uvl; *pr = u; return;
! 845: } else {
! 846: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
! 847: if ( n < n0 ) {
! 848: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
! 849: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
! 850: vl0 = uvl; p0 = u; n0 = n; d0 = d;
! 851: }
! 852: }
! 853: }
! 854: *pr = p0; *mvlp = vl0;
1.1 noro 855: }
856:
1.8 noro 857: void minlcdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 858: {
1.16 ! noro 859: P u,p0;
! 860: VL tvl,uvl,avl,vl0;
! 861: int d,d0;
! 862:
! 863: clctv(vl,p,&tvl); vl = tvl;
! 864: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
! 865: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
! 866: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 867: d = homdeg(COEF(DC(u)));
! 868: if ( d < d0 ) {
! 869: vl0 = uvl; p0 = u; d0 = d;
! 870: }
! 871: }
! 872: *pr = p0; *mvlp = vl0;
1.1 noro 873: }
874:
1.8 noro 875: void sort_by_deg(int n,P *p,P *pr)
1.1 noro 876: {
1.16 ! noro 877: int j,k,d,k0;
! 878: V v;
1.1 noro 879:
1.16 ! noro 880: for ( j = 0; j < n; j++ ) {
! 881: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 882: k < n; k++ )
! 883: if ( deg(v,p[k]) < d ) {
! 884: k0 = k;
! 885: d = deg(v,p[k]);
! 886: }
! 887: pr[j] = p[k0];
! 888: if ( j != k0 )
! 889: p[k0] = p[j];
! 890: }
1.1 noro 891: }
892:
1.8 noro 893: void sort_by_deg_rev(int n,P *p,P *pr)
1.1 noro 894: {
1.16 ! noro 895: int j,k,d,k0;
! 896: V v;
1.1 noro 897:
1.16 ! noro 898: for ( j = 0; j < n; j++ ) {
! 899: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 900: k < n; k++ )
! 901: if ( deg(v,p[k]) > d ) {
! 902: k0 = k;
! 903: d = deg(v,p[k]);
! 904: }
! 905: pr[j] = p[k0];
! 906: if ( j != k0 )
! 907: p[k0] = p[j];
! 908: }
1.1 noro 909: }
910:
911:
1.8 noro 912: void getmindeg(V v,P p,Q *dp)
1.1 noro 913: {
1.16 ! noro 914: Q dt,d;
! 915: DCP dc;
1.1 noro 916:
1.16 ! noro 917: if ( !p || NUM(p) )
! 918: *dp = 0;
! 919: else if ( v == VR(p) ) {
! 920: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
! 921: *dp = DEG(dc);
! 922: } else {
! 923: dc = DC(p);
! 924: getmindeg(v,COEF(dc),&d);
! 925: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 926: getmindeg(v,COEF(dc),&dt);
! 927: if ( cmpq(dt,d) < 0 )
! 928: d = dt;
! 929: }
! 930: *dp = d;
! 931: }
1.1 noro 932: }
933:
1.8 noro 934: void minchdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 935: {
1.16 ! noro 936: P t;
! 937: VL tvl,nvl,avl;
! 938: int n,d,z;
! 939: V v;
! 940:
! 941: if ( NUM(p) ) {
! 942: *mvlp = vl;
! 943: *pr = p;
! 944: return;
! 945: }
! 946: clctv(vl,p,&nvl);
! 947: v = nvl->v;
! 948: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
! 949: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 950: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
! 951: if ( n < d ) {
! 952: v = avl->v; d = n;
! 953: }
! 954: }
! 955: if ( v != nvl->v ) {
! 956: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 957: *pr = t; *mvlp = tvl;
! 958: } else {
! 959: *pr = p; *mvlp = nvl;
! 960: }
1.1 noro 961: }
962:
1.8 noro 963: int getchomdeg(V v,P p)
1.1 noro 964: {
1.16 ! noro 965: int m,m1;
! 966: DCP dc;
1.1 noro 967:
1.16 ! noro 968: if ( !p || NUM(p) )
! 969: return ( 0 );
! 970: else if ( VR(p) == v )
! 971: return ( dbound(v,p) );
! 972: else {
! 973: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 974: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
! 975: m = MAX(m,m1);
! 976: }
! 977: return ( m );
! 978: }
1.1 noro 979: }
980:
1.8 noro 981: int getlchomdeg(V v,P p,int *d)
1.1 noro 982: {
1.16 ! noro 983: int m0,m1,d0,d1;
! 984: DCP dc;
1.1 noro 985:
1.16 ! noro 986: if ( !p || NUM(p) ) {
! 987: *d = 0;
! 988: return ( 0 );
! 989: } else if ( VR(p) == v ) {
! 990: *d = QTOS(DEG(DC(p)));
! 991: return ( homdeg(LC(p)) );
! 992: } else {
! 993: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
! 994: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
! 995: if ( d1 > d0 ) {
! 996: m0 = m1;
! 997: d0 = d1;
! 998: } else if ( d1 == d0 )
! 999: m0 = MAX(m0,m1);
! 1000: }
! 1001: *d = d0;
! 1002: return ( m0 );
! 1003: }
1.1 noro 1004: }
1005:
1.8 noro 1006: int nmonop(P p)
1.1 noro 1007: {
1.16 ! noro 1008: int s;
! 1009: DCP dc;
1.1 noro 1010:
1.16 ! noro 1011: if ( !p )
! 1012: return 0;
! 1013: else
! 1014: switch ( OID(p) ) {
! 1015: case O_N:
! 1016: return 1; break;
! 1017: case O_P:
! 1018: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
! 1019: s += nmonop(COEF(dc));
! 1020: return s; break;
! 1021: default:
! 1022: return 0;
! 1023: }
1.1 noro 1024: }
1025:
1.8 noro 1026: int qpcheck(Obj p)
1.1 noro 1027: {
1.16 ! noro 1028: DCP dc;
1.1 noro 1029:
1.16 ! noro 1030: if ( !p )
! 1031: return 1;
! 1032: else
! 1033: switch ( OID(p) ) {
! 1034: case O_N:
! 1035: return RATN((Num)p)?1:0;
! 1036: case O_P:
! 1037: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
! 1038: if ( !qpcheck((Obj)COEF(dc)) )
! 1039: return 0;
! 1040: return 1;
! 1041: default:
! 1042: return 0;
! 1043: }
1.1 noro 1044: }
1045:
1046: /* check if p is univariate and all coeffs are INT or LM */
1047:
1.8 noro 1048: int uzpcheck(Obj p)
1.1 noro 1049: {
1.16 ! noro 1050: DCP dc;
! 1051: P c;
1.1 noro 1052:
1.16 ! noro 1053: if ( !p )
! 1054: return 1;
! 1055: else
! 1056: switch ( OID(p) ) {
! 1057: case O_N:
! 1058: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
! 1059: case O_P:
! 1060: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
! 1061: c = COEF(dc);
! 1062: if ( !NUM(c) || !uzpcheck((Obj)c) )
! 1063: return 0;
! 1064: }
! 1065: return 1;
! 1066: default:
! 1067: return 0;
! 1068: }
1.1 noro 1069: }
1070:
1.8 noro 1071: int p_mag(P p)
1.1 noro 1072: {
1.16 ! noro 1073: int s;
! 1074: DCP dc;
1.1 noro 1075:
1.16 ! noro 1076: if ( !p )
! 1077: return 0;
! 1078: else if ( OID(p) == O_N )
! 1079: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
! 1080: else {
! 1081: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
! 1082: s += p_mag(COEF(dc));
! 1083: return s;
! 1084: }
1.1 noro 1085: }
1.2 noro 1086:
1.8 noro 1087: int maxblenp(P p)
1.2 noro 1088: {
1.16 ! noro 1089: int s,t;
! 1090: DCP dc;
1.2 noro 1091:
1.16 ! noro 1092: if ( !p )
! 1093: return 0;
! 1094: else if ( OID(p) == O_N )
! 1095: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
! 1096: else {
! 1097: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 1098: t = maxblenp(COEF(dc));
! 1099: s = MAX(t,s);
! 1100: }
! 1101: return s;
! 1102: }
1.2 noro 1103: }
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 / asir2000 / engine