![[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.7 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/E.c,v 1.6 2001/05/02 09:03:53 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; ; ) {
1.7 ! noro 297: mod = get_lprime(index++);
1.1 noro 298: if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) )
299: continue;
300: ptomp(mod,g,&tg); ptomp(mod,h,&th);
301: srchump(mod,tg,th,&tr);
302: if ( !tr )
303: continue;
304: else
305: modres = CONT((MQ)tr);
306:
307: mptoum(tg,wg); mptoum(th,wh);
308: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
309: DEG(wm) = 0; COEF(wm)[0] = modres;
310: mulum(mod,ws,wm,wt);
311: for ( i = DEG(wt); i >= 0; i-- )
312: if ( ( COEF(wt)[i] * 2 ) > mod )
313: COEF(wt)[i] -= mod;
314: chnrem(mod,v,c,q,wt,&c1,&q1);
315: if ( !ucmpp(c,c1) ) {
316: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
317: if ( NUM(zzz) ) {
318: q = q1; c = c1;
319: break;
320: }
321: }
322: q = q1; c = c1;
323:
324: if ( cmpq(f,q) < 0 )
325: break;
326: }
327: ptozp(c,1,&s,&gc0);
328: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
329: mulp(CO,gc0,g,&z);
330: divsrp(CO,z,h,&zz,&zzz);
331: ptozp(zz,1,&s,(P *)&t);
332: if ( INT((Q)s) )
333: chsgnp(zz,&hc0);
334: else {
335: NTOQ(DN((Q)s),1,q); mulq(q,(Q)zzz,&q1); zzz = (P)q1;
336: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
337: }
338: if ( !INT((Q)zzz) ) {
339: NTOQ(DN((Q)zzz),1,q); NTOQ(NM((Q)zzz),1,q1); zzz = (P)q1;
340: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
341: }
342: for ( index = 0; ; ) {
1.7 ! noro 343: mod = get_lprime(index++);
1.1 noro 344: if ( !rem(NM((Q)zzz),mod) ||
345: !rem(NM((Q)LC(g)),mod) ||
346: !rem(NM((Q)LC(h)),mod) )
347: continue;
348: for ( STOQ(mod,q); cmpq(q,bound) < 0; ) {
349: mulq(q,q,&q1); q = q1;
350: }
351: *qp = q;
352: invl((Q)zzz,q,&q1);
353: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
354: return;
355: }
356: }
357:
358: void addm2p(vl,mod,mod2,n1,n2,nr)
359: VL vl;
360: Q mod,mod2;
361: P n1,n2,*nr;
362: {
363: P t;
364:
365: addp(vl,n1,n2,&t);
366: if ( !t )
367: *nr = 0;
368: else
369: cm2p(mod,mod2,t,nr);
370: }
371:
372: void subm2p(vl,mod,mod2,n1,n2,nr)
373: VL vl;
374: Q mod,mod2;
375: P n1,n2,*nr;
376: {
377: P t;
378:
379: subp(vl,n1,n2,&t);
380: if ( !t )
381: *nr = 0;
382: else
383: cm2p(mod,mod2,t,nr);
384: }
385:
386: void mulm2p(vl,mod,mod2,n1,n2,nr)
387: VL vl;
388: Q mod,mod2;
389: P n1,n2,*nr;
390: {
391: P t;
392:
393: mulp(vl,n1,n2,&t);
394: if ( !t )
395: *nr = 0;
396: else
397: cm2p(mod,mod2,t,nr);
398: }
399:
400: void cmp(mod,p,pr)
401: Q mod;
402: P p,*pr;
403: {
404: P t;
405: DCP dc,dcr,dcr0;
406:
407: if ( !p )
408: *pr = 0;
409: else if ( NUM(p) )
410: remq((Q)p,mod,(Q *)pr);
411: else {
412: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
413: cmp(mod,COEF(dc),&t);
414: if ( t ) {
415: NEXTDC(dcr0,dcr);
416: DEG(dcr) = DEG(dc);
417: COEF(dcr) = t;
418: }
419: }
420: if ( !dcr0 )
421: *pr = 0;
422: else {
423: NEXT(dcr) = 0;
424: MKP(VR(p),dcr0,*pr);
425: }
426: }
427: }
428:
429: void cm2p(mod,m,p,pr)
430: Q mod,m;
431: P p,*pr;
432: {
433: P t;
434: DCP dc,dcr,dcr0;
435:
436: if ( !p )
437: *pr = 0;
438: else if ( NUM(p) )
439: rem2q((Q)p,mod,m,(Q *)pr);
440: else {
441: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
442: cm2p(mod,m,COEF(dc),&t);
443: if ( t ) {
444: NEXTDC(dcr0,dcr);
445: DEG(dcr) = DEG(dc);
446: COEF(dcr) = t;
447: }
448: }
449: if ( !dcr0 )
450: *pr = 0;
451: else {
452: NEXT(dcr) = 0;
453: MKP(VR(p),dcr0,*pr);
454: }
455: }
456: }
457:
458: void addm2q(mod,mod2,n1,n2,nr)
459: Q mod,mod2;
460: Q n1,n2,*nr;
461: {
462: Q t;
463:
464: addq(n1,n2,&t);
465: if ( !t )
466: *nr = 0;
467: else
468: rem2q(t,mod,mod2,nr);
469: }
470:
471: void subm2q(mod,mod2,n1,n2,nr)
472: Q mod,mod2;
473: Q n1,n2,*nr;
474: {
475: Q t;
476:
477: subq(n1,n2,&t);
478: if ( !t )
479: *nr = 0;
480: else
481: rem2q(t,mod,mod2,nr);
482: }
483:
484: void mulm2q(mod,mod2,n1,n2,nr)
485: Q mod,mod2;
486: Q n1,n2,*nr;
487: {
488: Q t;
489:
490: mulq(n1,n2,&t);
491: if ( !t )
492: *nr = 0;
493: else
494: rem2q(t,mod,mod2,nr);
495: }
496:
497: void rem2q(n,m,m2,nr)
498: Q n,m,m2,*nr;
499: {
500: N q,r,s;
501: int sgn;
502:
503: divn(NM(n),NM(m),&q,&r);
504: if ( !r )
505: *nr = 0;
506: else {
507: sgn = cmpn(r,NM(m2));
508: if ( sgn > 0 ) {
509: subn(NM(m),r,&s); NTOQ(s,-1*SGN(n),*nr);
510: } else
511: NTOQ(r,SGN(n),*nr);
1.6 noro 512: }
513: }
514:
515: /*
516: extract d-homogeneous part with respect to vl - {v}
517: */
518:
519: void exthpc_generic(vl,p,d,v,pr)
520: VL vl;
521: P p;
522: int d;
523: V v;
524: P *pr;
525: {
526: P w,x,t,t1,a,xd;
527: V v0;
528: DCP dc;
529:
530: if ( d < 0 || !p )
531: *pr = 0;
532: else if ( NUM(p) )
533: if ( d == 0 )
534: *pr = p;
535: else
536: *pr = 0;
537: else if ( v == VR(p) )
538: exthpc(vl,v,p,d,pr);
539: else {
540: v0 = VR(p);
541: for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
542: exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t);
543: pwrp(vl,x,DEG(dc),&xd);
544: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
545: }
546: *pr = w;
1.1 noro 547: }
548: }
549:
550: void exthp(vl,p,d,pr)
551: VL vl;
552: P p;
553: int d;
554: P *pr;
555: {
556: P t,t1,a,w,x,xd;
557: DCP dc;
558:
559: if ( d < 0 )
560: *pr = 0;
561: else if ( NUM(p) )
562: if ( d == 0 )
563: *pr = p;
564: else
565: *pr = 0;
566: else {
567: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
568: exthp(vl,COEF(dc),d - QTOS(DEG(dc)),&t);
569: pwrp(vl,x,DEG(dc),&xd);
570: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
571: }
572: *pr = w;
573: }
574: }
575:
576: void exthpc(vl,v,p,d,pr)
577: VL vl;
578: V v;
579: P p;
580: int d;
581: P *pr;
582: {
583: P t,t1,a,w,x,xd;
584: DCP dc;
585:
586: if ( v != VR(p) )
587: exthp(vl,p,d,pr);
588: else if ( d < 0 )
589: *pr = 0;
590: else {
591: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
592: exthp(vl,COEF(dc),d,&t);
593: pwrp(vl,x,DEG(dc),&xd);
594: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
595: }
596: *pr = w;
597: }
598: }
599:
600: void cbound(vl,p,b)
601: VL vl;
602: P p;
603: Q *b;
604: {
605: Q n,e,t,m;
606: int k;
607:
608: cmax(p,&n);
609: addq(n,n,&m);
610:
611: k = geldb(vl,p);
612: STOQ(3,t); STOQ(k,e);
613:
614: pwrq(t,e,&n);
615: mulq(m,n,b);
616: }
617:
618: int geldb(vl,p)
619: VL vl;
620: P p;
621: {
622: int m;
623:
624: for ( m = 0; vl; vl = NEXT(vl) )
625: m += getdeg(vl->v,p);
626: return ( m );
627: }
628:
629: int getdeg(v,p)
630: V v;
631: P p;
632: {
1.4 noro 633: int m,t;
1.1 noro 634: DCP dc;
635:
636: if ( !p || NUM(p) )
637: return ( 0 );
638: else if ( v == VR(p) )
639: return ( deg(v,p) );
640: else {
1.4 noro 641: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
642: t = getdeg(v,COEF(dc));
643: m = MAX(m,t);
644: }
1.1 noro 645: return ( m );
646: }
647: }
648:
649: void cmax(p,b)
650: P p;
651: Q *b;
652: {
653: DCP dc;
654: Q m,m1;
655: N tn;
656:
657: if ( NUM(p) ) {
658: tn = NM((Q)p);
659: NTOQ(tn,1,*b);
660: } else {
661: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
662: cmax(COEF(dc),&m1);
663: if ( cmpq(m1,m) > 0 )
664: m = m1;
665: }
666: *b = m;
667: }
668: }
669:
670: int nextbin(vn,n)
671: VN vn;
672: int n;
673: {
674: int tmp,i,carry;
675:
676: if ( n == 0 )
677: return ( 1 );
678:
679: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
680: tmp = vn[i].n + carry;
681: vn[i].n = tmp % 2;
682: carry = tmp / 2;
683: }
684: return ( carry );
685: }
686:
687: void mulsgn(vn,vnt,n,vn1)
688: VN vn,vnt,vn1;
689: int n;
690: {
691: int i;
692:
693: for ( i = 0; vn[i].v; i++ )
694: vn1[i].n = vn[i].n;
695: for ( i = 0; i < n; i++ )
696: if ( vnt[i].n )
697: vn1[(int)vnt[i].v].n *= -1;
698: }
699:
700: void next(vn)
701: VN vn;
702: {
703: int i,m,n,tmp,carry;
704:
705: for ( m = 0, i = 0; vn[i].v; i++ )
706: m = MAX(m,ABS(vn[i].n));
707: if ( m == 0 ) {
708: vn[--i].n = 1;
709: return;
710: }
711: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
712: tmp = vn[i].n + carry;
713: vn[i].n = tmp % m;
714: carry = tmp / m;
715: }
716: if ( ( i == -1 ) && carry ) {
717: for ( i = 0; vn[i].v; i++ )
718: vn[i].n = 0;
719: vn[--i].n = m;
720: } else {
721: for ( n = 0, i = 0; vn[i].v; i++ )
722: n = MAX(n,ABS(vn[i].n));
723: if ( n < m - 1 )
724: vn[--i].n = m - 1;
725: }
726: }
727:
728: void clctv(vl,p,nvlp)
729: VL vl;
730: P p;
731: VL *nvlp;
732: {
733: int i,n;
734: VL tvl;
735: VN tvn;
736:
737: if ( !p || NUM(p) ) {
738: *nvlp = 0;
739: return;
740: }
741:
742: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
743: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
744: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
745: tvn[i].v = tvl->v;
746: tvn[i].n = 0;
747: }
748:
749: markv(tvn,n,p);
750: vntovl(tvn,n,nvlp);
751: }
752:
753: void markv(vn,n,p)
754: VN vn;
755: int n;
756: P p;
757: {
758: V v;
759: DCP dc;
760: int i;
761:
762: if ( NUM(p) )
763: return;
764: v = VR(p);
765: for ( i = 0, v = VR(p); i < n; i++ )
766: if ( v == vn[i].v ) {
767: vn[i].n = 1;
768: break;
769: }
770: for ( dc = DC(p); dc; dc = NEXT(dc) )
771: markv(vn,n,COEF(dc));
772: }
773:
774: void vntovl(vn,n,vlp)
775: VN vn;
776: int n;
777: VL *vlp;
778: {
779: int i;
780: VL tvl,tvl0;
781:
782: for ( i = 0, tvl0 = 0; ; ) {
783: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
784: if ( i == n )
785: break;
786: else {
787: if ( !tvl0 ) {
788: NEWVL(tvl0);
789: tvl = tvl0;
790: } else {
791: NEWVL(NEXT(tvl));
792: tvl = NEXT(tvl);
793: }
794: tvl->v = vn[i++].v;
795: }
796: }
797: if ( tvl0 )
798: NEXT(tvl) = 0;
799: *vlp = tvl0;
800: }
801:
802: int dbound(v,f)
803: V v;
804: P f;
805: {
1.5 noro 806: int m,t;
1.1 noro 807: DCP dc;
808:
809: if ( !f )
810: return ( -1 );
811: else if ( v != VR(f) )
812: return homdeg(f);
813: else {
1.5 noro 814: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
815: t = homdeg(COEF(dc));
816: m = MAX(m,t);
817: }
1.1 noro 818: return ( m );
819: }
820: }
821:
822: int homdeg(f)
823: P f;
824: {
1.5 noro 825: int m,t;
1.1 noro 826: DCP dc;
827:
828: if ( !f )
829: return ( -1 );
830: else if ( NUM(f) )
831: return( 0 );
832: else {
1.5 noro 833: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
834: t = QTOS(DEG(dc))+homdeg(COEF(dc));
835: m = MAX(m,t);
836: }
1.1 noro 837: return ( m );
838: }
839: }
840:
841: int minhomdeg(f)
842: P f;
843: {
1.5 noro 844: int m,t;
1.1 noro 845: DCP dc;
846:
847: if ( !f )
848: return ( -1 );
849: else if ( NUM(f) )
850: return( 0 );
851: else {
1.5 noro 852: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
853: t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
854: m = MIN(m,t);
855: }
1.1 noro 856: return ( m );
857: }
858: }
859:
860: void adjc(vl,f,a,lc0,q,fr,ar)
861: VL vl;
862: P f,a,lc0;
863: Q q;
864: P *fr,*ar;
865: {
866: P m,m1;
867: Q t;
868:
869: invl((Q)LC(f),q,&t);
870: mulq((Q)lc0,t,(Q *)&m);
871: mulpq(f,m,&m1); cmp(q,m1,fr);
872: invl((Q)m,q,&t);
873: mulpq(a,(P)t,&m1);
874: cmp(q,m1,ar);
875: }
876:
877: #if 1
878: void affinemain(vl,p,v0,n,pl,pr)
879: VL vl;
880: V v0;
881: int n;
882: P *pl;
883: P p;
884: P *pr;
885: {
886: P x,t,m,c,s,a;
887: DCP dc;
888: Q d;
889:
890: if ( !p )
891: *pr = 0;
892: else if ( NUM(p) )
893: *pr = p;
894: else if ( VR(p) != v0 ) {
895: MKV(VR(p),x);
896: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
897: affinemain(vl,COEF(dc),v0,n,pl,&t);
898: if ( DEG(dc) ) {
899: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
900: addp(vl,m,c,&a); c = a;
901: } else {
902: addp(vl,t,c,&a); c = a;
903: }
904: }
905: *pr = c;
906: } else {
907: dc = DC(p);
908: c = COEF(dc);
909: for ( d = DEG(dc), dc = NEXT(dc);
910: dc; d = DEG(dc), dc = NEXT(dc) ) {
911: mulp(vl,pl[QTOS(d)-QTOS(DEG(dc))],c,&m);
912: addp(vl,m,COEF(dc),&c);
913: }
914: if ( d ) {
915: mulp(vl,pl[QTOS(d)],c,&m); c = m;
916: }
917: *pr = c;
918: }
919: }
920: #endif
921:
922: #if 0
923: affine(vl,p,vn,r)
924: VL vl;
925: P p;
926: VN vn;
927: P *r;
928: {
929: int n,d,d1,i;
930: Q *q;
931: Q **bc;
932:
933: if ( !p || NUM(p) )
934: *r = p;
935: else {
936: for ( i = 0, d = 0; vn[i].v; i++ )
937: d1 = getdeg(vn[i].v,p), d = MAX(d,d1);
938: W_CALLOC(d+1,Q *,bc);
939: for ( i = 0; i <= d; i++ )
940: W_CALLOC(i+1,Q,q), mkbc(i,q), bc[i] = q;
941: afmain(vl,bc,p,vn,r);
942: }
943: }
944:
945: afmain(vl,bc,p,vn,r)
946: VL vl;
947: Q **bc;
948: P p;
949: VN vn;
950: P *r;
951: {
952: P t,s,u;
953: P *c,*rc;
954: Q *q;
955: DCP dc;
956: int n,i,j;
957:
958: if ( !p || NUM(p) || !vn[0].v )
959: *r = p;
960: else if ( vn[0].v != VR(p) ) {
961: for ( i = 0; vn[i].v && (vn[i].v != VR(p)); i++ );
962: if ( vn[i].v )
963: afmain(vl,bc,p,vn+i,r);
964: else {
965: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
966: for ( dc = DC(p); dc; dc = NEXT(dc) )
967: afmain(vl,bc,COEF(dc),vn,&c[QTOS(DEG(dc))]);
968: plisttop(c,VR(p),n,r);
969: }
970: } else {
971: n = QTOS(DEG(DC(p))); W_CALLOC(n+1,P,c);
972: W_CALLOC(n+1,P,rc); W_CALLOC(n+1,Q,q);
973: for ( dc = DC(p); dc; dc = NEXT(dc) )
974: afmain(vl,bc,COEF(dc),vn+1,&c[QTOS(DEG(dc))]);
975: if ( !vn[0].n )
976: bcopy(c,rc,sizeof(P)*(n+1));
977: else {
978: for ( q[0] = ONE, STOQ(vn[0].n,q[1]), i = 2; i <= n; i++ )
979: mulq(q[i-1],q[1],&q[i]);
980: for ( j = 0; j <= n; rc[j] = t, j++ )
981: for ( i = j, t = 0; i <= n; i++ )
982: if ( c[i] )
983: mulq(bc[i][j],q[i-j],&s), mulp(vl,s,c[i],&u),
984: addp(CO,u,t,&s), t = s;
985: }
986: plisttop(rc,VR(p),n,r);
987: }
988: }
989: #endif
990:
991: void restore(vl,f,vn,fr)
992: VL vl;
993: P f;
994: VN vn;
995: P *fr;
996: {
997: int i;
998: P vv,g,g1,t;
999: Q s;
1000:
1001: g = f;
1002: for ( i = 0; vn[i].v; i++ ) {
1003: MKV(vn[i].v,t);
1004: if ( vn[i].n ) {
1005: STOQ(-vn[i].n,s);
1006: addp(vl,t,(P)s,&vv);
1007: } else
1008: vv = t;
1009:
1010: substp(vl,g,vn[i].v,vv,&g1); g = g1;
1011: }
1012: *fr = g;
1013: }
1014:
1015: void mergev(vl,vl1,vl2,nvlp)
1016: VL vl,vl1,vl2,*nvlp;
1017: {
1018: int i,n;
1019: VL tvl;
1020: VN vn;
1021:
1022: if ( !vl1 ) {
1023: *nvlp = vl2; return;
1024: } else if ( !vl2 ) {
1025: *nvlp = vl1; return;
1026: }
1027: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
1028: n = i;
1029: W_CALLOC(n,struct oVN,vn);
1030: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
1031: vn[i].v = tvl->v;
1032: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
1033: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
1034: i++;
1035: if ( i == n )
1036: break;
1037: else
1038: vn[i].n = 1;
1039: }
1040: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
1041: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
1042: i++;
1043: if ( i == n )
1044: break;
1045: else
1046: vn[i].n = 1;
1047: }
1048: vntovl(vn,n,nvlp);
1049: }
1050:
1051: #if 0
1052: void substvp(vl,f,vn,g)
1053: VL vl;
1054: P f;
1055: VN vn;
1056: P *g;
1057: {
1058: V v;
1059: int i;
1060: P h,h1;
1061: Q t;
1062:
1063: h = f;
1064: for ( i = 0; v = vn[i].v; i++ ) {
1065: STOQ(vn[i].n,t);
1066: substp(vl,h,v,(P)t,&h1); h = h1;
1067: }
1068: *g = h;
1069: }
1070:
1071: void affine(vl,f,vn,fr)
1072: VL vl;
1073: P f;
1074: VN vn;
1075: P *fr;
1076: {
1077: int i,j,n;
1078: P vv,g,g1,t,u;
1079: Q s;
1080: int *dlist;
1081: P **plist;
1082:
1083: for ( n = 0; vn[n].v; n++);
1084: dlist = (int *)ALLOCA((n+1)*sizeof(int));
1085: plist = (P **)ALLOCA((n+1)*sizeof(P *));
1086: for ( i = 0; vn[i].v; i++ ) {
1087: if ( !vn[i].n )
1088: continue;
1089: dlist[i] = getdeg(vn[i].v,f);
1090: plist[i] = (P *)ALLOCA((dlist[i]+1)*sizeof(P));
1091:
1092: MKV(vn[i].v,t);
1093: if ( vn[i].n ) {
1094: STOQ(vn[i].n,s); addp(vl,t,(P)s,&vv);
1095: } else
1096: vv = t;
1097:
1098: for ( j = 0, t = (P)ONE; j < dlist[i]; j++ ) {
1099: plist[i][j] = t;
1100: mulp(vl,t,vv,&u);
1101: t = u;
1102: }
1103: plist[i][j] = t;
1104: }
1105:
1106: g = f;
1107: for ( i = 0; vn[i].v; i++ ) {
1108: if ( !vn[i].n )
1109: continue;
1110: affinemain(vl,g,vn[i].v,dlist[i],plist[i],&g1); g = g1;
1111: }
1112: *fr = g;
1113: }
1114: #endif