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