![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Mgfs.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 Mgfs.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.4 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.3 2001/06/25 01:35:21 noro Exp $ */
1.1 noro 2:
3: #include "ca.h"
4:
5: void mulssfum(UM,int,UM);
6:
7: void addsfum(p1,p2,pr)
8: UM p1,p2,pr;
9: {
10: int *c1,*c2,*cr,i,dmax,dmin;
11:
12: if ( DEG(p1) == -1 ) {
13: cpyum(p2,pr);
14: return;
15: }
16: if ( DEG(p2) == -1 ) {
17: cpyum(p1,pr);
18: return;
19: }
20: if ( DEG(p1) >= DEG(p2) ) {
21: c1 = COEF(p1); c2 = COEF(p2); dmax = DEG(p1); dmin = DEG(p2);
22: } else {
23: c1 = COEF(p2); c2 = COEF(p1); dmax = DEG(p2); dmin = DEG(p1);
24: }
25: for ( i = 0, cr = COEF(pr); i <= dmin; i++ )
26: cr[i] = _addsf(c1[i],c2[i]);
27: for ( ; i <= dmax; i++ )
28: cr[i] = c1[i];
29: if ( dmax == dmin )
30: degum(pr,dmax);
31: else
32: DEG(pr) = dmax;
33: }
34:
35: void subsfum(p1,p2,pr)
36: UM p1,p2,pr;
37: {
38: int *c1,*c2,*cr,i;
39: int dmax,dmin;
40:
41: if ( DEG(p1) == -1 ) {
42: for ( i = DEG(pr) = DEG(p2), c2 = COEF(p2), cr = COEF(pr);
43: i >= 0; i-- )
44: cr[i] = _chsgnsf(c2[i]);
45: return;
46: }
47: if ( DEG(p2) == -1 ) {
48: cpyum(p1,pr);
49: return;
50: }
51: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
52: if ( DEG(p1) >= DEG(p2) ) {
53: dmax = DEG(p1); dmin = DEG(p2);
54: for ( i = 0; i <= dmin; i++ )
55: cr[i] = _subsf(c1[i],c2[i]);
56: for ( ; i <= dmax; i++ )
57: cr[i] = c1[i];
58: } else {
59: dmax = DEG(p2); dmin = DEG(p1);
60: for ( i = 0; i <= dmin; i++ )
61: cr[i] = _subsf(c1[i],c2[i]);
62: for ( ; i <= dmax; i++ )
63: cr[i] = _chsgnsf(c2[i]);
64: }
65: if ( dmax == dmin )
66: degum(pr,dmax);
67: else
68: DEG(pr) = dmax;
69: }
70:
71: void gcdsfum(p1,p2,pr)
72: UM p1,p2,pr;
73: {
74: int inv;
75: UM t1,t2,q,tum;
76: int drem;
77:
78: if ( DEG(p1) < 0 )
79: cpyum(p2,pr);
80: else if ( DEG(p2) < 0 )
81: cpyum(p1,pr);
82: else {
83: if ( DEG(p1) >= DEG(p2) ) {
84: t1 = p1; t2 = p2;
85: } else {
86: t1 = p2; t2 = p1;
87: }
88: q = W_UMALLOC(DEG(t1));
89: while ( ( drem = divsfum(t1,t2,q) ) >= 0 ) {
90: tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;
91: }
92: inv = _invsf(COEF(t2)[DEG(t2)]);
93: mulssfum(t2,inv,pr);
94: }
95: }
96: void mulsfum(p1,p2,pr)
97: UM p1,p2,pr;
98: {
99: int *pc1,*pcr;
100: int *c1,*c2,*cr;
101: int mul;
102: int i,j,d1,d2;
103:
104: if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
105: DEG(pr) = -1;
106: return;
107: }
108: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
109: bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
110: for ( i = 0; i <= d2; i++, cr++ )
111: if ( mul = *c2++ )
112: for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )
113: *pcr = _addsf(_mulsf(*pc1,mul),*pcr);
114: DEG(pr) = d1 + d2;
115: }
116:
117: void mulssfum(p,n,pr)
118: int n;
119: UM p,pr;
120: {
121: int *sp,*dp;
122: int i;
123:
124: for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
125: i >= 0; i--, dp--, sp-- )
126: *dp = _mulsf(*sp,n);
127: }
128:
129: int divsfum(p1,p2,pq)
130: UM p1,p2,pq;
131: {
132: int *pc1,*pct;
133: int *c1,*c2,*ct;
134: int inv,hd,tmp;
135: int i,j, d1,d2,dd;
136:
137: if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
138: DEG(pq) = -1;
139: return d1;
140: }
141: c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
142: if ( ( hd = c2[d2] ) != _onesf() ) {
143: inv = _invsf(hd);
144: for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
145: *pc1 = _mulsf(*pc1,inv);
146: } else
147: inv = _onesf();
148: for ( i = dd, ct = c1+d1; i >= 0; i-- )
149: if ( tmp = *ct-- ) {
150: tmp = _chsgnsf(tmp);
151: for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
152: *pct = _addsf(_mulsf(*pc1,tmp),*pct);
153: }
154: if ( inv != _onesf() ) {
155: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
156: *pc1 = _mulsf(*pc1,inv);
157: for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )
158: *pc1 = _mulsf(*pc1,hd);
159: }
160: for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
161: for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
162: *pct-- = *pc1--;
163: return i;
164: }
165:
166: void diffsfum(f,fd)
167: UM f,fd;
168: {
169: int *dp,*sp;
170: int i;
171:
172: for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
173: i >= 1; i--, dp--, sp-- ) {
174: *dp = _mulsf(*sp,_itosf(i));
175: }
176: degum(fd,DEG(f) - 1);
177: }
178:
179: void monicsfum(f)
180: UM f;
181: {
182: int *sp;
183: int i,inv;
184:
185: i = DEG(f); sp = COEF(f)+i;
186: inv = _invsf(*sp);
187: for ( ; i >= 0; i--, sp-- )
188: *sp = _mulsf(*sp,inv);
1.2 noro 189: }
190:
1.3 noro 191: void addsfarray(int,int *,int *);
192: void mulsfarray_trunc(int,int *,int *,int *);
1.2 noro 193:
194: void mulsflum(n,f1,f2,fr)
195: int n;
196: LUM f1,f2,fr;
197: {
198: int max;
199: int i,j,**p1,**p2,*px;
200: int *w,*w1,*w2;
201:
202: p1 = (int **)COEF(f1); p2 = (int **)COEF(f2);
203: w = W_ALLOC(2*(n+1)); w1 = W_ALLOC(DEG(f1)); w2 = W_ALLOC(DEG(f2));
204: for ( i = DEG(f1); i >= 0; i-- ) {
205: for ( j = n - 1, px = p1[i]; ( j >= 0 ) && ( px[j] == 0 ); j-- );
206: w1[i] = ( j == -1 ? 0 : 1 );
207: }
208: for ( i = DEG(f2); i >= 0; i-- ) {
209: for ( j = n - 1, px = p2[i]; ( j >= 0 ) && ( px[j] == 0 ); j-- );
210: w2[i] = ( j == -1 ? 0 : 1 );
211: }
212: for ( j = DEG(fr) = DEG(f1) + DEG(f2); j >= 0; j-- ) {
213: for ( i = n - 1, px = COEF(fr)[j]; i >= 0; i-- )
214: px[i] = 0;
215: for ( max = MIN(DEG(f1),j), i = MAX(0,j-DEG(f2)); i <= max; i++ )
216: if ( w1[i] != 0 && w2[j - i] != 0 ) {
1.3 noro 217: mulsfarray_trunc(n,p1[i],p2[j - i],w); addsfarray(n,w,px);
1.2 noro 218: }
219: }
220: }
221:
1.4 ! noro 222: /* f1 = f1->c[0]+f1->c[1]*y+..., f2 = f2->c[0]+f2->c[1]*y+... mod y^n */
! 223:
! 224: void mulsfbm(bound,f1,f2,fr)
! 225: int bound;
! 226: BM f1,f2,fr;
! 227: {
! 228: UM mul,t,s;
! 229: int i,j,h,d1,d2;
! 230:
! 231: if ( DEG(f1) < bound || DEG(f2) < bound )
! 232: error("mulsfbm : invalid input");
! 233:
! 234: d1 = DEG(COEF(f1)[0]);
! 235: for ( i = 1; i < bound; i++ )
! 236: d1 = MAX(DEG(COEF(f1)[i]),d1);
! 237: d2 = DEG(COEF(f2)[0]);
! 238: for ( i = 1; i < bound; i++ )
! 239: d2 = MAX(DEG(COEF(f2)[i]),d2);
! 240: t = W_UMALLOC(d1+d2);
! 241: s = W_UMALLOC(d1+d2);
! 242:
! 243: for ( i = 0; i < bound; i++ ) {
! 244: mul = COEF(f2)[i];
! 245: if ( DEG(mul) >= 0 )
! 246: for ( j = 0; i+j < bound; j++ ) {
! 247: if ( COEF(f1)[j] ) {
! 248: mulsfum(COEF(f1)[j],mul,t);
! 249: addsfum(t,COEF(fr)[i+j],s);
! 250: cpyum(s,COEF(fr)[i+j]);
! 251: }
! 252: }
! 253: }
! 254: DEG(fr) = bound;
! 255: }
! 256:
! 257: /* g += f */
! 258:
! 259: void addtosfbm(bound,f,g)
! 260: int bound;
! 261: BM f,g;
! 262: {
! 263: int i,d1,d2;
! 264: UM t;
! 265:
! 266: d1 = DEG(COEF(f)[0]);
! 267: for ( i = 1; i < bound; i++ )
! 268: d1 = MAX(DEG(COEF(f)[i]),d1);
! 269: d2 = DEG(COEF(g)[0]);
! 270: for ( i = 1; i < bound; i++ )
! 271: d2 = MAX(DEG(COEF(g)[i]),d2);
! 272: t = W_UMALLOC(MAX(d1,d2));
! 273: for ( i = 0; i < bound; i++ ) {
! 274: addsfum(COEF(f)[i],COEF(g)[i],t);
! 275: cpyum(t,COEF(g)[i]);
! 276: }
! 277: }
! 278:
1.3 noro 279: void addsfarray(n,a1,a2)
1.2 noro 280: int n;
281: int *a1,*a2;
282: {
1.3 noro 283: int i;
284:
285: for ( i = 0; i < n; i++, a1++, a2++ )
286: if ( *a1 )
287: *a2 = _addsf(*a1,*a2);
1.2 noro 288: }
289:
1.3 noro 290: void mulsfarray_trunc(n,a1,a2,r)
1.2 noro 291: int n;
292: int *a1,*a2,*r;
293: {
1.3 noro 294: int mul,i,j;
295: int *c1,*c2,*cr;
296: int *pc1,*pc2,*pcr;
297:
298: bzero(r,n*sizeof(int));
299: c1 = a1; c2 = a2; cr = r;
300: for ( i = 0; i < n; i++, cr++ ) {
301: if ( mul = *c2++ )
302: for ( j = 0, pc1 = c1, pcr = cr; j+i < n; j++, pc1++, pcr++ )
303: if ( *pc1 )
304: *pcr = _addsf(_mulsf(*pc1,mul),*pcr);
305: }
1.2 noro 306: }
307:
308: void eucsfum(f1,f2,a,b)
309: UM f1,f2,a,b;
310: {
311: UM g1,g2,a1,a2,a3,wm,q,tum;
312: int d,dr;
313:
314: /* XXX */
315: d = DEG(f1) + DEG(f2) + 10;
316: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);
317: a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);
318: q = W_UMALLOC(d);
1.3 noro 319: DEG(a1) = 0; COEF(a1)[0] = _onesf(); DEG(a2) = -1;
1.2 noro 320: cpyum(f1,g1); cpyum(f2,g2);
321: while ( 1 ) {
322: dr = divsfum(g1,g2,q); tum = g1; g1 = g2; g2 = tum;
323: if ( ( DEG(g2) = dr ) == -1 )
324: break;
325: mulsfum(a2,q,wm); subsfum(a1,wm,a3); dr = divsfum(a3,f2,q);
326: tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;
327: }
1.3 noro 328: if ( COEF(g1)[0] != _onesf() )
1.2 noro 329: mulssfum(a2,_invsf(COEF(g1)[0]),a);
330: else
331: cpyum(a2,a);
332: mulsfum(a,f1,wm);
333: if ( DEG(wm) >= 0 )
334: COEF(wm)[0] = _subsf(COEF(wm)[0],_onesf());
335: else {
336: DEG(wm) = 0; COEF(wm)[0] = _chsgnsf(_onesf());
337: }
338: divsfum(wm,f2,q); mulssfum(q,_chsgnsf(_onesf()),b);
1.3 noro 339: }
340:
341: void shiftsfum(UM,int,UM);
342:
343: void shiftsflum(n,f,ev)
344: int n;
345: LUM f;
346: int ev;
347: {
348: int d,i,j;
349: UM w,w1;
350: int *coef;
351: int **c;
352:
353: d = f->d;
354: c = f->c;
355: w = W_UMALLOC(n);
356: w1 = W_UMALLOC(n);
357: for ( i = 0; i <= d; i++ ) {
358: coef = c[i];
359: for ( j = n-1; j >= 0 && coef[j] == 0; j-- );
360: if ( j <= 0 )
361: continue;
362: DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));
363: shiftsfum(w,ev,w1);
364: bcopy(COEF(w1),coef,(j+1)*sizeof(int));
365: }
366: }
367:
368: /* f(x) -> g(x) = f(x+a) */
369:
370: void shiftsfum(f,a,g)
371: UM f;
372: int a;
373: UM g;
374: {
375: int n,i;
376: UM pwr,xa,w,t;
377: int *c;
378:
379: n = DEG(f);
380: if ( n <= 0 )
381: cpyum(f,g);
382: else {
383: c = COEF(f);
384:
385: /* g = 0 */
386: DEG(g) = -1;
387:
388: /* pwr = 1 */
389: pwr = W_UMALLOC(n); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
390:
391: /* xa = x+a */
392: xa = W_UMALLOC(n); DEG(xa) = 1;
393: COEF(xa)[0] = a; COEF(xa)[1] = _onesf();
394:
395: w = W_UMALLOC(n);
396: t = W_UMALLOC(n);
397:
398: /* g += pwr*c[0] */
399: for ( i = 0; i <= n; i++ ) {
400: mulssfum(pwr,c[i],w); addsfum(g,w,t); cpyum(t,g);
401: if ( i < n ) {
402: mulsfum(pwr,xa,t); cpyum(t,pwr);
403: }
404: }
405: }
406: }
407:
1.4 ! noro 408: /* f(y) -> f(y+a) */
! 409:
! 410: void shiftsfbm(bound,f,a)
! 411: int bound;
! 412: BM f;
! 413: int a;
! 414: {
! 415: int i,j,d;
! 416: UM pwr,ya,w,t,s;
! 417: UM *c;
! 418:
! 419: if ( bound <= 0 )
! 420: return;
! 421: else {
! 422: c = COEF(f);
! 423: d = DEG(c[0]);
! 424: for ( i = 1; i < bound; i++ )
! 425: d = MAX(DEG(c[i]),d);
! 426:
! 427: w = W_UMALLOC(d);
! 428: t = W_UMALLOC(d);
! 429: s = W_UMALLOC(bound);
! 430:
! 431: /* pwr = 1 */
! 432: pwr = W_UMALLOC(bound); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
! 433:
! 434: /* ya = y+a */
! 435: ya = W_UMALLOC(1); DEG(ya) = 1;
! 436: COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
! 437:
! 438: for ( i = 0; i < bound; i++ ) {
! 439: /* c[i] does not change */
! 440: for ( j = 0; j < i; j++ ) {
! 441: mulssfum(c[i],COEF(pwr)[j],w);
! 442: addsfum(w,c[j],t); cpyum(t,c[j]);
! 443: }
! 444: if ( i < bound-1 ) {
! 445: mulsfum(pwr,ya,s); cpyum(s,pwr);
! 446: }
! 447: }
! 448: }
! 449: }
! 450:
1.3 noro 451: /* clear the body of f */
452:
453: void clearsflum(bound,n,f)
454: int bound,n;
455: LUM f;
456: {
457: int **c;
458: int i;
459:
460: DEG(f) = 0;
461: for ( c = COEF(f), i = 0; i <= n; i++ )
462: bzero(c[i],(bound+1)*sizeof(int));
1.4 ! noro 463: }
! 464:
! 465: void clearsfbm(bound,n,f)
! 466: int bound,n;
! 467: BM f;
! 468: {
! 469: int i;
! 470: UM *c;
! 471:
! 472: DEG(f) = bound;
! 473: for ( c = COEF(f), i = 0; i < bound; i++ )
! 474: clearum(c[i],n);
1.3 noro 475: }
476:
477: /* f += g */
478:
479: void addtosflum(bound,g,f)
480: int bound;
481: LUM g,f;
482: {
483: int dg,i,j;
484: int **cg,**cf;
485:
486: dg = DEG(g);
487: cg = COEF(g);
488: cf = COEF(f);
489: for ( i = 0; i <= dg; i++ )
490: for ( j = 0; j <= bound; j++ )
491: cf[i][j] = _addsf(cf[i][j],cg[i][j]);
1.1 noro 492: }