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