Annotation of OpenXM_contrib2/asir2018/engine/E.c, Revision 1.2
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
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.2 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/E.c,v 1.1 2018/09/19 05:45:06 noro Exp $
1.1 noro 49: */
50: #include "ca.h"
51:
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,Z q,int k,P *gp,P *hp)
53: {
54: P g1,h1,a1,b1;
55: Z q2,r2,two;
56:
1.2 ! noro 57: STOZ(2,two);
1.1 noro 58: divqrz(q,two,&q2,&r2);
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:
63: void henmvmain(VL vl,VN vn,P f,P fi0,P fi1,P gi0,P gi1,P l0,P l1,Z mod,Z mod2,int k,P *fr0,P *fr1)
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: invz((Z)LC(fi1),mod,(Z *)&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:
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,Z *qrp)
156: {
157: Z q,qq,q2,two,r2;
158: int n,i;
159: UM wg0,wg1,wf0,wf1;
160: P f0,f1,g0,g1,m,m1,s,c,d,q1,r1,rm,rm1,a,a1,a2;
161:
162: n = UDEG(f);
163: wf0 = W_UMALLOC(n); wf1 = W_UMALLOC(n);
164: wg0 = W_UMALLOC(n); wg1 = W_UMALLOC(n);
165: cpyum(fi0,wf0); cpyum(fi1,wf1);
166: eucum(p,wf0,wf1,wg1,wg0);
167: umtop(VR(f),fi0,&f0); umtop(VR(f),fi1,&f1);
168: umtop(VR(f),wg0,&g0); umtop(VR(f),wg1,&g1);
169:
1.2 ! noro 170: STOZ(2,two); STOZ(p,q); divqrz(q,two,&q2,&r2);
1.1 noro 171: for ( i = 1; i < k; i++ ) {
172: /* c = ((f - f0*f1)/q) mod q;
173: q1 = (c*g1) / f1;
174: r1 = (c*g1) mod f1;
175: f1 += (r1 mod q)*q;
176: f0 += ((c*g0 + q1*f0) mod q)*q;
177:
178: d = ((1 - (f1*g0 + f0*g1))/q) mod q;
179: q1 = (d*g0) / f1;
180: r1 = (d*g0) mod f1;
181: g1 += (r1 mod q)*q;
182: g0 += ((c*g0 + q1*f0) mod q)*q;
183: q = q^2;
184: */
185:
186: /* c = ((f - f0*f1)/q) mod q */
187: mulp(CO,f0,f1,&m); subp(CO,f,m,&s);
188: divcp(s,(Q)q,&m); cm2p(q,q2,m,&c);
189:
190: /* q1 = (c*g1) / f1; r1 = (c*g1) mod f1; */
191: mulp(CO,c,g1,&m); cm2p(q,q2,m,&s);
192: udivpwm(q,s,f1,&q1,&r1);
193:
194: /* f1 = f1 + (r1 mod q)*q; */
195: cm2p(q,q2,r1,&rm);
196: mulpq(rm,(P)q,&rm1); addp(CO,f1,rm1,&a);
197: f1 = a;
198:
199: /* a1 = (c*g0 + q1*f0) mod q; */
200: mulp(CO,c,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
201: cm2p(q,q2,a,&a1);
202:
203: /* f0 = f0 + a1*q; */
204: mulpq(a1,(P)q,&a2);
205: addp(CO,f0,a2,&a);
206: f0 = a;
207:
208: /* d = ((1 - (f1*g0 + f0*g1))/q) mod q; */
209: mulp(CO,f1,g0,&m); mulp(CO,f0,g1,&m1); addp(CO,m,m1,&a);
210: subp(CO,(P)ONE,a,&s);
211: divcp(s,(Q)q,&m); cm2p(q,q2,m,&d);
212:
213: /* q1 = (d*g0) / f1; r1 = (d*g0) mod f1; */
214: mulp(CO,d,g1,&m); udivpwm(q,m,f1,&q1,&r1);
215:
216: /* g1 = g1 + (r1 mod q )*q; */
217: cm2p(q,q2,r1,&rm);
218: mulpq(rm,(P)q,&rm1); addp(CO,g1,rm1,&a);
219: g1 = a;
220:
221: /* g0 = g0 + (c*g0 + q1*f0 mod q)*q; */
222: mulp(CO,d,g0,&m); mulp(CO,q1,f0,&m1); addp(CO,m,m1,&a);
223: cm2p(q,q2,a,&a1); mulpq(a1,(P)q,&a2);
224: addp(CO,g0,a2,&a);
225: g0 = a;
226:
227: /* q = q^2; */
228: mulz(q,q,&qq); q = qq; divqrz(q,two,&q2,&r2);
229: }
230: *fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
231: }
232:
233: void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Z *qp)
234: {
235: V v;
236: Z w,q,q1;
237: Q f;
238: Q u,t,a,b,s,ts;
239: P c,c1;
240: P tg,th,tr;
241: UM wg,wh,ws,wt,wm;
242: int n,m,modres,mod,index,i;
243: P gc0,hc0;
244: P z,zz,zzz;
245:
246:
247: v = VR(g); n=deg(v,g); m=deg(v,h);
248: norm(g,&a); norm(h,&b);
1.2 ! noro 249: STOZ(m,w); pwrq(a,(Q)w,&t);
! 250: STOZ(n,w); pwrq(b,(Q)w,&s);
1.1 noro 251: mulq(t,s,&ts);
252:
253: factorialz(n+m,&w); mulq(ts,(Q)w,&s);
254: addq(s,s,&f);
255:
256: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
257: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
258: wm = W_UMALLOC(m+n);
259:
260: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
261: mod = get_lprime(index++);
262: if ( !remqi((Q)LC(g),mod) || !remqi((Q)LC(h),mod) )
263: continue;
264: ptomp(mod,g,&tg); ptomp(mod,h,&th);
265: srchump(mod,tg,th,&tr);
266: if ( !tr )
267: continue;
268: else
269: modres = CONT((MQ)tr);
270:
271: mptoum(tg,wg); mptoum(th,wh);
272: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
273: DEG(wm) = 0; COEF(wm)[0] = modres;
274: mulum(mod,ws,wm,wt);
275: for ( i = DEG(wt); i >= 0; i-- )
276: if ( ( COEF(wt)[i] * 2 ) > mod )
277: COEF(wt)[i] -= mod;
278: chnrem(mod,v,c,q,wt,&c1,&q1);
279: if ( !ucmpp(c,c1) ) {
280: mulp(CO,c,g,&z); divsrp(CO,z,h,&zz,&zzz);
281: if ( NUM(zzz) ) {
282: q = q1; c = c1;
283: break;
284: }
285: }
286: q = q1; c = c1;
287:
288: if ( cmpq(f,(Q)q) < 0 )
289: break;
290: }
291: ptozp(c,1,&s,&gc0);
292: /* gc0*g = zz*h + zzz -> c*g + (-zz)*h = zzz */
293: mulp(CO,gc0,g,&z);
294: divsrp(CO,z,h,&zz,&zzz);
295: ptozp(zz,1,&s,(P *)&t);
296: if ( INT(s) )
297: chsgnp(zz,&hc0);
298: else {
299: MPZTOZ(mpq_denref(BDY(s)),q); mulq((Q)q,(Q)zzz,(Q *)&q1); zzz = (P)q1;
300: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,zz,(P)q,&c); chsgnp(c,&hc0);
301: }
302: if ( !INT((Q)zzz) ) {
303: MPZTOZ(mpq_denref(BDY((Q)zzz)),q); MPZTOZ(mpq_numref(BDY((Q)zzz)),q1); zzz = (P)q1;
304: mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
305: }
306: for ( index = 0; ; ) {
307: mod = get_lprime(index++);
308: if ( !remqi((Q)zzz,mod) ||
309: !remqi((Q)LC(g),mod) ||
310: !remqi((Q)LC(h),mod) )
311: continue;
1.2 ! noro 312: for ( STOZ(mod,q); cmpq((Q)q,bound) < 0; ) {
1.1 noro 313: mulz(q,q,&q1); q = q1;
314: }
315: *qp = q;
316: invz((Z)zzz,q,&q1);
317: mulp(CO,gc0,(P)q1,hcp); mulp(CO,hc0,(P)q1,gcp);
318: return;
319: }
320: }
321:
322: void addm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
323: {
324: P t;
325:
326: addp(vl,n1,n2,&t);
327: if ( !t )
328: *nr = 0;
329: else
330: cm2p(mod,mod2,t,nr);
331: }
332:
333: void subm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
334: {
335: P t;
336:
337: subp(vl,n1,n2,&t);
338: if ( !t )
339: *nr = 0;
340: else
341: cm2p(mod,mod2,t,nr);
342: }
343:
344: void mulm2p(VL vl,Z mod,Z mod2,P n1,P n2,P *nr)
345: {
346: P t;
347:
348: mulp(vl,n1,n2,&t);
349: if ( !t )
350: *nr = 0;
351: else
352: cm2p(mod,mod2,t,nr);
353: }
354:
355: void cmp(Z mod,P p,P *pr)
356: {
357: P t;
358: DCP dc,dcr,dcr0;
359:
360: if ( !p )
361: *pr = 0;
362: else if ( NUM(p) )
363: remz((Z)p,mod,(Z *)pr);
364: else {
365: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
366: cmp(mod,COEF(dc),&t);
367: if ( t ) {
368: NEXTDC(dcr0,dcr);
369: DEG(dcr) = DEG(dc);
370: COEF(dcr) = t;
371: }
372: }
373: if ( !dcr0 )
374: *pr = 0;
375: else {
376: NEXT(dcr) = 0;
377: MKP(VR(p),dcr0,*pr);
378: }
379: }
380: }
381:
382: void cm2p(Z mod,Z m,P p,P *pr)
383: {
384: P t;
385: DCP dc,dcr,dcr0;
386:
387: if ( !p )
388: *pr = 0;
389: else if ( NUM(p) )
390: rem2q((Z)p,mod,m,(Z *)pr);
391: else {
392: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
393: cm2p(mod,m,COEF(dc),&t);
394: if ( t ) {
395: NEXTDC(dcr0,dcr);
396: DEG(dcr) = DEG(dc);
397: COEF(dcr) = t;
398: }
399: }
400: if ( !dcr0 )
401: *pr = 0;
402: else {
403: NEXT(dcr) = 0;
404: MKP(VR(p),dcr0,*pr);
405: }
406: }
407: }
408:
409: void addm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
410: {
411: Z t;
412:
413: addz(n1,n2,&t);
414: if ( !t )
415: *nr = 0;
416: else
417: rem2q(t,mod,mod2,nr);
418: }
419:
420: void subm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
421: {
422: Z t;
423:
424: subz(n1,n2,&t);
425: if ( !t )
426: *nr = 0;
427: else
428: rem2q(t,mod,mod2,nr);
429: }
430:
431: void mulm2q(Z mod,Z mod2,Z n1,Z n2,Z *nr)
432: {
433: Z t;
434:
435: mulz(n1,n2,&t);
436: if ( !t )
437: *nr = 0;
438: else
439: rem2q(t,mod,mod2,nr);
440: }
441:
442: void rem2q(Z n,Z m,Z m2,Z *nr)
443: {
444: Z r;
445: int sgn;
446:
447: remz(n,m,&r);
448: if ( !r )
449: *nr = 0;
450: else {
451: sgn = cmpz(r,m2);
452: if ( sgn > 0 )
453: subz(r,m,nr);
454: else
455: *nr = r;
456: }
457: }
458:
459: /*
460: extract d-homogeneous part with respect to vl - {v}
461: */
462:
463: void exthpc_generic(VL vl,P p,int d,V v,P *pr)
464: {
465: P w,x,t,t1,a,xd;
466: V v0;
467: DCP dc;
468:
469: if ( d < 0 || !p )
470: *pr = 0;
471: else if ( NUM(p) )
472: if ( d == 0 )
473: *pr = p;
474: else
475: *pr = 0;
476: else if ( v == VR(p) )
477: exthpc(vl,v,p,d,pr);
478: else {
479: v0 = VR(p);
480: for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
1.2 ! noro 481: exthpc_generic(vl,COEF(dc),d-ZTOS(DEG(dc)),v,&t);
1.1 noro 482: pwrp(vl,x,DEG(dc),&xd);
483: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
484: }
485: *pr = w;
486: }
487: }
488:
489: void exthp(VL vl,P p,int d,P *pr)
490: {
491: P t,t1,a,w,x,xd;
492: DCP dc;
493:
494: if ( d < 0 )
495: *pr = 0;
496: else if ( NUM(p) )
497: if ( d == 0 )
498: *pr = p;
499: else
500: *pr = 0;
501: else {
502: for ( MKV(VR(p),x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
1.2 ! noro 503: exthp(vl,COEF(dc),d - ZTOS(DEG(dc)),&t);
1.1 noro 504: pwrp(vl,x,DEG(dc),&xd);
505: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
506: }
507: *pr = w;
508: }
509: }
510:
511: void exthpc(VL vl,V v,P p,int d,P *pr)
512: {
513: P t,t1,a,w,x,xd;
514: DCP dc;
515:
516: if ( v != VR(p) )
517: exthp(vl,p,d,pr);
518: else if ( d < 0 )
519: *pr = 0;
520: else {
521: for ( MKV(v,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
522: exthp(vl,COEF(dc),d,&t);
523: pwrp(vl,x,DEG(dc),&xd);
524: mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
525: }
526: *pr = w;
527: }
528: }
529:
530: void cbound(VL vl,P p,Q *b)
531: {
532: Q n,m;
533: Z t,e;
534: int k;
535:
536: cmax(p,&n);
537: addq(n,n,&m);
538:
539: k = geldb(vl,p);
1.2 ! noro 540: STOZ(3,t); STOZ(k,e);
1.1 noro 541:
542: pwrq((Q)t,(Q)e,&n);
543: mulq(m,n,b);
544: }
545:
546: int geldb(VL vl,P p)
547: {
548: int m;
549:
550: for ( m = 0; vl; vl = NEXT(vl) )
551: m += getdeg(vl->v,p);
552: return ( m );
553: }
554:
555: int getdeg(V v,P p)
556: {
557: int m,t;
558: DCP dc;
559:
560: if ( !p || NUM(p) )
561: return ( 0 );
562: else if ( v == VR(p) )
563: return ( deg(v,p) );
564: else {
565: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
566: t = getdeg(v,COEF(dc));
567: m = MAX(m,t);
568: }
569: return ( m );
570: }
571: }
572:
573: /* YYY */
574:
575: void cmax(P p,Q *b)
576: {
577: DCP dc;
578: Q m,m1;
579:
580: if ( NUM(p) )
581: absq((Q)p,b);
582: else {
583: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
584: cmax(COEF(dc),&m1);
585: if ( cmpq(m1,m) > 0 )
586: m = m1;
587: }
588: *b = m;
589: }
590: }
591:
592: int nextbin(VN vn,int n)
593: {
594: int tmp,i,carry;
595:
596: if ( n == 0 )
597: return ( 1 );
598:
599: for ( i = n - 1, carry = 1; i >= 0; i-- ) {
600: tmp = vn[i].n + carry;
601: vn[i].n = tmp % 2;
602: carry = tmp / 2;
603: }
604: return ( carry );
605: }
606:
607: void mulsgn(VN vn,VN vnt,int n,VN vn1)
608: {
609: int i;
610:
611: for ( i = 0; vn[i].v; i++ )
612: vn1[i].n = vn[i].n;
613: for ( i = 0; i < n; i++ )
614: if ( vnt[i].n )
615: vn1[(long)vnt[i].v].n *= -1;
616: }
617:
618: void next(VN vn)
619: {
620: int i,m,n,tmp,carry;
621:
622: for ( m = 0, i = 0; vn[i].v; i++ )
623: m = MAX(m,ABS(vn[i].n));
624: if ( m == 0 ) {
625: vn[--i].n = 1;
626: return;
627: }
628: for ( carry = 1, m++, i--; ( i >= 0 ) && carry; i-- ) {
629: tmp = vn[i].n + carry;
630: vn[i].n = tmp % m;
631: carry = tmp / m;
632: }
633: if ( ( i == -1 ) && carry ) {
634: for ( i = 0; vn[i].v; i++ )
635: vn[i].n = 0;
636: vn[--i].n = m;
637: } else {
638: for ( n = 0, i = 0; vn[i].v; i++ )
639: n = MAX(n,ABS(vn[i].n));
640: if ( n < m - 1 )
641: vn[--i].n = m - 1;
642: }
643: }
644:
645: void clctv(VL vl,P p,VL *nvlp)
646: {
647: int i,n;
648: VL tvl;
649: VN tvn;
650:
651: if ( !p || NUM(p) ) {
652: *nvlp = 0;
653: return;
654: }
655:
656: for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
657: tvn = (VN) ALLOCA((n+1)*sizeof(struct oVN));
658: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ ) {
659: tvn[i].v = tvl->v;
660: tvn[i].n = 0;
661: }
662:
663: markv(tvn,n,p);
664: vntovl(tvn,n,nvlp);
665: }
666:
667: void markv(VN vn,int n,P p)
668: {
669: V v;
670: DCP dc;
671: int i;
672:
673: if ( NUM(p) )
674: return;
675: v = VR(p);
676: for ( i = 0, v = VR(p); i < n; i++ )
677: if ( v == vn[i].v ) {
678: vn[i].n = 1;
679: break;
680: }
681: for ( dc = DC(p); dc; dc = NEXT(dc) )
682: markv(vn,n,COEF(dc));
683: }
684:
685: void vntovl(VN vn,int n,VL *vlp)
686: {
687: int i;
688: VL tvl,tvl0;
689:
690: for ( i = 0, tvl0 = 0; ; ) {
691: while ( ( i < n ) && ( vn[i].n == 0 ) ) i++;
692: if ( i == n )
693: break;
694: else {
695: if ( !tvl0 ) {
696: NEWVL(tvl0);
697: tvl = tvl0;
698: } else {
699: NEWVL(NEXT(tvl));
700: tvl = NEXT(tvl);
701: }
702: tvl->v = vn[i++].v;
703: }
704: }
705: if ( tvl0 )
706: NEXT(tvl) = 0;
707: *vlp = tvl0;
708: }
709:
710: int dbound(V v,P f)
711: {
712: int m,t;
713: DCP dc;
714:
715: if ( !f )
716: return ( -1 );
717: else if ( v != VR(f) )
718: return homdeg(f);
719: else {
720: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
721: t = homdeg(COEF(dc));
722: m = MAX(m,t);
723: }
724: return ( m );
725: }
726: }
727:
728: int homdeg(P f)
729: {
730: int m,t;
731: DCP dc;
732:
733: if ( !f )
734: return ( -1 );
735: else if ( NUM(f) )
736: return( 0 );
737: else {
738: for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
1.2 ! noro 739: t = ZTOS(DEG(dc))+homdeg(COEF(dc));
1.1 noro 740: m = MAX(m,t);
741: }
742: return ( m );
743: }
744: }
745:
746: int minhomdeg(P f)
747: {
748: int m,t;
749: DCP dc;
750:
751: if ( !f )
752: return ( -1 );
753: else if ( NUM(f) )
754: return( 0 );
755: else {
756: for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
1.2 ! noro 757: t = ZTOS(DEG(dc))+minhomdeg(COEF(dc));
1.1 noro 758: m = MIN(m,t);
759: }
760: return ( m );
761: }
762: }
763:
764: void adjc(VL vl,P f,P a,P lc0,Z q,P *fr,P *ar)
765: {
766: Z m,t;
767: P m1;
768:
769: invz((Z)LC(f),q,&t);
770: mulz((Z)lc0,t,&m);
771: mulpq(f,(P)m,&m1); cmp(q,m1,fr);
772: invz(m,q,&t);
773: mulpq(a,(P)t,&m1);
774: cmp(q,m1,ar);
775: }
776:
777: void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr)
778: {
779: P x,t,m,c,s,a;
780: DCP dc;
781: Z d;
782:
783: if ( !p )
784: *pr = 0;
785: else if ( NUM(p) )
786: *pr = p;
787: else if ( VR(p) != v0 ) {
788: MKV(VR(p),x);
789: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
790: affinemain(vl,COEF(dc),v0,n,pl,&t);
791: if ( DEG(dc) ) {
792: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
793: addp(vl,m,c,&a); c = a;
794: } else {
795: addp(vl,t,c,&a); c = a;
796: }
797: }
798: *pr = c;
799: } else {
800: dc = DC(p);
801: c = COEF(dc);
802: for ( d = DEG(dc), dc = NEXT(dc);
803: dc; d = DEG(dc), dc = NEXT(dc) ) {
1.2 ! noro 804: mulp(vl,pl[ZTOS(d)-ZTOS(DEG(dc))],c,&m);
1.1 noro 805: addp(vl,m,COEF(dc),&c);
806: }
807: if ( d ) {
1.2 ! noro 808: mulp(vl,pl[ZTOS(d)],c,&m); c = m;
1.1 noro 809: }
810: *pr = c;
811: }
812: }
813:
814: void restore(VL vl,P f,VN vn,P *fr)
815: {
816: int i;
817: P vv,g,g1,t;
818: Z s;
819:
820: g = f;
821: for ( i = 0; vn[i].v; i++ ) {
822: MKV(vn[i].v,t);
823: if ( vn[i].n ) {
1.2 ! noro 824: STOZ(-vn[i].n,s);
1.1 noro 825: addp(vl,t,(P)s,&vv);
826: } else
827: vv = t;
828:
829: substp(vl,g,vn[i].v,vv,&g1); g = g1;
830: }
831: *fr = g;
832: }
833:
834: void mergev(VL vl,VL vl1,VL vl2,VL *nvlp)
835: {
836: int i,n;
837: VL tvl;
838: VN vn;
839:
840: if ( !vl1 ) {
841: *nvlp = vl2; return;
842: } else if ( !vl2 ) {
843: *nvlp = vl1; return;
844: }
845: for ( i = 0, tvl = vl; tvl; i++, tvl = NEXT(tvl) );
846: n = i;
847: W_CALLOC(n,struct oVN,vn);
848: for ( i = 0, tvl = vl; tvl; tvl = NEXT(tvl), i++ )
849: vn[i].v = tvl->v;
850: for ( i = 0, tvl = vl1; (i < n) && tvl; i++, tvl = NEXT(tvl)) {
851: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
852: i++;
853: if ( i == n )
854: break;
855: else
856: vn[i].n = 1;
857: }
858: for ( i = 0, tvl = vl2; (i < n) && tvl; i++, tvl = NEXT(tvl) ) {
859: while ( ( i < n ) && ( vn[i].v != tvl->v ) )
860: i++;
861: if ( i == n )
862: break;
863: else
864: vn[i].n = 1;
865: }
866: vntovl(vn,n,nvlp);
867: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / engine