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