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