![[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.5 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.4 2001/04/20 02:27:52 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);
516: }
517: }
518:
519: void exthp(vl,p,d,pr)
520: VL vl;
521: P p;
522: int d;
523: P *pr;
524: {
525: P t,t1,a,w,x,xd;
526: DCP dc;
527:
528: if ( d < 0 )
529: *pr = 0;
530: else if ( NUM(p) )
531: if ( d == 0 )
532: *pr = p;
533: else
534: *pr = 0;
535: else {
536: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
537: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
538: pwrp(vl,x,DEG(dc),&xd);
539: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
540: }
541: *pr = w;
542: }
543: }
544:
545: void exthpc(vl,v,p,d,pr)
546: VL vl;
547: V v;
548: P p;
549: int d;
550: P *pr;
551: {
552: P t,t1,a,w,x,xd;
553: DCP dc;
554:
555: if ( v != VR(p) )
556: exthp(vl,p,d,pr);
557: else if ( d < 0 )
558: *pr = 0;
559: else {
560: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
561: exthp(vl,COEF(dc),d,&t);
562: pwrp(vl,x,DEG(dc),&xd);
563: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
564: }
565: *pr = w;
566: }
567: }
568:
569: void cbound(vl,p,b)
570: VL vl;
571: P p;
572: Q *b;
573: {
574: Q n,e,t,m;
575: int k;
576:
577: cmax(p,&n);
578: addq(n,n,&m);
579:
580: k = geldb(vl,p);
581: STOQ(3,t); STOQ(k,e);
582:
583: pwrq(t,e,&n);
584: mulq(m,n,b);
585: }
586:
587: int geldb(vl,p)
588: VL vl;
589: P p;
590: {
591: int m;
592:
593: for ( m = 0; vl; vl = NEXT(vl) )
594: m += getdeg(vl->v,p);
595: return ( m );
596: }
597:
598: int getdeg(v,p)
599: V v;
600: P p;
601: {
1.4 noro 602: int m,t;
1.1 noro 603: DCP dc;
604:
605: if ( !p || NUM(p) )
606: return ( 0 );
607: else if ( v == VR(p) )
608: return ( deg(v,p) );
609: else {
1.4 noro 610: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
611: t = getdeg(v,COEF(dc));
612: m = MAX(m,t);
613: }
1.1 noro 614: return ( m );
615: }
616: }
617:
618: void cmax(p,b)
619: P p;
620: Q *b;
621: {
622: DCP dc;
623: Q m,m1;
624: N tn;
625:
626: if ( NUM(p) ) {
627: tn = NM((Q)p);
628: NTOQ(tn,1,*b);
629: } else {
630: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
631: cmax(COEF(dc),&m1);
632: if ( cmpq(m1,m) > 0 )
633: m = m1;
634: }
635: *b = m;
636: }
637: }
638:
639: int nextbin(vn,n)
640: VN vn;
641: int n;
642: {
643: int tmp,i,carry;
644:
645: if ( n == 0 )
646: return ( 1 );
647:
648: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
649: tmp = vn[i].n + carry;
650: vn[i].n = tmp % 2;
651: carry = tmp / 2;
652: }
653: return ( carry );
654: }
655:
656: void mulsgn(vn,vnt,n,vn1)
657: VN vn,vnt,vn1;
658: int n;
659: {
660: int i;
661:
662: for ( i = 0; vn[i].v; i++ )
663: vn1[i].n = vn[i].n;
664: for ( i = 0; i < n; i++ )
665: if ( vnt[i].n )
666: vn1[(int)vnt[i].v].n *= -1;
667: }
668:
669: void next(vn)
670: VN vn;
671: {
672: int i,m,n,tmp,carry;
673:
674: for ( m = 0, i = 0; vn[i].v; i++ )
675: m = MAX(m,ABS(vn[i].n));
676: if ( m == 0 ) {
677: vn[--i].n = 1;
678: return;
679: }
680: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
681: tmp = vn[i].n + carry;
682: vn[i].n = tmp % m;
683: carry = tmp / m;
684: }
685: if ( ( i == -1 ) && carry ) {
686: for ( i = 0; vn[i].v; i++ )
687: vn[i].n = 0;
688: vn[--i].n = m;
689: } else {
690: for ( n = 0, i = 0; vn[i].v; i++ )
691: n = MAX(n,ABS(vn[i].n));
692: if ( n < m - 1 )
693: vn[--i].n = m - 1;
694: }
695: }
696:
697: void clctv(vl,p,nvlp)
698: VL vl;
699: P p;
700: VL *nvlp;
701: {
702: int i,n;
703: VL tvl;
704: VN tvn;
705:
706: if ( !p || NUM(p) ) {
707: *nvlp = 0;
708: return;
709: }
710:
711: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
712: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
713: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
714: tvn[i].v = tvl->v;
715: tvn[i].n = 0;
716: }
717:
718: markv(tvn,n,p);
719: vntovl(tvn,n,nvlp);
720: }
721:
722: void markv(vn,n,p)
723: VN vn;
724: int n;
725: P p;
726: {
727: V v;
728: DCP dc;
729: int i;
730:
731: if ( NUM(p) )
732: return;
733: v = VR(p);
734: for ( i = 0, v = VR(p); i < n; i++ )
735: if ( v == vn[i].v ) {
736: vn[i].n = 1;
737: break;
738: }
739: for ( dc = DC(p); dc; dc = NEXT(dc) )
740: markv(vn,n,COEF(dc));
741: }
742:
743: void vntovl(vn,n,vlp)
744: VN vn;
745: int n;
746: VL *vlp;
747: {
748: int i;
749: VL tvl,tvl0;
750:
751: for ( i = 0, tvl0 = 0; ; ) {
752: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
753: if ( i == n )
754: break;
755: else {
756: if ( !tvl0 ) {
757: NEWVL(tvl0);
758: tvl = tvl0;
759: } else {
760: NEWVL(NEXT(tvl));
761: tvl = NEXT(tvl);
762: }
763: tvl->v = vn[i++].v;
764: }
765: }
766: if ( tvl0 )
767: NEXT(tvl) = 0;
768: *vlp = tvl0;
769: }
770:
771: int dbound(v,f)
772: V v;
773: P f;
774: {
1.5 ! noro 775: int m,t;
1.1 noro 776: DCP dc;
777:
778: if ( !f )
779: return ( -1 );
780: else if ( v != VR(f) )
781: return homdeg(f);
782: else {
1.5 ! noro 783: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 784: t = homdeg(COEF(dc));
! 785: m = MAX(m,t);
! 786: }
1.1 noro 787: return ( m );
788: }
789: }
790:
791: int homdeg(f)
792: P f;
793: {
1.5 ! noro 794: int m,t;
1.1 noro 795: DCP dc;
796:
797: if ( !f )
798: return ( -1 );
799: else if ( NUM(f) )
800: return( 0 );
801: else {
1.5 ! noro 802: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
! 803: t = QTOS(DEG(dc))+homdeg(COEF(dc));
! 804: m = MAX(m,t);
! 805: }
1.1 noro 806: return ( m );
807: }
808: }
809:
810: int minhomdeg(f)
811: P f;
812: {
1.5 ! noro 813: int m,t;
1.1 noro 814: DCP dc;
815:
816: if ( !f )
817: return ( -1 );
818: else if ( NUM(f) )
819: return( 0 );
820: else {
1.5 ! noro 821: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
! 822: t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
! 823: m = MIN(m,t);
! 824: }
1.1 noro 825: return ( m );
826: }
827: }
828:
829: void adjc(vl,f,a,lc0,q,fr,ar)
830: VL vl;
831: P f,a,lc0;
832: Q q;
833: P *fr,*ar;
834: {
835: P m,m1;
836: Q t;
837:
838: invl((Q)LC(f),q,&t);
839: mulq((Q)lc0,t,(Q *)&m);
840: mulpq(f,m,&m1); cmp(q,m1,fr);
841: invl((Q)m,q,&t);
842: mulpq(a,(P)t,&m1);
843: cmp(q,m1,ar);
844: }
845:
846: #if 1
847: void affinemain(vl,p,v0,n,pl,pr)
848: VL vl;
849: V v0;
850: int n;
851: P *pl;
852: P p;
853: P *pr;
854: {
855: P x,t,m,c,s,a;
856: DCP dc;
857: Q d;
858:
859: if ( !p )
860: *pr = 0;
861: else if ( NUM(p) )
862: *pr = p;
863: else if ( VR(p) != v0 ) {
864: MKV(VR(p),x);
865: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
866: affinemain(vl,COEF(dc),v0,n,pl,&t);
867: if ( DEG(dc) ) {
868: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
869: addp(vl,m,c,&a); c = a;
870: } else {
871: addp(vl,t,c,&a); c = a;
872: }
873: }
874: *pr = c;
875: } else {
876: dc = DC(p);
877: c = COEF(dc);
878: for ( d = DEG(dc), dc = NEXT(dc);
879: dc; d = DEG(dc), dc = NEXT(dc) ) {
880: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
881: addp(vl,m,COEF(dc),&c);
882: }
883: if ( d ) {
884: mulp(vl,pl[QTOS(d)],c,&m); c = m;
885: }
886: *pr = c;
887: }
888: }
889: #endif
890:
891: #if 0
892: affine(vl,p,vn,r)
893: VL vl;
894: P p;
895: VN vn;
896: P *r;
897: {
898: int n,d,d1,i;
899: Q *q;
900: Q **bc;
901:
902: if ( !p || NUM(p) )
903: *r = p;
904: else {
905: for ( i = 0, d = 0; vn[i].v; i++ )
906: d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
907: W_CALLOC(d+1,Q *,bc);
908: for ( i = 0; i <= d; i++ )
909: W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
910: afmain(vl,bc,p,vn,r);
911: }
912: }
913:
914: afmain(vl,bc,p,vn,r)
915: VL vl;
916: Q **bc;
917: P p;
918: VN vn;
919: P *r;
920: {
921: P t,s,u;
922: P *c,*rc;
923: Q *q;
924: DCP dc;
925: int n,i,j;
926:
927: if ( !p || NUM(p) || !vn[0].v )
928: *r = p;
929: else if ( vn[0].v != VR(p) ) {
930: for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
931: if ( vn[i].v )
932: afmain(vl,bc,p,vn+i,r);
933: else {
934: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
935: for ( dc = DC(p); dc; dc = NEXT(dc) )
936: afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
937: plisttop(c,VR(p),n,r);
938: }
939: } else {
940: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
941: W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
942: for ( dc = DC(p); dc; dc = NEXT(dc) )
943: afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
944: if ( !vn[0].n )
945: bcopy(c,rc,sizeof(P)*(n+1));
946: else {
947: for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
948: mulq(q[i-1],q[1],&q[i]);
949: for ( j = 0; j <= n; rc[j] = t, j++ )
950: for ( i = j, t = 0; i <= n; i++ )
951: if ( c[i] )
952: mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
953: addp(CO,u,t,&s), t = s;
954: }
955: plisttop(rc,VR(p),n,r);
956: }
957: }
958: #endif
959:
960: void restore(vl,f,vn,fr)
961: VL vl;
962: P f;
963: VN vn;
964: P *fr;
965: {
966: int i;
967: P vv,g,g1,t;
968: Q s;
969:
970: g = f;
971: for ( i = 0; vn[i].v; i++ ) {
972: MKV(vn[i].v,t);
973: if ( vn[i].n ) {
974: STOQ(-vn[i].n,s);
975: addp(vl,t,(P)s,&vv);
976: } else
977: vv = t;
978:
979: substp(vl,g,vn[i].v,vv,&g1); g = g1;
980: }
981: *fr = g;
982: }
983:
984: void mergev(vl,vl1,vl2,nvlp)
985: VL vl,vl1,vl2,*nvlp;
986: {
987: int i,n;
988: VL tvl;
989: VN vn;
990:
991: if ( !vl1 ) {
992: *nvlp = vl2; return;
993: } else if ( !vl2 ) {
994: *nvlp = vl1; return;
995: }
996: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
997: n = i;
998: W_CALLOC(n,struct oVN,vn);
999: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
1000: vn[i].v = tvl->v;
1001: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
1002: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
1003: i++;
1004: if ( i == n )
1005: break;
1006: else
1007: vn[i].n = 1;
1008: }
1009: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
1010: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
1011: i++;
1012: if ( i == n )
1013: break;
1014: else
1015: vn[i].n = 1;
1016: }
1017: vntovl(vn,n,nvlp);
1018: }
1019:
1020: #if 0
1021: void substvp(vl,f,vn,g)
1022: VL vl;
1023: P f;
1024: VN vn;
1025: P *g;
1026: {
1027: V v;
1028: int i;
1029: P h,h1;
1030: Q t;
1031:
1032: h = f;
1033: for ( i = 0; v = vn[i].v; i++ ) {
1034: STOQ(vn[i].n,t);
1035: substp(vl,h,v,(P)t,&h1); h = h1;
1036: }
1037: *g = h;
1038: }
1039:
1040: void affine(vl,f,vn,fr)
1041: VL vl;
1042: P f;
1043: VN vn;
1044: P *fr;
1045: {
1046: int i,j,n;
1047: P vv,g,g1,t,u;
1048: Q s;
1049: int *dlist;
1050: P **plist;
1051:
1052: for ( n = 0; vn[n].v; n++);
1053: dlist = (int *)ALLOCA((n+1)*sizeof(int));
1054: plist = (P **)ALLOCA((n+1)*sizeof(P *));
1055: for ( i = 0; vn[i].v; i++ ) {
1056: if ( !vn[i].n )
1057: continue;
1058: dlist[i] = getdeg(vn[i].v,f);
1059: plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
1060:
1061: MKV(vn[i].v,t);
1062: if ( vn[i].n ) {
1063: STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
1064: } else
1065: vv = t;
1066:
1067: for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
1068: plist[i][j] = t;
1069: mulp(vl,t,vv,&u);
1070: t = u;
1071: }
1072: plist[i][j] = t;
1073: }
1074:
1075: g = f;
1076: for ( i = 0; vn[i].v; i++ ) {
1077: if ( !vn[i].n )
1078: continue;
1079: affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
1080: }
1081: *fr = g;
1082: }
1083: #endif