![[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.12 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.11 2001/10/09 01:36:10 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.12 ! noro 540: }
! 541:
! 542: static void create_bmono(P c,V x,int i,V y,int j,P *mono);
! 543:
! 544: struct lb {
! 545: int pos,len;
! 546: int *r;
! 547: int *hist;
! 548: };
! 549:
! 550: static NODE insert_lb(NODE g,struct lb *a)
! 551: {
! 552: NODE prev,cur,n;
! 553:
! 554: prev = 0; cur = g;
! 555: while ( cur ) {
! 556: if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
! 557: MKNODE(n,a,cur);
! 558: if ( !prev )
! 559: return n;
! 560: else {
! 561: NEXT(prev) = n;
! 562: return g;
! 563: }
! 564: } else {
! 565: prev = cur;
! 566: cur = NEXT(cur);
! 567: }
! 568: }
! 569: MKNODE(n,a,0);
! 570: NEXT(prev) = n;
! 571: return g;
! 572: }
! 573:
! 574: static void lnf(int *r,int *h,int n,int len,NODE g)
! 575: {
! 576: struct lb *t;
! 577: int pos,i,j,len1,c;
! 578: int *r1,*h1;
! 579:
! 580: for ( ; g; g = NEXT(g) ) {
! 581: t = (struct lb *)BDY(g);
! 582: pos = t->pos;
! 583: if ( c = r[pos] ) {
! 584: c = _chsgnsf(c);
! 585: r1 = t->r;
! 586: h1 = t->hist;
! 587: len1 = t->len;
! 588: for ( i = pos; i < n; i++ )
! 589: if ( r1[i] )
! 590: r[i] = _ADDSF(r[i],_MULSF(r1[i],c));
! 591: for ( i = 0; i < len1; i++ )
! 592: if ( h1[i] )
! 593: h[i] = _ADDSF(h[i],_MULSF(h1[i],c));
! 594: }
! 595: }
! 596: for ( i = 0; i < n && !r[i]; i++ );
! 597: if ( i < n ) {
! 598: c = _invsf(r[i]);
! 599: for ( j = i; j < n; j++ )
! 600: if ( r[j] )
! 601: r[j] = _MULSF(r[j],c);
! 602: for ( j = 0; j < len; j++ )
! 603: if ( h[j] )
! 604: h[j] = _MULSF(h[j],c);
! 605: }
! 606: }
! 607:
! 608: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
! 609: {
! 610: V x,y;
! 611: int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
! 612: UP *rx;
! 613: DCP dc;
! 614: P t,f,mono,f1;
! 615: UP ut,h;
! 616: int ***nf;
! 617: int *r,*hist,*prev,*r1;
! 618: struct lb *lb;
! 619: GFS s;
! 620: NODE g;
! 621:
! 622: x = vl->v;
! 623: y = NEXT(vl)->v;
! 624: dx = getdeg(x,fx);
! 625: dxdy = dx*dy;
! 626: /* rx = -(fx-x^dx) */
! 627: rx = (UP *)CALLOC(dx,sizeof(UP));
! 628: for ( dc = DC(fx); dc; dc = NEXT(dc)) {
! 629: chsgnp(COEF(dc),&t);
! 630: ptoup(t,&ut);
! 631: rx[QTOS(DEG(dc))] = ut;
! 632: }
! 633: /* nf[d] = normal form table of monomials with total degree d */
! 634: nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
! 635: nf[0] = (int **)CALLOC(1,sizeof(int *));
! 636:
! 637: /* nf[0][0] = 1 */
! 638: r = (int *)CALLOC(dxdy,sizeof(int));
! 639: r[0] = _onesf();
! 640: nf[0][0] = r;
! 641:
! 642: hist = (int *)CALLOC(1,sizeof(int));
! 643: hist[0] = _onesf();
! 644:
! 645: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
! 646: lb->pos = 0;
! 647: lb->r = r;
! 648: lb->hist = hist;
! 649: lb->len = 1;
! 650:
! 651: /* g : table of normal form as linear form */
! 652: MKNODE(g,lb,0);
! 653:
! 654: len = 1;
! 655: h = UPALLOC(dy);
! 656: for ( d = 1; ; d++ ) {
! 657: if ( d > c ){
! 658: return;
! 659: }
! 660: nf[d] = (int **)CALLOC(d+1,sizeof(int *));
! 661: len0 = len;
! 662: len += d+1;
! 663:
! 664: for ( i = d; i >= 0; i-- ) {
! 665: /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
! 666: /* nf(x^d) = nf(nf(x^(d-1))*x) */
! 667: r = (int *)CALLOC(dxdy,sizeof(int));
! 668: if ( i == 0 ) {
! 669: prev = nf[d-1][0];
! 670: bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
! 671:
! 672: /* create the head coeff */
! 673: for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
! 674: if ( prev[k] ) {
! 675: u = IFTOF(prev[k]);
! 676: MKGFS(u,s);
! 677: } else
! 678: s = 0;
! 679: COEF(h)[l] = (Num)s;
! 680: }
! 681: for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
! 682: DEG(h) = l;
! 683:
! 684: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
! 685: tmulup(rx[k],h,dy,&ut);
! 686: if ( ut )
! 687: for ( l = 0; l < dy; l++ ) {
! 688: s = (GFS)COEF(ut)[l];
! 689: if ( s ) {
! 690: u = CONT(s);
! 691: r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));
! 692: }
! 693: }
! 694: }
! 695: } else {
! 696: prev = nf[d-1][i-1];
! 697: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
! 698: for ( l = 1; l < dy; l++ )
! 699: r[dyk+l] = prev[dyk+l-1];
! 700: }
! 701: }
! 702: nf[d][i] = r;
! 703: hist = (int *)CALLOC(len,sizeof(int));
! 704: hist[len0+i] = _onesf();
! 705: r1 = (int *)CALLOC(dxdy,sizeof(int));
! 706: bcopy(r,r1,dxdy*sizeof(int));
! 707: lnf(r1,hist,dxdy,len,g);
! 708: for ( k = 0; k < dxdy && !r1[k]; k++ );
! 709: if ( k == dxdy ) {
! 710: f = 0;
! 711: for ( k = j = 0; k <= d; k++ )
! 712: for ( i = 0; i <= k; i++, j++ )
! 713: if ( hist[j] ) {
! 714: u = IFTOF(hist[j]);
! 715: MKGFS(u,s);
! 716: /* mono = s*x^(k-i)*y^i */
! 717: create_bmono((P)s,x,k-i,y,i,&mono);
! 718: addp(vl,f,mono,&f1);
! 719: f = f1;
! 720: }
! 721: *fr = f;
! 722: return;
! 723: } else {
! 724: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
! 725: lb->pos = k;
! 726: lb->r = r1;
! 727: lb->hist = hist;
! 728: lb->len = len;
! 729: g = insert_lb(g,lb);
! 730: }
! 731: }
! 732: }
! 733: }
! 734:
! 735: static void create_bmono(P c,V x,int i,V y,int j,P *mono)
! 736: {
! 737: P t,s;
! 738:
! 739: if ( !i )
! 740: if ( !j )
! 741: t = c;
! 742: else {
! 743: /* c*y^j */
! 744: MKV(y,t);
! 745: COEF(DC(t)) = c;
! 746: STOQ(j,DEG(DC(t)));
! 747: }
! 748: else if ( !j ) {
! 749: /* c*x^i */
! 750: MKV(x,t);
! 751: COEF(DC(t)) = c;
! 752: STOQ(i,DEG(DC(t)));
! 753: } else {
! 754: MKV(y,s);
! 755: COEF(DC(s)) = c;
! 756: STOQ(j,DEG(DC(s)));
! 757: MKV(x,t);
! 758: COEF(DC(t)) = s;
! 759: STOQ(i,DEG(DC(t)));
! 760: }
! 761: *mono = t;
1.1 noro 762: }