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