![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to E.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 E.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/engine/E.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
2: #include "ca.h"
3:
4: void henmv(vl,vn,f,g0,h0,a0,b0,lg,lh,lg0,lh0,q,k,gp,hp)
5: VL vl;
6: VN vn;
7: P f,g0,h0,a0,b0,lg,lh,lg0,lh0;
8: Q q;
9: int k;
10: P *gp,*hp;
11: {
12: P g1,h1,a1,b1;
13: N qn;
14: Q q2;
15:
16: divin((NM(q)),2,&qn); NTOQ(qn,1,q2);
17: adjc(vl,g0,a0,lg0,q,&g1,&a1); adjc(vl,h0,b0,lh0,q,&h1,&b1);
18: henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);
19: }
20:
21: void henmvmain(vl,vn,f,fi0,fi1,gi0,gi1,l0,l1,mod,mod2,k,fr0,fr1)
22: VL vl;
23: VN vn;
24: P f,fi0,fi1,gi0,gi1,l0,l1;
25: Q mod,mod2;
26: int k;
27: P *fr0,*fr1;
28: {
29: V v;
30: int n,i,j;
31: int *md;
32: P x,m,m1,c,q,r,a,s,u,ff,f0,f1;
33: P w0,w1,cf,cfi,t,q1,dvr;
34: P *c0,*c1;
35: P *f0h,*f1h;
36:
37: v = VR(f); n = deg(v,f); MKV(v,x);
38: c0 = (P *)ALLOCA((n+1)*sizeof(P));
39: c1 = (P *)ALLOCA((n+1)*sizeof(P));
40: invl((Q)LC(fi1),mod,(Q *)&c); mulp(vl,fi1,c,&t); cmp(mod,t,&dvr);
41: cm2p(mod,mod2,gi0,&c0[0]); cm2p(mod,mod2,gi1,&c1[0]);
42: for ( i = 1; i <= n; i++ ) {
43: mulp(vl,x,c1[i-1],&m); divsrp(vl,m,dvr,&q,&r); mulp(vl,q,c,&q1);
44: cm2p(mod,mod2,r,&c1[i]);
45: mulp(vl,x,c0[i-1],&m); mulp(vl,q1,fi0,&m1); addp(vl,m,m1,&a);
46: cm2p(mod,mod2,a,&c0[i]);
47: }
48: affine(vl,f,vn,&t); cm2p(mod,mod2,t,&ff);
49: for ( i = 0; vn[i].v; i++ );
50: md = ( int *) ALLOCA((i+1)*sizeof(int));
51: for ( i = 0; vn[i].v; i++ )
52: md[i] = getdeg(vn[i].v,ff);
53: cm2p(mod,mod2,fi0,&f0); affine(vl,l0,vn,&t);
54: if ( NUM(f0) )
55: cm2p(mod,mod2,t,&f0);
56: else
57: cm2p(mod,mod2,t,&COEF(DC(f0)));
58: cm2p(mod,mod2,fi1,&f1); affine(vl,l1,vn,&t);
59: if ( NUM(f1) )
60: cm2p(mod,mod2,t,&f1);
61: else
62: cm2p(mod,mod2,t,&COEF(DC(f1)));
63: W_CALLOC(k,P,f0h); W_CALLOC(k,P,f1h);
64: for ( i = 0; i <= k; i++ ) {
65: exthpc(vl,v,f0,i,&f0h[i]); exthpc(vl,v,f1,i,&f1h[i]);
66: }
67: for ( j = 1; j <= k; j++ ) {
68: for ( i = 0; vn[i].v; i++ )
69: if ( getdeg(vn[i].v,f0)+getdeg(vn[i].v,f1) > md[i] )
70: goto END;
71: for ( i = 0, s = 0; i <= j; i++ ) {
72: mulp(vl,f0h[i],f1h[j-i],&t); addp(vl,s,t,&u); s = u;
73: }
74: cm2p(mod,mod2,s,&t);
75: exthpc(vl,v,ff,j,&u); subm2p(vl,mod,mod2,u,t,&cf);
76: for ( i = 0, w0 = 0, w1 = 0; i <= n; i++ ) {
77: if ( !cf )
78: cfi = 0;
79: else if ( VR(cf) == v )
80: coefp(cf,i,&cfi);
81: else if ( i == 0 )
82: cfi = cf;
83: else
84: cfi = 0;
85: if ( cfi ) {
86: mulp(vl,cfi,c0[i],&m); addp(vl,w0,m,&a); w0 = a;
87: mulp(vl,cfi,c1[i],&m); addp(vl,w1,m,&a); w1 = a;
88: }
89: }
90: cm2p(mod,mod2,w0,&t); addm2p(vl,mod,mod2,f0,t,&a);
91: addm2p(vl,mod,mod2,f0h[j],t,&s); f0h[j] = s; f0 = a;
92: cm2p(mod,mod2,w1,&u); addm2p(vl,mod,mod2,f1,u,&a);
93: addm2p(vl,mod,mod2,f1h[j],u,&s); f1h[j] = s; f1 = a;
94: if ( !t ) {
95: restore(vl,f0,vn,&s); cm2p(mod,mod2,s,&t);
96: if ( divtpz(vl,f,t,&s) ) {
97: *fr0 = t; *fr1 = s;
98: return;
99: }
100: }
101: if ( !u ) {
102: restore(vl,f1,vn,&s); cm2p(mod,mod2,s,&t);
103: if ( divtpz(vl,f,t,&s) ) {
104: *fr0 = s; *fr1 = t;
105: return;
106: }
107: }
108: }
109: END:
110: restore(vl,f0,vn,&t); cm2p(mod,mod2,t,fr0);
111: restore(vl,f1,vn,&t); cm2p(mod,mod2,t,fr1);
112: }
113:
114: /*
115: input : f, fi0, fi1, p, k; f = fi0 * fi1 mod p; ( p, k : integer )
116: output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )
117: */
118:
119: void henzq(f,i0,fi0,i1,fi1,p,k,fr0p,fr1p,gr0p,gr1p,qrp)
120: P f;
121: UM fi0,fi1;
122: int p,k;
123: P i0,i1;
124: P *fr0p,*fr1p,*gr0p,*gr1p;
125: Q *qrp;
126: {
127: N qn;
128: Q q,qq,q2;
129: int n,i;
130: UM wg0,wg1,wf0,wf1;
131: P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
132:
133: n = UDEG(f);
134: wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
135: wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
136: cpyum(fi0,wf0); cpyum(fi1,wf1);
137: eucum(p,wf0,wf1,wg1,wg0);
138: umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
139: umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
140:
141: STOQ(p,q); divin(NM(q),2,&qn); NTOQ(qn,1,q2);
142: for ( i = 1; i < k; i++ ) {
143: #if 0
144: mulp(CO,i1,g0,&m); mulp(CO,i0,g1,&m1); addp(CO,m,m1,&a);
145: if ( NUM(a) ) {
146: for ( STOQ(p,q), j = 1; j < k; j++ ) {
147: mulq(q,q,&qq); q = qq;
148: }
149: f0 = i0; f1 = i1;
150: invl(a,q,&qq);
151: mulp(CO,g0,qq,&s); g0 = s; mulp(CO,g1,qq,&s); g1 = s;
152: break;
153: }
154: #endif
155: /* c = ((f - f0*f1)/q) mod q;
156: q1 = (c*g1) / f1;
157: r1 = (c*g1) mod f1;
158: f1 += (r1 mod q)*q;
159: f0 += ((c*g0 + q1*f0) mod q)*q;
160:
161: d = ((1 - (f1*g0 + f0*g1))/q) mod q;
162: q1 = (d*g0) / f1;
163: r1 = (d*g0) mod f1;
164: g1 += (r1 mod q)*q;
165: g0 += ((c*g0 + q1*f0) mod q)*q;
166: q = q^2;
167: */
168:
169: /* c = ((f - f0*f1)/q) mod q */
170: mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
171: divcp(s,q,&m); cm2p(q,q2,m,&c);
172:
173: /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
174: mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
175: udivpwm(q,s,f1,&q1,&r1);
176:
177: /* f1 = f1 + (r1 mod q)*q; */
178: cm2p(q,q2,r1,&rm);
179: mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
180: f1 = a;
181:
182: /* a1 = (c*g0 + q1*f0) mod q; */
183: mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
184: cm2p(q,q2,a,&a1);
185:
186: /* f0 = f0 + a1*q; */
187: mulpq(a1,(P)q,&a2);
188: addp(CO,f0,a2,&a);
189: f0 = a;
190:
191: /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
192: mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
193: subp(CO,(P)ONE,a,&s);
194: divcp(s,q,&m); cm2p(q,q2,m,&d);
195:
196: /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
197: mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
198:
199: /* g1 = g1 + (r1 mod q )*q; */
200: cm2p(q,q2,r1,&rm);
201: mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
202: g1 = a;
203:
204: /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
205: mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
206: cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
207: addp(CO,g0,a2,&a);
208: g0 = a;
209:
210: /* q = q^2; */
211: mulq(q,q,&qq);
212: q = qq;
213: divin(NM(q),2,&qn); NTOQ(qn,1,q2);
214: }
215: *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
216: }
217:
218: void henzq1(g,h,bound,gcp,hcp,qp)
219: P g,h;
220: Q bound;
221: P *gcp,*hcp;
222: Q *qp;
223: {
224: V v;
225: Q f,q,q1;
226: Q u,t,a,b,s;
227: P c,c1;
228: P tg,th,tr;
229: UM wg,wh,ws,wt,wm;
230: int n,m,modres,mod,index,i;
231: P gc0,hc0;
232: P z,zz,zzz;
233:
234:
235: v = VR(g); n=deg(v,g); m=deg(v,h);
236: norm(g,&a); norm(h,&b);
237: STOQ(m,u); pwrq(a,u,&t);
238: STOQ(n,u); pwrq(b,u,&s);
239: mulq(t,s,&u);
240:
241: factorial(n+m,&t); mulq(u,t,&s);
242: addq(s,s,&f);
243:
244: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
245: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
246: wm = W_UMALLOC(m+n);
247:
248: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
249: mod = lprime[index++];
250: if ( !mod )
251: error("henzq1 : lprime[] exhausted.");
252: if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )
253: continue;
254: ptomp(mod,g,&tg); ptomp(mod,h,&th);
255: srchump(mod,tg,th,&tr);
256: if ( !tr )
257: continue;
258: else
259: modres = CONT((MQ)tr);
260:
261: mptoum(tg,wg); mptoum(th,wh);
262: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
263: DEG(wm) = 0; COEF(wm)[0] = modres;
264: mulum(mod,ws,wm,wt);
265: for ( i = DEG(wt); i >= 0; i-- )
266: if ( ( COEF(wt)[i] * 2 ) > mod )
267: COEF(wt)[i] -= mod;
268: chnrem(mod,v,c,q,wt,&c1,&q1);
269: if ( !ucmpp(c,c1) ) {
270: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
271: if ( NUM(zzz) ) {
272: q = q1; c = c1;
273: break;
274: }
275: }
276: q = q1; c = c1;
277:
278: if ( cmpq(f,q) < 0 )
279: break;
280: }
281: ptozp(c,1,&s,&gc0);
282: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
283: mulp(CO,gc0,g,&z);
284: divsrp(CO,z,h,&zz,&zzz);
285: ptozp(zz,1,&s,(P *)&t);
286: if ( INT((Q)s) )
287: chsgnp(zz,&hc0);
288: else {
289: NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;
290: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
291: }
292: if ( !INT((Q)zzz) ) {
293: NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;
294: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
295: }
296: for ( index = 0; ; ) {
297: mod = lprime[index++];
298: if ( !mod )
299: error("henzq1 : lprime[] exhausted.");
300: if ( !rem(NM((Q)zzz),mod) ||
301: !rem(NM((Q)LC(g)),mod) ||
302: !rem(NM((Q)LC(h)),mod) )
303: continue;
304: for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {
305: mulq(q,q,&q1); q = q1;
306: }
307: *qp = q;
308: invl((Q)zzz,q,&q1);
309: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
310: return;
311: }
312: }
313:
314: void addm2p(vl,mod,mod2,n1,n2,nr)
315: VL vl;
316: Q mod,mod2;
317: P n1,n2,*nr;
318: {
319: P t;
320:
321: addp(vl,n1,n2,&t);
322: if ( !t )
323: *nr = 0;
324: else
325: cm2p(mod,mod2,t,nr);
326: }
327:
328: void subm2p(vl,mod,mod2,n1,n2,nr)
329: VL vl;
330: Q mod,mod2;
331: P n1,n2,*nr;
332: {
333: P t;
334:
335: subp(vl,n1,n2,&t);
336: if ( !t )
337: *nr = 0;
338: else
339: cm2p(mod,mod2,t,nr);
340: }
341:
342: void mulm2p(vl,mod,mod2,n1,n2,nr)
343: VL vl;
344: Q mod,mod2;
345: P n1,n2,*nr;
346: {
347: P t;
348:
349: mulp(vl,n1,n2,&t);
350: if ( !t )
351: *nr = 0;
352: else
353: cm2p(mod,mod2,t,nr);
354: }
355:
356: void cmp(mod,p,pr)
357: Q mod;
358: P p,*pr;
359: {
360: P t;
361: DCP dc,dcr,dcr0;
362:
363: if ( !p )
364: *pr = 0;
365: else if ( NUM(p) )
366: remq((Q)p,mod,(Q *)pr);
367: else {
368: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
369: cmp(mod,COEF(dc),&t);
370: if ( t ) {
371: NEXTDC(dcr0,dcr);
372: DEG(dcr) = DEG(dc);
373: COEF(dcr) = t;
374: }
375: }
376: if ( !dcr0 )
377: *pr = 0;
378: else {
379: NEXT(dcr) = 0;
380: MKP(VR(p),dcr0,*pr);
381: }
382: }
383: }
384:
385: void cm2p(mod,m,p,pr)
386: Q mod,m;
387: P p,*pr;
388: {
389: P t;
390: DCP dc,dcr,dcr0;
391:
392: if ( !p )
393: *pr = 0;
394: else if ( NUM(p) )
395: rem2q((Q)p,mod,m,(Q *)pr);
396: else {
397: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
398: cm2p(mod,m,COEF(dc),&t);
399: if ( t ) {
400: NEXTDC(dcr0,dcr);
401: DEG(dcr) = DEG(dc);
402: COEF(dcr) = t;
403: }
404: }
405: if ( !dcr0 )
406: *pr = 0;
407: else {
408: NEXT(dcr) = 0;
409: MKP(VR(p),dcr0,*pr);
410: }
411: }
412: }
413:
414: void addm2q(mod,mod2,n1,n2,nr)
415: Q mod,mod2;
416: Q n1,n2,*nr;
417: {
418: Q t;
419:
420: addq(n1,n2,&t);
421: if ( !t )
422: *nr = 0;
423: else
424: rem2q(t,mod,mod2,nr);
425: }
426:
427: void subm2q(mod,mod2,n1,n2,nr)
428: Q mod,mod2;
429: Q n1,n2,*nr;
430: {
431: Q t;
432:
433: subq(n1,n2,&t);
434: if ( !t )
435: *nr = 0;
436: else
437: rem2q(t,mod,mod2,nr);
438: }
439:
440: void mulm2q(mod,mod2,n1,n2,nr)
441: Q mod,mod2;
442: Q n1,n2,*nr;
443: {
444: Q t;
445:
446: mulq(n1,n2,&t);
447: if ( !t )
448: *nr = 0;
449: else
450: rem2q(t,mod,mod2,nr);
451: }
452:
453: void rem2q(n,m,m2,nr)
454: Q n,m,m2,*nr;
455: {
456: N q,r,s;
457: int sgn;
458:
459: divn(NM(n),NM(m),&q,&r);
460: if ( !r )
461: *nr = 0;
462: else {
463: sgn = cmpn(r,NM(m2));
464: if ( sgn > 0 ) {
465: subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);
466: } else
467: NTOQ(r,SGN(n),*nr);
468: }
469: }
470:
471: void exthp(vl,p,d,pr)
472: VL vl;
473: P p;
474: int d;
475: P *pr;
476: {
477: P t,t1,a,w,x,xd;
478: DCP dc;
479:
480: if ( d < 0 )
481: *pr = 0;
482: else if ( NUM(p) )
483: if ( d == 0 )
484: *pr = p;
485: else
486: *pr = 0;
487: else {
488: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
489: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
490: pwrp(vl,x,DEG(dc),&xd);
491: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
492: }
493: *pr = w;
494: }
495: }
496:
497: void exthpc(vl,v,p,d,pr)
498: VL vl;
499: V v;
500: P p;
501: int d;
502: P *pr;
503: {
504: P t,t1,a,w,x,xd;
505: DCP dc;
506:
507: if ( v != VR(p) )
508: exthp(vl,p,d,pr);
509: else if ( d < 0 )
510: *pr = 0;
511: else {
512: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
513: exthp(vl,COEF(dc),d,&t);
514: pwrp(vl,x,DEG(dc),&xd);
515: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
516: }
517: *pr = w;
518: }
519: }
520:
521: void cbound(vl,p,b)
522: VL vl;
523: P p;
524: Q *b;
525: {
526: Q n,e,t,m;
527: int k;
528:
529: cmax(p,&n);
530: addq(n,n,&m);
531:
532: k = geldb(vl,p);
533: STOQ(3,t); STOQ(k,e);
534:
535: pwrq(t,e,&n);
536: mulq(m,n,b);
537: }
538:
539: int geldb(vl,p)
540: VL vl;
541: P p;
542: {
543: int m;
544:
545: for ( m = 0; vl; vl = NEXT(vl) )
546: m += getdeg(vl->v,p);
547: return ( m );
548: }
549:
550: int getdeg(v,p)
551: V v;
552: P p;
553: {
554: int m;
555: DCP dc;
556:
557: if ( !p || NUM(p) )
558: return ( 0 );
559: else if ( v == VR(p) )
560: return ( deg(v,p) );
561: else {
562: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) )
563: m = MAX(m,getdeg(v,COEF(dc)));
564: return ( m );
565: }
566: }
567:
568: void cmax(p,b)
569: P p;
570: Q *b;
571: {
572: DCP dc;
573: Q m,m1;
574: N tn;
575:
576: if ( NUM(p) ) {
577: tn = NM((Q)p);
578: NTOQ(tn,1,*b);
579: } else {
580: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
581: cmax(COEF(dc),&m1);
582: if ( cmpq(m1,m) > 0 )
583: m = m1;
584: }
585: *b = m;
586: }
587: }
588:
589: int nextbin(vn,n)
590: VN vn;
591: int n;
592: {
593: int tmp,i,carry;
594:
595: if ( n == 0 )
596: return ( 1 );
597:
598: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
599: tmp = vn[i].n + carry;
600: vn[i].n = tmp % 2;
601: carry = tmp / 2;
602: }
603: return ( carry );
604: }
605:
606: void mulsgn(vn,vnt,n,vn1)
607: VN vn,vnt,vn1;
608: int n;
609: {
610: int i;
611:
612: for ( i = 0; vn[i].v; i++ )
613: vn1[i].n = vn[i].n;
614: for ( i = 0; i < n; i++ )
615: if ( vnt[i].n )
616: vn1[(int)vnt[i].v].n *= -1;
617: }
618:
619: void next(vn)
620: VN vn;
621: {
622: int i,m,n,tmp,carry;
623:
624: for ( m = 0, i = 0; vn[i].v; i++ )
625: m = MAX(m,ABS(vn[i].n));
626: if ( m == 0 ) {
627: vn[--i].n = 1;
628: return;
629: }
630: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
631: tmp = vn[i].n + carry;
632: vn[i].n = tmp % m;
633: carry = tmp / m;
634: }
635: if ( ( i == -1 ) && carry ) {
636: for ( i = 0; vn[i].v; i++ )
637: vn[i].n = 0;
638: vn[--i].n = m;
639: } else {
640: for ( n = 0, i = 0; vn[i].v; i++ )
641: n = MAX(n,ABS(vn[i].n));
642: if ( n < m - 1 )
643: vn[--i].n = m - 1;
644: }
645: }
646:
647: void clctv(vl,p,nvlp)
648: VL vl;
649: P p;
650: VL *nvlp;
651: {
652: int i,n;
653: VL tvl;
654: VN tvn;
655:
656: if ( !p || NUM(p) ) {
657: *nvlp = 0;
658: return;
659: }
660:
661: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
662: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
663: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
664: tvn[i].v = tvl->v;
665: tvn[i].n = 0;
666: }
667:
668: markv(tvn,n,p);
669: vntovl(tvn,n,nvlp);
670: }
671:
672: void markv(vn,n,p)
673: VN vn;
674: int n;
675: P p;
676: {
677: V v;
678: DCP dc;
679: int i;
680:
681: if ( NUM(p) )
682: return;
683: v = VR(p);
684: for ( i = 0, v = VR(p); i < n; i++ )
685: if ( v == vn[i].v ) {
686: vn[i].n = 1;
687: break;
688: }
689: for ( dc = DC(p); dc; dc = NEXT(dc) )
690: markv(vn,n,COEF(dc));
691: }
692:
693: void vntovl(vn,n,vlp)
694: VN vn;
695: int n;
696: VL *vlp;
697: {
698: int i;
699: VL tvl,tvl0;
700:
701: for ( i = 0, tvl0 = 0; ; ) {
702: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
703: if ( i == n )
704: break;
705: else {
706: if ( !tvl0 ) {
707: NEWVL(tvl0);
708: tvl = tvl0;
709: } else {
710: NEWVL(NEXT(tvl));
711: tvl = NEXT(tvl);
712: }
713: tvl->v = vn[i++].v;
714: }
715: }
716: if ( tvl0 )
717: NEXT(tvl) = 0;
718: *vlp = tvl0;
719: }
720:
721: int dbound(v,f)
722: V v;
723: P f;
724: {
725: int m;
726: DCP dc;
727:
728: if ( !f )
729: return ( -1 );
730: else if ( v != VR(f) )
731: return homdeg(f);
732: else {
733: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )
734: m = MAX(m,homdeg(COEF(dc)));
735: return ( m );
736: }
737: }
738:
739: int homdeg(f)
740: P f;
741: {
742: int m;
743: DCP dc;
744:
745: if ( !f )
746: return ( -1 );
747: else if ( NUM(f) )
748: return( 0 );
749: else {
750: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) )
751: m = MAX(m,QTOS(DEG(dc))+homdeg(COEF(dc)));
752: return ( m );
753: }
754: }
755:
756: int minhomdeg(f)
757: P f;
758: {
759: int m;
760: DCP dc;
761:
762: if ( !f )
763: return ( -1 );
764: else if ( NUM(f) )
765: return( 0 );
766: else {
767: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) )
768: m = MIN(m,QTOS(DEG(dc))+minhomdeg(COEF(dc)));
769: return ( m );
770: }
771: }
772:
773: void adjc(vl,f,a,lc0,q,fr,ar)
774: VL vl;
775: P f,a,lc0;
776: Q q;
777: P *fr,*ar;
778: {
779: P m,m1;
780: Q t;
781:
782: invl((Q)LC(f),q,&t);
783: mulq((Q)lc0,t,(Q *)&m);
784: mulpq(f,m,&m1); cmp(q,m1,fr);
785: invl((Q)m,q,&t);
786: mulpq(a,(P)t,&m1);
787: cmp(q,m1,ar);
788: }
789:
790: #if 1
791: void affinemain(vl,p,v0,n,pl,pr)
792: VL vl;
793: V v0;
794: int n;
795: P *pl;
796: P p;
797: P *pr;
798: {
799: P x,t,m,c,s,a;
800: DCP dc;
801: Q d;
802:
803: if ( !p )
804: *pr = 0;
805: else if ( NUM(p) )
806: *pr = p;
807: else if ( VR(p) != v0 ) {
808: MKV(VR(p),x);
809: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
810: affinemain(vl,COEF(dc),v0,n,pl,&t);
811: if ( DEG(dc) ) {
812: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
813: addp(vl,m,c,&a); c = a;
814: } else {
815: addp(vl,t,c,&a); c = a;
816: }
817: }
818: *pr = c;
819: } else {
820: dc = DC(p);
821: c = COEF(dc);
822: for ( d = DEG(dc), dc = NEXT(dc);
823: dc; d = DEG(dc), dc = NEXT(dc) ) {
824: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
825: addp(vl,m,COEF(dc),&c);
826: }
827: if ( d ) {
828: mulp(vl,pl[QTOS(d)],c,&m); c = m;
829: }
830: *pr = c;
831: }
832: }
833: #endif
834:
835: #if 0
836: affine(vl,p,vn,r)
837: VL vl;
838: P p;
839: VN vn;
840: P *r;
841: {
842: int n,d,d1,i;
843: Q *q;
844: Q **bc;
845:
846: if ( !p || NUM(p) )
847: *r = p;
848: else {
849: for ( i = 0, d = 0; vn[i].v; i++ )
850: d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
851: W_CALLOC(d+1,Q *,bc);
852: for ( i = 0; i <= d; i++ )
853: W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
854: afmain(vl,bc,p,vn,r);
855: }
856: }
857:
858: afmain(vl,bc,p,vn,r)
859: VL vl;
860: Q **bc;
861: P p;
862: VN vn;
863: P *r;
864: {
865: P t,s,u;
866: P *c,*rc;
867: Q *q;
868: DCP dc;
869: int n,i,j;
870:
871: if ( !p || NUM(p) || !vn[0].v )
872: *r = p;
873: else if ( vn[0].v != VR(p) ) {
874: for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
875: if ( vn[i].v )
876: afmain(vl,bc,p,vn+i,r);
877: else {
878: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
879: for ( dc = DC(p); dc; dc = NEXT(dc) )
880: afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
881: plisttop(c,VR(p),n,r);
882: }
883: } else {
884: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
885: W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
886: for ( dc = DC(p); dc; dc = NEXT(dc) )
887: afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
888: if ( !vn[0].n )
889: bcopy(c,rc,sizeof(P)*(n+1));
890: else {
891: for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
892: mulq(q[i-1],q[1],&q[i]);
893: for ( j = 0; j <= n; rc[j] = t, j++ )
894: for ( i = j, t = 0; i <= n; i++ )
895: if ( c[i] )
896: mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
897: addp(CO,u,t,&s), t = s;
898: }
899: plisttop(rc,VR(p),n,r);
900: }
901: }
902: #endif
903:
904: void restore(vl,f,vn,fr)
905: VL vl;
906: P f;
907: VN vn;
908: P *fr;
909: {
910: int i;
911: P vv,g,g1,t;
912: Q s;
913:
914: g = f;
915: for ( i = 0; vn[i].v; i++ ) {
916: MKV(vn[i].v,t);
917: if ( vn[i].n ) {
918: STOQ(-vn[i].n,s);
919: addp(vl,t,(P)s,&vv);
920: } else
921: vv = t;
922:
923: substp(vl,g,vn[i].v,vv,&g1); g = g1;
924: }
925: *fr = g;
926: }
927:
928: void mergev(vl,vl1,vl2,nvlp)
929: VL vl,vl1,vl2,*nvlp;
930: {
931: int i,n;
932: VL tvl;
933: VN vn;
934:
935: if ( !vl1 ) {
936: *nvlp = vl2; return;
937: } else if ( !vl2 ) {
938: *nvlp = vl1; return;
939: }
940: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
941: n = i;
942: W_CALLOC(n,struct oVN,vn);
943: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
944: vn[i].v = tvl->v;
945: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
946: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
947: i++;
948: if ( i == n )
949: break;
950: else
951: vn[i].n = 1;
952: }
953: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
954: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
955: i++;
956: if ( i == n )
957: break;
958: else
959: vn[i].n = 1;
960: }
961: vntovl(vn,n,nvlp);
962: }
963:
964: #if 0
965: void substvp(vl,f,vn,g)
966: VL vl;
967: P f;
968: VN vn;
969: P *g;
970: {
971: V v;
972: int i;
973: P h,h1;
974: Q t;
975:
976: h = f;
977: for ( i = 0; v = vn[i].v; i++ ) {
978: STOQ(vn[i].n,t);
979: substp(vl,h,v,(P)t,&h1); h = h1;
980: }
981: *g = h;
982: }
983:
984: void affine(vl,f,vn,fr)
985: VL vl;
986: P f;
987: VN vn;
988: P *fr;
989: {
990: int i,j,n;
991: P vv,g,g1,t,u;
992: Q s;
993: int *dlist;
994: P **plist;
995:
996: for ( n = 0; vn[n].v; n++);
997: dlist = (int *)ALLOCA((n+1)*sizeof(int));
998: plist = (P **)ALLOCA((n+1)*sizeof(P *));
999: for ( i = 0; vn[i].v; i++ ) {
1000: if ( !vn[i].n )
1001: continue;
1002: dlist[i] = getdeg(vn[i].v,f);
1003: plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
1004:
1005: MKV(vn[i].v,t);
1006: if ( vn[i].n ) {
1007: STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
1008: } else
1009: vv = t;
1010:
1011: for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
1012: plist[i][j] = t;
1013: mulp(vl,t,vv,&u);
1014: t = u;
1015: }
1016: plist[i][j] = t;
1017: }
1018:
1019: g = f;
1020: for ( i = 0; vn[i].v; i++ ) {
1021: if ( !vn[i].n )
1022: continue;
1023: affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
1024: }
1025: *fr = g;
1026: }
1027: #endif