![[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|
Return to PU.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
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.15 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.14 2010/01/31 03:25:54 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: {
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:
1.8 noro 76: void substp(VL vl,P p,V v0,P p0,P *pr)
1.1 noro 77: {
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: }
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: {
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;
1.14 noro 134: if ( svect[i] && OID(svect[i]) < 0 ) {
1.13 noro 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: }
166: }
167:
1.8 noro 168: void detp(VL vl,P **rmat,int n,P *dp)
1.1 noro 169: {
1.11 noro 170: int i,j,k,l,sgn,nmin,kmin,lmin,ntmp;
1.1 noro 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: }
1.11 noro 184: nmin = nmonop(mat[i][j]);
1.12 noro 185: kmin=i; lmin=j;
1.11 noro 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;
1.1 noro 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: {
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: {
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:
1.8 noro 283: void gcdprsp(VL vl,P p1,P p2,P *pr)
1.1 noro 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:
1.8 noro 321: void gcdcp(VL vl,P p,P *pr)
1.1 noro 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:
1.8 noro 339: void sprs(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 340: {
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;
375: }
376:
1.8 noro 377: void resultp(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 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: 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: }
467: }
468:
1.8 noro 469: void srch2(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 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: 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: }
556: }
557:
1.8 noro 558: void srcr(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 559: {
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; ) {
1.6 noro 587: mod = get_lprime(index++);
1.1 noro 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;
600: }
601:
1.8 noro 602: void res_ch_det(VL vl,V v,P p1,P p2,P *pr)
1.1 noro 603: {
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; ) {
1.6 noro 631: mod = get_lprime(index++);
1.1 noro 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;
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: {
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);
689: }
690:
691: #if 0
1.8 noro 692: showmat(VL vl,P **mat,int n)
1.1 noro 693: {
694: int i,j;
695: P t;
696:
697: for ( i = 0; i < n; i++ ) {
698: for ( j = 0; j < n; j++ ) {
1.9 noro 699: mptop(mat[i][j],&t); asir_printp(vl,t); fprintf(out," ");
1.1 noro 700: }
701: fprintf(out,"\n");
702: }
703: fflush(out);
704: }
705:
1.8 noro 706: showmp(VL vl,P p)
1.1 noro 707: {
708: P t;
709:
1.9 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: {
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.1 noro 736: pw = (P *)ALLOCA((n1+1)*sizeof(P));
737: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
738:
739: for ( dc = DC(p1); dc; dc = NEXT(dc) )
740: pw[QTOS(DEG(dc))] = COEF(dc);
741:
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: }
769: }
770:
1.8 noro 771: void ptozp0(P p,P *pr)
1.1 noro 772: {
773: Q c;
774:
1.5 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: {
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: }
801: }
802:
1.8 noro 803: void maxdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 804: {
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: }
823: }
824:
1.8 noro 825: void min_common_vars_in_coefp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 826: {
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;
855: }
856:
1.8 noro 857: void minlcdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 858: {
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;
873: }
874:
1.8 noro 875: void sort_by_deg(int n,P *p,P *pr)
1.1 noro 876: {
877: int j,k,d,k0;
878: V v;
879:
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: }
891: }
892:
1.8 noro 893: void sort_by_deg_rev(int n,P *p,P *pr)
1.1 noro 894: {
895: int j,k,d,k0;
896: V v;
897:
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: }
909: }
910:
911:
1.8 noro 912: void getmindeg(V v,P p,Q *dp)
1.1 noro 913: {
914: Q dt,d;
915: DCP dc;
916:
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: }
932: }
933:
1.8 noro 934: void minchdegp(VL vl,P p,VL *mvlp,P *pr)
1.1 noro 935: {
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: }
961: }
962:
1.8 noro 963: int getchomdeg(V v,P p)
1.1 noro 964: {
965: int m,m1;
966: DCP dc;
967:
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: }
979: }
980:
1.8 noro 981: int getlchomdeg(V v,P p,int *d)
1.1 noro 982: {
983: int m0,m1,d0,d1;
984: DCP dc;
985:
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: }
1004: }
1005:
1.8 noro 1006: int nmonop(P p)
1.1 noro 1007: {
1008: int s;
1009: DCP dc;
1010:
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: }
1024: }
1025:
1.8 noro 1026: int qpcheck(Obj p)
1.1 noro 1027: {
1028: DCP dc;
1029:
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: }
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: {
1050: DCP dc;
1051: P c;
1052:
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);
1.8 noro 1062: if ( !NUM(c) || !uzpcheck((Obj)c) )
1.1 noro 1063: return 0;
1064: }
1065: return 1;
1066: default:
1067: return 0;
1068: }
1069: }
1070:
1.8 noro 1071: int p_mag(P p)
1.1 noro 1072: {
1073: int s;
1074: DCP dc;
1075:
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: }
1085: }
1.2 noro 1086:
1.8 noro 1087: int maxblenp(P p)
1.2 noro 1088: {
1089: int s,t;
1090: DCP dc;
1091:
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) ) {
1.10 noro 1098: t = maxblenp(COEF(dc));
1.2 noro 1099: s = MAX(t,s);
1100: }
1101: return s;
1102: }
1103: }