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