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