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