![[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.15 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.14 2002/09/27 08:40:49 noro Exp $ */
1.1 noro 2:
3: #include "ca.h"
1.13 noro 4: #include "inline.h"
1.1 noro 5:
1.15 ! noro 6: extern int current_gfs_p;
1.5 noro 7: extern int up_kara_mag, current_gfs_q1;
8: extern int *current_gfs_plus1;
9:
1.9 noro 10: #if defined(__GNUC__)
11: #define INLINE inline
12: #elif defined(VISUAL)
13: #define INLINE __inline
14: #else
15: #define INLINE
16: #endif
17:
1.11 noro 18: INLINE int _ADDSF(int a,int b)
1.5 noro 19: {
20: if ( !a )
21: return b;
22: else if ( !b )
23: return a;
24:
25: a = IFTOF(a); b = IFTOF(b);
1.13 noro 26:
27: if ( !current_gfs_ntoi ) {
28: a = a+b-current_gfs_q;
29: if ( a == 0 )
30: return 0;
31: else {
32: if ( a < 0 )
33: a += current_gfs_q;
34: return FTOIF(a);
35: }
36: }
37:
1.5 noro 38: if ( a > b ) {
39: /* tab[a]+tab[b] = tab[b](tab[a-b]+1) */
40: a = current_gfs_plus1[a-b];
41: if ( a < 0 )
42: return 0;
43: else {
44: a += b;
45: if ( a >= current_gfs_q1 )
46: a -= current_gfs_q1;
47: return FTOIF(a);
48: }
49: } else {
50: /* tab[a]+tab[b] = tab[a](tab[b-a]+1) */
51: b = current_gfs_plus1[b-a];
52: if ( b < 0 )
53: return 0;
54: else {
55: b += a;
56: if ( b >= current_gfs_q1 )
57: b -= current_gfs_q1;
58: return FTOIF(b);
59: }
60: }
61: }
62:
1.11 noro 63: INLINE int _MULSF(int a,int b)
1.5 noro 64: {
1.13 noro 65: int c;
66:
1.5 noro 67: if ( !a || !b )
68: return 0;
1.13 noro 69: else if ( !current_gfs_ntoi ) {
70: a = IFTOF(a); b = IFTOF(b);
71: DMAR(a,b,0,current_gfs_q,c);
72: return FTOIF(c);
73: } else {
1.5 noro 74: a = IFTOF(a) + IFTOF(b);
75: if ( a >= current_gfs_q1 )
76: a -= current_gfs_q1;
77: return FTOIF(a);
78: }
79: }
1.1 noro 80:
1.11 noro 81: void addsfum(UM p1,UM p2,UM pr)
1.1 noro 82: {
83: int *c1,*c2,*cr,i,dmax,dmin;
84:
85: if ( DEG(p1) == -1 ) {
86: cpyum(p2,pr);
87: return;
88: }
89: if ( DEG(p2) == -1 ) {
90: cpyum(p1,pr);
91: return;
92: }
93: if ( DEG(p1) >= DEG(p2) ) {
94: c1 = COEF(p1); c2 = COEF(p2); dmax = DEG(p1); dmin = DEG(p2);
95: } else {
96: c1 = COEF(p2); c2 = COEF(p1); dmax = DEG(p2); dmin = DEG(p1);
97: }
98: for ( i = 0, cr = COEF(pr); i <= dmin; i++ )
1.5 noro 99: cr[i] = _ADDSF(c1[i],c2[i]);
1.1 noro 100: for ( ; i <= dmax; i++ )
101: cr[i] = c1[i];
102: if ( dmax == dmin )
103: degum(pr,dmax);
104: else
105: DEG(pr) = dmax;
106: }
107:
1.11 noro 108: void subsfum(UM p1,UM p2,UM pr)
1.1 noro 109: {
110: int *c1,*c2,*cr,i;
111: int dmax,dmin;
112:
113: if ( DEG(p1) == -1 ) {
114: for ( i = DEG(pr) = DEG(p2), c2 = COEF(p2), cr = COEF(pr);
115: i >= 0; i-- )
116: cr[i] = _chsgnsf(c2[i]);
117: return;
118: }
119: if ( DEG(p2) == -1 ) {
120: cpyum(p1,pr);
121: return;
122: }
123: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
124: if ( DEG(p1) >= DEG(p2) ) {
125: dmax = DEG(p1); dmin = DEG(p2);
126: for ( i = 0; i <= dmin; i++ )
127: cr[i] = _subsf(c1[i],c2[i]);
128: for ( ; i <= dmax; i++ )
129: cr[i] = c1[i];
130: } else {
131: dmax = DEG(p2); dmin = DEG(p1);
132: for ( i = 0; i <= dmin; i++ )
133: cr[i] = _subsf(c1[i],c2[i]);
134: for ( ; i <= dmax; i++ )
135: cr[i] = _chsgnsf(c2[i]);
136: }
137: if ( dmax == dmin )
138: degum(pr,dmax);
139: else
140: DEG(pr) = dmax;
141: }
142:
1.11 noro 143: void gcdsfum(UM p1,UM p2,UM pr)
1.1 noro 144: {
145: int inv;
146: UM t1,t2,q,tum;
147: int drem;
148:
149: if ( DEG(p1) < 0 )
150: cpyum(p2,pr);
151: else if ( DEG(p2) < 0 )
152: cpyum(p1,pr);
153: else {
154: if ( DEG(p1) >= DEG(p2) ) {
155: t1 = p1; t2 = p2;
156: } else {
157: t1 = p2; t2 = p1;
158: }
159: q = W_UMALLOC(DEG(t1));
160: while ( ( drem = divsfum(t1,t2,q) ) >= 0 ) {
161: tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;
162: }
163: inv = _invsf(COEF(t2)[DEG(t2)]);
164: mulssfum(t2,inv,pr);
165: }
166: }
1.11 noro 167:
168: void mulsfum(UM p1,UM p2,UM pr)
1.1 noro 169: {
170: int *pc1,*pcr;
171: int *c1,*c2,*cr;
172: int mul;
173: int i,j,d1,d2;
174:
175: if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
176: DEG(pr) = -1;
177: return;
178: }
179: c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
180: bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
181: for ( i = 0; i <= d2; i++, cr++ )
182: if ( mul = *c2++ )
183: for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )
1.5 noro 184: if ( *pc1 )
185: *pcr = _ADDSF(_MULSF(*pc1,mul),*pcr);
1.1 noro 186: DEG(pr) = d1 + d2;
187: }
188:
1.11 noro 189: void mulssfum(UM p,int n,UM pr)
1.1 noro 190: {
191: int *sp,*dp;
192: int i;
193:
194: for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
195: i >= 0; i--, dp--, sp-- )
1.5 noro 196: *dp = _MULSF(*sp,n);
1.1 noro 197: }
198:
1.11 noro 199: void kmulsfum(UM n1,UM n2,UM nr)
1.5 noro 200: {
201: UM n,t,s,m,carry;
202: int d,d1,d2,len,i,l;
203: unsigned int *r,*r0;
204:
205: if ( !n1 || !n2 ) {
206: nr->d = -1; return;
207: }
208: d1 = DEG(n1)+1; d2 = DEG(n2)+1;
209: if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {
210: mulsfum(n1,n2,nr); return;
211: }
212: if ( d1 < d2 ) {
213: n = n1; n1 = n2; n2 = n;
214: d = d1; d1 = d2; d2 = d;
215: }
216: if ( d2 > (d1+1)/2 ) {
217: kmulsfummain(n1,n2,nr); return;
218: }
219: d = (d1/d2)+((d1%d2)!=0);
220: len = (d+1)*d2;
221: r0 = (unsigned int *)ALLOCA(len*sizeof(int));
222: bzero((char *)r0,len*sizeof(int));
223: m = W_UMALLOC(d2+1);
224: carry = W_UMALLOC(d2+1);
225: t = W_UMALLOC(d1+d2+1);
226: s = W_UMALLOC(d1+d2+1);
227: for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {
228: extractum(n1,i*d2,d2,m);
229: if ( m ) {
230: kmulsfum(m,n2,t);
231: addsfum(t,carry,s);
232: c_copyum(s,d2,r);
233: extractum(s,d2,d2,carry);
234: } else {
235: c_copyum(carry,d2,r);
236: carry = 0;
237: }
238: }
239: c_copyum(carry,d2,r);
240: for ( l = len - 1; !r0[l]; l-- );
241: l++;
242: DEG(nr) = l-1;
243: bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
244: }
245:
1.11 noro 246: void kmulsfummain(UM n1,UM n2,UM nr)
1.5 noro 247: {
248: int d1,d2,h,len;
249: UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
250:
251: d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
252: n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);
253: n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);
254: lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);
255: mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);
256: mid = W_UMALLOC(d1+d2+1);
257: s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);
258: extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
259: extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
260: kmulsfum(n1hi,n2hi,hi); kmulsfum(n1lo,n2lo,lo);
261: len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
262: bzero((char *)COEF(t1),len*sizeof(int));
263: if ( lo )
264: bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
265: if ( hi )
266: bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
267:
268: addsfum(hi,lo,mid1);
269: subsfum(n1hi,n1lo,s1); subsfum(n2lo,n2hi,s2);
270: kmulsfum(s1,s2,mid2); addsfum(mid1,mid2,mid);
271: if ( mid ) {
272: len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
273: bzero((char *)COEF(t2),len*sizeof(int));
274: bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
275: addsfum(t1,t2,nr);
276: } else
277: copyum(t1,nr);
278: }
279:
1.11 noro 280: int divsfum(UM p1,UM p2,UM pq)
1.1 noro 281: {
282: int *pc1,*pct;
283: int *c1,*c2,*ct;
284: int inv,hd,tmp;
285: int i,j, d1,d2,dd;
286:
287: if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
288: DEG(pq) = -1;
289: return d1;
290: }
291: c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
292: if ( ( hd = c2[d2] ) != _onesf() ) {
293: inv = _invsf(hd);
294: for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
1.5 noro 295: *pc1 = _MULSF(*pc1,inv);
1.1 noro 296: } else
297: inv = _onesf();
298: for ( i = dd, ct = c1+d1; i >= 0; i-- )
299: if ( tmp = *ct-- ) {
300: tmp = _chsgnsf(tmp);
301: for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
1.5 noro 302: *pct = _ADDSF(_MULSF(*pc1,tmp),*pct);
1.1 noro 303: }
304: if ( inv != _onesf() ) {
305: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
1.5 noro 306: *pc1 = _MULSF(*pc1,inv);
1.1 noro 307: for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )
1.5 noro 308: *pc1 = _MULSF(*pc1,hd);
1.1 noro 309: }
310: for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
311: for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
312: *pct-- = *pc1--;
313: return i;
314: }
315:
1.11 noro 316: void diffsfum(UM f,UM fd)
1.1 noro 317: {
318: int *dp,*sp;
319: int i;
320:
321: for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
322: i >= 1; i--, dp--, sp-- ) {
1.15 ! noro 323: *dp = _MULSF(*sp,_itosf(i%current_gfs_p));
1.1 noro 324: }
325: degum(fd,DEG(f) - 1);
326: }
327:
1.11 noro 328: void monicsfum(UM f)
1.1 noro 329: {
330: int *sp;
331: int i,inv;
332:
333: i = DEG(f); sp = COEF(f)+i;
334: inv = _invsf(*sp);
335: for ( ; i >= 0; i--, sp-- )
1.5 noro 336: *sp = _MULSF(*sp,inv);
1.10 noro 337: }
338:
1.11 noro 339: int compsfum(UM a,UM b)
1.10 noro 340: {
341: int i,da,db;
342:
343: if ( !a )
344: return !b?0:1;
345: else if ( !b )
346: return 1;
347: else if ( (da = DEG(a)) > (db = DEG(b)) )
348: return 1;
349: else if ( da < db )
350: return -1;
351: else {
352: for ( i = da; i >= 0 && COEF(a)[i] == COEF(b)[i]; i-- );
353: if ( i < 0 )
354: return 0;
355: else if ( (unsigned int)COEF(a)[i] > (unsigned int)COEF(b)[i] )
356: return 1;
357: else
358: return -1;
359: }
1.2 noro 360: }
361:
1.3 noro 362: void addsfarray(int,int *,int *);
363: void mulsfarray_trunc(int,int *,int *,int *);
1.2 noro 364:
1.4 noro 365: /* f1 = f1->c[0]+f1->c[1]*y+..., f2 = f2->c[0]+f2->c[1]*y+... mod y^n */
366:
1.11 noro 367: void mulsfbm(BM f1,BM f2,BM fr)
1.4 noro 368: {
369: UM mul,t,s;
1.11 noro 370: int i,j,d1,d2,dy;
1.4 noro 371:
1.8 noro 372: dy = MIN(DEG(f1),DEG(f2));
1.4 noro 373:
1.8 noro 374: d1 = degbm(f1);
375: d2 = degbm(f2);
1.4 noro 376: t = W_UMALLOC(d1+d2);
377: s = W_UMALLOC(d1+d2);
378:
1.8 noro 379: for ( i = 0; i <= dy; i++ ) {
1.4 noro 380: mul = COEF(f2)[i];
381: if ( DEG(mul) >= 0 )
1.8 noro 382: for ( j = 0; i+j <= dy; j++ ) {
1.4 noro 383: if ( COEF(f1)[j] ) {
1.7 noro 384: kmulsfum(COEF(f1)[j],mul,t);
1.4 noro 385: addsfum(t,COEF(fr)[i+j],s);
386: cpyum(s,COEF(fr)[i+j]);
387: }
388: }
389: }
1.8 noro 390: DEG(fr) = dy;
1.4 noro 391: }
392:
1.11 noro 393: int degbm(BM f)
1.6 noro 394: {
1.8 noro 395: int d,i,dy;
1.6 noro 396:
1.8 noro 397: dy = DEG(f);
1.6 noro 398: d = DEG(COEF(f)[0]);
1.8 noro 399: for ( i = 1; i <= dy; i++ )
1.6 noro 400: d = MAX(DEG(COEF(f)[i]),d);
401: return d;
402: }
403:
1.4 noro 404: /* g += f */
405:
1.11 noro 406: void addtosfbm(BM f,BM g)
1.4 noro 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.11 noro 424: void eucsfum(UM f1,UM f2,UM a,UM b)
1.2 noro 425: {
426: UM g1,g2,a1,a2,a3,wm,q,tum;
427: int d,dr;
428:
429: /* XXX */
430: d = DEG(f1) + DEG(f2) + 10;
431: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);
432: a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);
433: q = W_UMALLOC(d);
1.3 noro 434: DEG(a1) = 0; COEF(a1)[0] = _onesf(); DEG(a2) = -1;
1.2 noro 435: cpyum(f1,g1); cpyum(f2,g2);
436: while ( 1 ) {
437: dr = divsfum(g1,g2,q); tum = g1; g1 = g2; g2 = tum;
438: if ( ( DEG(g2) = dr ) == -1 )
439: break;
440: mulsfum(a2,q,wm); subsfum(a1,wm,a3); dr = divsfum(a3,f2,q);
441: tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;
442: }
1.3 noro 443: if ( COEF(g1)[0] != _onesf() )
1.2 noro 444: mulssfum(a2,_invsf(COEF(g1)[0]),a);
445: else
446: cpyum(a2,a);
447: mulsfum(a,f1,wm);
448: if ( DEG(wm) >= 0 )
449: COEF(wm)[0] = _subsf(COEF(wm)[0],_onesf());
450: else {
451: DEG(wm) = 0; COEF(wm)[0] = _chsgnsf(_onesf());
452: }
453: divsfum(wm,f2,q); mulssfum(q,_chsgnsf(_onesf()),b);
1.3 noro 454: }
455:
456: void shiftsfum(UM,int,UM);
457:
1.11 noro 458: void shiftsflum(int n,LUM f,int ev)
1.3 noro 459: {
460: int d,i,j;
461: UM w,w1;
462: int *coef;
463: int **c;
464:
465: d = f->d;
466: c = f->c;
467: w = W_UMALLOC(n);
468: w1 = W_UMALLOC(n);
469: for ( i = 0; i <= d; i++ ) {
470: coef = c[i];
471: for ( j = n-1; j >= 0 && coef[j] == 0; j-- );
472: if ( j <= 0 )
473: continue;
474: DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));
475: shiftsfum(w,ev,w1);
476: bcopy(COEF(w1),coef,(j+1)*sizeof(int));
477: }
478: }
479:
480: /* f(x) -> g(x) = f(x+a) */
481:
1.11 noro 482: void shiftsfum(UM f,int a,UM g)
1.3 noro 483: {
484: int n,i;
485: UM pwr,xa,w,t;
486: int *c;
487:
488: n = DEG(f);
489: if ( n <= 0 )
490: cpyum(f,g);
491: else {
492: c = COEF(f);
493:
494: /* g = 0 */
495: DEG(g) = -1;
496:
497: /* pwr = 1 */
498: pwr = W_UMALLOC(n); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
499:
500: /* xa = x+a */
501: xa = W_UMALLOC(n); DEG(xa) = 1;
502: COEF(xa)[0] = a; COEF(xa)[1] = _onesf();
503:
504: w = W_UMALLOC(n);
505: t = W_UMALLOC(n);
506:
507: /* g += pwr*c[0] */
508: for ( i = 0; i <= n; i++ ) {
509: mulssfum(pwr,c[i],w); addsfum(g,w,t); cpyum(t,g);
510: if ( i < n ) {
511: mulsfum(pwr,xa,t); cpyum(t,pwr);
512: }
513: }
514: }
515: }
516:
1.4 noro 517: /* f(y) -> f(y+a) */
518:
1.11 noro 519: void shiftsfbm(BM f,int a)
1.4 noro 520: {
1.8 noro 521: int i,j,d,dy;
1.4 noro 522: UM pwr,ya,w,t,s;
523: UM *c;
524:
1.8 noro 525: dy = DEG(f);
526: c = COEF(f);
527: d = degbm(f);
528:
529: w = W_UMALLOC(d);
530: t = W_UMALLOC(d);
531: s = W_UMALLOC(dy);
532:
533: /* pwr = 1 */
534: pwr = W_UMALLOC(dy); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
535:
536: /* ya = y+a */
537: ya = W_UMALLOC(1); DEG(ya) = 1;
538: COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
539:
540: for ( i = 0; i <= dy; i++ ) {
541: /* c[i] does not change */
542: for ( j = 0; j < i; j++ ) {
543: mulssfum(c[i],COEF(pwr)[j],w);
544: addsfum(w,c[j],t); cpyum(t,c[j]);
545: }
546: if ( i < dy ) {
547: mulsfum(pwr,ya,s); cpyum(s,pwr);
1.4 noro 548: }
549: }
550: }
551:
1.11 noro 552: void clearbm(int n,BM f)
1.4 noro 553: {
1.8 noro 554: int i,dy;
1.4 noro 555: UM *c;
556:
1.8 noro 557: dy = DEG(f);
558: for ( c = COEF(f), i = 0; i <= dy; i++ )
1.4 noro 559: clearum(c[i],n);
1.12 noro 560: }
561:
562: static void create_bmono(P c,V x,int i,V y,int j,P *mono);
563:
564: struct lb {
565: int pos,len;
566: int *r;
567: int *hist;
568: };
569:
570: static NODE insert_lb(NODE g,struct lb *a)
571: {
572: NODE prev,cur,n;
573:
574: prev = 0; cur = g;
575: while ( cur ) {
576: if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
577: MKNODE(n,a,cur);
578: if ( !prev )
579: return n;
580: else {
581: NEXT(prev) = n;
582: return g;
583: }
584: } else {
585: prev = cur;
586: cur = NEXT(cur);
587: }
588: }
589: MKNODE(n,a,0);
590: NEXT(prev) = n;
591: return g;
592: }
593:
594: static void lnf(int *r,int *h,int n,int len,NODE g)
595: {
596: struct lb *t;
597: int pos,i,j,len1,c;
598: int *r1,*h1;
599:
600: for ( ; g; g = NEXT(g) ) {
601: t = (struct lb *)BDY(g);
602: pos = t->pos;
603: if ( c = r[pos] ) {
604: c = _chsgnsf(c);
605: r1 = t->r;
606: h1 = t->hist;
607: len1 = t->len;
608: for ( i = pos; i < n; i++ )
609: if ( r1[i] )
610: r[i] = _ADDSF(r[i],_MULSF(r1[i],c));
611: for ( i = 0; i < len1; i++ )
612: if ( h1[i] )
613: h[i] = _ADDSF(h[i],_MULSF(h1[i],c));
614: }
615: }
616: for ( i = 0; i < n && !r[i]; i++ );
617: if ( i < n ) {
618: c = _invsf(r[i]);
619: for ( j = i; j < n; j++ )
620: if ( r[j] )
621: r[j] = _MULSF(r[j],c);
622: for ( j = 0; j < len; j++ )
623: if ( h[j] )
624: h[j] = _MULSF(h[j],c);
625: }
626: }
627:
628: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
629: {
630: V x,y;
631: int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
632: UP *rx;
633: DCP dc;
634: P t,f,mono,f1;
635: UP ut,h;
636: int ***nf;
637: int *r,*hist,*prev,*r1;
638: struct lb *lb;
639: GFS s;
640: NODE g;
641:
642: x = vl->v;
643: y = NEXT(vl)->v;
644: dx = getdeg(x,fx);
645: dxdy = dx*dy;
646: /* rx = -(fx-x^dx) */
647: rx = (UP *)CALLOC(dx,sizeof(UP));
648: for ( dc = DC(fx); dc; dc = NEXT(dc)) {
649: chsgnp(COEF(dc),&t);
650: ptoup(t,&ut);
651: rx[QTOS(DEG(dc))] = ut;
652: }
653: /* nf[d] = normal form table of monomials with total degree d */
654: nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
655: nf[0] = (int **)CALLOC(1,sizeof(int *));
656:
657: /* nf[0][0] = 1 */
658: r = (int *)CALLOC(dxdy,sizeof(int));
659: r[0] = _onesf();
660: nf[0][0] = r;
661:
662: hist = (int *)CALLOC(1,sizeof(int));
663: hist[0] = _onesf();
664:
665: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
666: lb->pos = 0;
667: lb->r = r;
668: lb->hist = hist;
669: lb->len = 1;
670:
671: /* g : table of normal form as linear form */
672: MKNODE(g,lb,0);
673:
674: len = 1;
675: h = UPALLOC(dy);
676: for ( d = 1; ; d++ ) {
677: if ( d > c ){
678: return;
679: }
680: nf[d] = (int **)CALLOC(d+1,sizeof(int *));
681: len0 = len;
682: len += d+1;
683:
684: for ( i = d; i >= 0; i-- ) {
685: /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
686: /* nf(x^d) = nf(nf(x^(d-1))*x) */
687: r = (int *)CALLOC(dxdy,sizeof(int));
688: if ( i == 0 ) {
689: prev = nf[d-1][0];
690: bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
691:
692: /* create the head coeff */
693: for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
1.14 noro 694: iftogfs(prev[k],&s);
1.12 noro 695: COEF(h)[l] = (Num)s;
696: }
697: for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
698: DEG(h) = l;
699:
700: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
701: tmulup(rx[k],h,dy,&ut);
702: if ( ut )
703: for ( l = 0; l < dy; l++ ) {
704: s = (GFS)COEF(ut)[l];
705: if ( s ) {
706: u = CONT(s);
707: r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));
708: }
709: }
710: }
711: } else {
712: prev = nf[d-1][i-1];
713: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
714: for ( l = 1; l < dy; l++ )
715: r[dyk+l] = prev[dyk+l-1];
716: }
717: }
718: nf[d][i] = r;
719: hist = (int *)CALLOC(len,sizeof(int));
720: hist[len0+i] = _onesf();
721: r1 = (int *)CALLOC(dxdy,sizeof(int));
722: bcopy(r,r1,dxdy*sizeof(int));
723: lnf(r1,hist,dxdy,len,g);
724: for ( k = 0; k < dxdy && !r1[k]; k++ );
725: if ( k == dxdy ) {
726: f = 0;
727: for ( k = j = 0; k <= d; k++ )
728: for ( i = 0; i <= k; i++, j++ )
729: if ( hist[j] ) {
1.14 noro 730: iftogfs(hist[j],&s);
1.12 noro 731: /* mono = s*x^(k-i)*y^i */
732: create_bmono((P)s,x,k-i,y,i,&mono);
733: addp(vl,f,mono,&f1);
734: f = f1;
735: }
736: *fr = f;
737: return;
738: } else {
739: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
740: lb->pos = k;
741: lb->r = r1;
742: lb->hist = hist;
743: lb->len = len;
744: g = insert_lb(g,lb);
745: }
746: }
747: }
748: }
749:
750: static void create_bmono(P c,V x,int i,V y,int j,P *mono)
751: {
752: P t,s;
753:
754: if ( !i )
755: if ( !j )
756: t = c;
757: else {
758: /* c*y^j */
759: MKV(y,t);
760: COEF(DC(t)) = c;
761: STOQ(j,DEG(DC(t)));
762: }
763: else if ( !j ) {
764: /* c*x^i */
765: MKV(x,t);
766: COEF(DC(t)) = c;
767: STOQ(i,DEG(DC(t)));
768: } else {
769: MKV(y,s);
770: COEF(DC(s)) = c;
771: STOQ(j,DEG(DC(s)));
772: MKV(x,t);
773: COEF(DC(t)) = s;
774: STOQ(i,DEG(DC(t)));
775: }
776: *mono = t;
1.1 noro 777: }