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