![[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.19 ! fujimoto 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.18 2015/08/06 10:01:52 fujimoto 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__)
1.17 ohara 11: #define INLINE static inline
1.19 ! fujimoto 12: #elif defined(VISUAL) || defined(__MINGW32__)
1.9 noro 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));
1.16 noro 223: m = W_UMALLOC(d1+d2+1);
224: carry = W_UMALLOC(d1+d2+1);
1.5 noro 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: {
1.16 noro 248: int d1,d2,h,len,d;
1.5 noro 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;
1.16 noro 252: d = d1+d2+1;
253: n1lo = W_UMALLOC(d); n1hi = W_UMALLOC(d);
254: n2lo = W_UMALLOC(d); n2hi = W_UMALLOC(d);
255: lo = W_UMALLOC(d); hi = W_UMALLOC(d);
256: mid1 = W_UMALLOC(d); mid2 = W_UMALLOC(d);
257: mid = W_UMALLOC(d);
258: s1 = W_UMALLOC(d); s2 = W_UMALLOC(d);
1.5 noro 259: extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
260: extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
261: kmulsfum(n1hi,n2hi,hi); kmulsfum(n1lo,n2lo,lo);
262: len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
263: bzero((char *)COEF(t1),len*sizeof(int));
264: if ( lo )
265: bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
266: if ( hi )
267: bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
268:
269: addsfum(hi,lo,mid1);
270: subsfum(n1hi,n1lo,s1); subsfum(n2lo,n2hi,s2);
271: kmulsfum(s1,s2,mid2); addsfum(mid1,mid2,mid);
272: if ( mid ) {
273: len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
274: bzero((char *)COEF(t2),len*sizeof(int));
275: bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
276: addsfum(t1,t2,nr);
277: } else
278: copyum(t1,nr);
279: }
280:
1.11 noro 281: int divsfum(UM p1,UM p2,UM pq)
1.1 noro 282: {
283: int *pc1,*pct;
284: int *c1,*c2,*ct;
285: int inv,hd,tmp;
286: int i,j, d1,d2,dd;
287:
288: if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
289: DEG(pq) = -1;
290: return d1;
291: }
292: c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
293: if ( ( hd = c2[d2] ) != _onesf() ) {
294: inv = _invsf(hd);
295: for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
1.5 noro 296: *pc1 = _MULSF(*pc1,inv);
1.1 noro 297: } else
298: inv = _onesf();
299: for ( i = dd, ct = c1+d1; i >= 0; i-- )
300: if ( tmp = *ct-- ) {
301: tmp = _chsgnsf(tmp);
302: for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
1.5 noro 303: *pct = _ADDSF(_MULSF(*pc1,tmp),*pct);
1.1 noro 304: }
305: if ( inv != _onesf() ) {
306: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
1.5 noro 307: *pc1 = _MULSF(*pc1,inv);
1.1 noro 308: for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )
1.5 noro 309: *pc1 = _MULSF(*pc1,hd);
1.1 noro 310: }
311: for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
312: for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
313: *pct-- = *pc1--;
314: return i;
315: }
316:
1.11 noro 317: void diffsfum(UM f,UM fd)
1.1 noro 318: {
319: int *dp,*sp;
320: int i;
321:
322: for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
323: i >= 1; i--, dp--, sp-- ) {
1.15 noro 324: *dp = _MULSF(*sp,_itosf(i%current_gfs_p));
1.1 noro 325: }
326: degum(fd,DEG(f) - 1);
327: }
328:
1.11 noro 329: void monicsfum(UM f)
1.1 noro 330: {
331: int *sp;
332: int i,inv;
333:
334: i = DEG(f); sp = COEF(f)+i;
335: inv = _invsf(*sp);
336: for ( ; i >= 0; i--, sp-- )
1.5 noro 337: *sp = _MULSF(*sp,inv);
1.10 noro 338: }
339:
1.11 noro 340: int compsfum(UM a,UM b)
1.10 noro 341: {
342: int i,da,db;
343:
344: if ( !a )
345: return !b?0:1;
346: else if ( !b )
347: return 1;
348: else if ( (da = DEG(a)) > (db = DEG(b)) )
349: return 1;
350: else if ( da < db )
351: return -1;
352: else {
353: for ( i = da; i >= 0 && COEF(a)[i] == COEF(b)[i]; i-- );
354: if ( i < 0 )
355: return 0;
356: else if ( (unsigned int)COEF(a)[i] > (unsigned int)COEF(b)[i] )
357: return 1;
358: else
359: return -1;
360: }
1.2 noro 361: }
362:
1.3 noro 363: void addsfarray(int,int *,int *);
364: void mulsfarray_trunc(int,int *,int *,int *);
1.2 noro 365:
1.4 noro 366: /* f1 = f1->c[0]+f1->c[1]*y+..., f2 = f2->c[0]+f2->c[1]*y+... mod y^n */
367:
1.11 noro 368: void mulsfbm(BM f1,BM f2,BM fr)
1.4 noro 369: {
370: UM mul,t,s;
1.11 noro 371: int i,j,d1,d2,dy;
1.4 noro 372:
1.8 noro 373: dy = MIN(DEG(f1),DEG(f2));
1.4 noro 374:
1.8 noro 375: d1 = degbm(f1);
376: d2 = degbm(f2);
1.4 noro 377: t = W_UMALLOC(d1+d2);
378: s = W_UMALLOC(d1+d2);
379:
1.8 noro 380: for ( i = 0; i <= dy; i++ ) {
1.4 noro 381: mul = COEF(f2)[i];
382: if ( DEG(mul) >= 0 )
1.8 noro 383: for ( j = 0; i+j <= dy; j++ ) {
1.4 noro 384: if ( COEF(f1)[j] ) {
1.7 noro 385: kmulsfum(COEF(f1)[j],mul,t);
1.4 noro 386: addsfum(t,COEF(fr)[i+j],s);
387: cpyum(s,COEF(fr)[i+j]);
388: }
389: }
390: }
1.8 noro 391: DEG(fr) = dy;
1.4 noro 392: }
393:
1.11 noro 394: int degbm(BM f)
1.6 noro 395: {
1.8 noro 396: int d,i,dy;
1.6 noro 397:
1.8 noro 398: dy = DEG(f);
1.6 noro 399: d = DEG(COEF(f)[0]);
1.8 noro 400: for ( i = 1; i <= dy; i++ )
1.6 noro 401: d = MAX(DEG(COEF(f)[i]),d);
402: return d;
403: }
404:
1.4 noro 405: /* g += f */
406:
1.11 noro 407: void addtosfbm(BM f,BM g)
1.4 noro 408: {
1.8 noro 409: int i,d1,d2,dy;
1.4 noro 410: UM t;
411:
1.8 noro 412: dy = DEG(g);
413: if ( DEG(f) > dy )
414: error("addtosfbm : invalid input");
415: dy = MIN(DEG(f),dy);
416: d1 = degbm(f);
417: d2 = degbm(g);
1.4 noro 418: t = W_UMALLOC(MAX(d1,d2));
1.8 noro 419: for ( i = 0; i <= dy; i++ ) {
1.4 noro 420: addsfum(COEF(f)[i],COEF(g)[i],t);
421: cpyum(t,COEF(g)[i]);
422: }
423: }
424:
1.11 noro 425: void eucsfum(UM f1,UM f2,UM a,UM b)
1.2 noro 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:
1.11 noro 459: void shiftsflum(int n,LUM f,int ev)
1.3 noro 460: {
461: int d,i,j;
462: UM w,w1;
463: int *coef;
464: int **c;
465:
466: d = f->d;
467: c = f->c;
468: w = W_UMALLOC(n);
469: w1 = W_UMALLOC(n);
470: for ( i = 0; i <= d; i++ ) {
471: coef = c[i];
472: for ( j = n-1; j >= 0 && coef[j] == 0; j-- );
473: if ( j <= 0 )
474: continue;
475: DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));
476: shiftsfum(w,ev,w1);
477: bcopy(COEF(w1),coef,(j+1)*sizeof(int));
478: }
479: }
480:
481: /* f(x) -> g(x) = f(x+a) */
482:
1.11 noro 483: void shiftsfum(UM f,int a,UM g)
1.3 noro 484: {
485: int n,i;
486: UM pwr,xa,w,t;
487: int *c;
488:
489: n = DEG(f);
490: if ( n <= 0 )
491: cpyum(f,g);
492: else {
493: c = COEF(f);
494:
495: /* g = 0 */
496: DEG(g) = -1;
497:
498: /* pwr = 1 */
499: pwr = W_UMALLOC(n); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
500:
501: /* xa = x+a */
502: xa = W_UMALLOC(n); DEG(xa) = 1;
503: COEF(xa)[0] = a; COEF(xa)[1] = _onesf();
504:
505: w = W_UMALLOC(n);
506: t = W_UMALLOC(n);
507:
508: /* g += pwr*c[0] */
509: for ( i = 0; i <= n; i++ ) {
510: mulssfum(pwr,c[i],w); addsfum(g,w,t); cpyum(t,g);
511: if ( i < n ) {
512: mulsfum(pwr,xa,t); cpyum(t,pwr);
513: }
514: }
515: }
516: }
517:
1.4 noro 518: /* f(y) -> f(y+a) */
519:
1.11 noro 520: void shiftsfbm(BM f,int a)
1.4 noro 521: {
1.8 noro 522: int i,j,d,dy;
1.4 noro 523: UM pwr,ya,w,t,s;
524: UM *c;
525:
1.8 noro 526: dy = DEG(f);
527: c = COEF(f);
528: d = degbm(f);
529:
530: w = W_UMALLOC(d);
531: t = W_UMALLOC(d);
532: s = W_UMALLOC(dy);
533:
534: /* pwr = 1 */
535: pwr = W_UMALLOC(dy); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
536:
537: /* ya = y+a */
538: ya = W_UMALLOC(1); DEG(ya) = 1;
539: COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
540:
541: for ( i = 0; i <= dy; i++ ) {
542: /* c[i] does not change */
543: for ( j = 0; j < i; j++ ) {
544: mulssfum(c[i],COEF(pwr)[j],w);
545: addsfum(w,c[j],t); cpyum(t,c[j]);
546: }
547: if ( i < dy ) {
548: mulsfum(pwr,ya,s); cpyum(s,pwr);
1.4 noro 549: }
550: }
551: }
552:
1.11 noro 553: void clearbm(int n,BM f)
1.4 noro 554: {
1.8 noro 555: int i,dy;
1.4 noro 556: UM *c;
557:
1.8 noro 558: dy = DEG(f);
559: for ( c = COEF(f), i = 0; i <= dy; i++ )
1.4 noro 560: clearum(c[i],n);
1.12 noro 561: }
562:
563: static void create_bmono(P c,V x,int i,V y,int j,P *mono);
564:
565: struct lb {
566: int pos,len;
567: int *r;
568: int *hist;
569: };
570:
571: static NODE insert_lb(NODE g,struct lb *a)
572: {
573: NODE prev,cur,n;
574:
575: prev = 0; cur = g;
576: while ( cur ) {
577: if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
578: MKNODE(n,a,cur);
579: if ( !prev )
580: return n;
581: else {
582: NEXT(prev) = n;
583: return g;
584: }
585: } else {
586: prev = cur;
587: cur = NEXT(cur);
588: }
589: }
590: MKNODE(n,a,0);
591: NEXT(prev) = n;
592: return g;
593: }
594:
595: static void lnf(int *r,int *h,int n,int len,NODE g)
596: {
597: struct lb *t;
598: int pos,i,j,len1,c;
599: int *r1,*h1;
600:
601: for ( ; g; g = NEXT(g) ) {
602: t = (struct lb *)BDY(g);
603: pos = t->pos;
604: if ( c = r[pos] ) {
605: c = _chsgnsf(c);
606: r1 = t->r;
607: h1 = t->hist;
608: len1 = t->len;
609: for ( i = pos; i < n; i++ )
610: if ( r1[i] )
611: r[i] = _ADDSF(r[i],_MULSF(r1[i],c));
612: for ( i = 0; i < len1; i++ )
613: if ( h1[i] )
614: h[i] = _ADDSF(h[i],_MULSF(h1[i],c));
615: }
616: }
617: for ( i = 0; i < n && !r[i]; i++ );
618: if ( i < n ) {
619: c = _invsf(r[i]);
620: for ( j = i; j < n; j++ )
621: if ( r[j] )
622: r[j] = _MULSF(r[j],c);
623: for ( j = 0; j < len; j++ )
624: if ( h[j] )
625: h[j] = _MULSF(h[j],c);
626: }
627: }
628:
629: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
630: {
631: V x,y;
632: int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
633: UP *rx;
634: DCP dc;
635: P t,f,mono,f1;
636: UP ut,h;
637: int ***nf;
638: int *r,*hist,*prev,*r1;
639: struct lb *lb;
640: GFS s;
641: NODE g;
642:
643: x = vl->v;
644: y = NEXT(vl)->v;
645: dx = getdeg(x,fx);
646: dxdy = dx*dy;
647: /* rx = -(fx-x^dx) */
648: rx = (UP *)CALLOC(dx,sizeof(UP));
649: for ( dc = DC(fx); dc; dc = NEXT(dc)) {
650: chsgnp(COEF(dc),&t);
651: ptoup(t,&ut);
652: rx[QTOS(DEG(dc))] = ut;
653: }
654: /* nf[d] = normal form table of monomials with total degree d */
655: nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
656: nf[0] = (int **)CALLOC(1,sizeof(int *));
657:
658: /* nf[0][0] = 1 */
659: r = (int *)CALLOC(dxdy,sizeof(int));
660: r[0] = _onesf();
661: nf[0][0] = r;
662:
663: hist = (int *)CALLOC(1,sizeof(int));
664: hist[0] = _onesf();
665:
666: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
667: lb->pos = 0;
668: lb->r = r;
669: lb->hist = hist;
670: lb->len = 1;
671:
672: /* g : table of normal form as linear form */
673: MKNODE(g,lb,0);
674:
675: len = 1;
676: h = UPALLOC(dy);
677: for ( d = 1; ; d++ ) {
678: if ( d > c ){
679: return;
680: }
681: nf[d] = (int **)CALLOC(d+1,sizeof(int *));
682: len0 = len;
683: len += d+1;
684:
685: for ( i = d; i >= 0; i-- ) {
686: /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
687: /* nf(x^d) = nf(nf(x^(d-1))*x) */
688: r = (int *)CALLOC(dxdy,sizeof(int));
689: if ( i == 0 ) {
690: prev = nf[d-1][0];
691: bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
692:
693: /* create the head coeff */
694: for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
1.14 noro 695: iftogfs(prev[k],&s);
1.12 noro 696: COEF(h)[l] = (Num)s;
697: }
698: for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
699: DEG(h) = l;
700:
701: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
702: tmulup(rx[k],h,dy,&ut);
703: if ( ut )
704: for ( l = 0; l < dy; l++ ) {
705: s = (GFS)COEF(ut)[l];
706: if ( s ) {
707: u = CONT(s);
708: r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));
709: }
710: }
711: }
712: } else {
713: prev = nf[d-1][i-1];
714: for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
715: for ( l = 1; l < dy; l++ )
716: r[dyk+l] = prev[dyk+l-1];
717: }
718: }
719: nf[d][i] = r;
720: hist = (int *)CALLOC(len,sizeof(int));
721: hist[len0+i] = _onesf();
722: r1 = (int *)CALLOC(dxdy,sizeof(int));
723: bcopy(r,r1,dxdy*sizeof(int));
724: lnf(r1,hist,dxdy,len,g);
725: for ( k = 0; k < dxdy && !r1[k]; k++ );
726: if ( k == dxdy ) {
727: f = 0;
728: for ( k = j = 0; k <= d; k++ )
729: for ( i = 0; i <= k; i++, j++ )
730: if ( hist[j] ) {
1.14 noro 731: iftogfs(hist[j],&s);
1.12 noro 732: /* mono = s*x^(k-i)*y^i */
733: create_bmono((P)s,x,k-i,y,i,&mono);
734: addp(vl,f,mono,&f1);
735: f = f1;
736: }
737: *fr = f;
738: return;
739: } else {
740: lb = (struct lb *)CALLOC(1,sizeof(struct lb));
741: lb->pos = k;
742: lb->r = r1;
743: lb->hist = hist;
744: lb->len = len;
745: g = insert_lb(g,lb);
746: }
747: }
748: }
749: }
750:
751: static void create_bmono(P c,V x,int i,V y,int j,P *mono)
752: {
753: P t,s;
754:
755: if ( !i )
756: if ( !j )
757: t = c;
758: else {
759: /* c*y^j */
760: MKV(y,t);
761: COEF(DC(t)) = c;
762: STOQ(j,DEG(DC(t)));
763: }
764: else if ( !j ) {
765: /* c*x^i */
766: MKV(x,t);
767: COEF(DC(t)) = c;
768: STOQ(i,DEG(DC(t)));
769: } else {
770: MKV(y,s);
771: COEF(DC(s)) = c;
772: STOQ(j,DEG(DC(s)));
773: MKV(x,t);
774: COEF(DC(t)) = s;
775: STOQ(i,DEG(DC(t)));
776: }
777: *mono = t;
1.1 noro 778: }