Annotation of OpenXM_contrib2/asir2000/engine/Mgfs.c, Revision 1.20
1.20 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.19 2015/08/14 13:51:54 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: {
1.20 ! noro 20: if ( !a )
! 21: return b;
! 22: else if ( !b )
! 23: return a;
! 24:
! 25: a = IFTOF(a); b = IFTOF(b);
! 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:
! 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: }
1.5 noro 61: }
62:
1.11 noro 63: INLINE int _MULSF(int a,int b)
1.5 noro 64: {
1.20 ! noro 65: int c;
1.13 noro 66:
1.20 ! noro 67: if ( !a || !b )
! 68: return 0;
! 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 {
! 74: a = IFTOF(a) + IFTOF(b);
! 75: if ( a >= current_gfs_q1 )
! 76: a -= current_gfs_q1;
! 77: return FTOIF(a);
! 78: }
1.5 noro 79: }
1.1 noro 80:
1.11 noro 81: void addsfum(UM p1,UM p2,UM pr)
1.1 noro 82: {
1.20 ! noro 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++ )
! 99: cr[i] = _ADDSF(c1[i],c2[i]);
! 100: for ( ; i <= dmax; i++ )
! 101: cr[i] = c1[i];
! 102: if ( dmax == dmin )
! 103: degum(pr,dmax);
! 104: else
! 105: DEG(pr) = dmax;
1.1 noro 106: }
107:
1.11 noro 108: void subsfum(UM p1,UM p2,UM pr)
1.1 noro 109: {
1.20 ! noro 110: int *c1,*c2,*cr,i;
! 111: int dmax,dmin;
1.1 noro 112:
1.20 ! noro 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;
1.1 noro 141: }
1.20 ! noro 142:
1.11 noro 143: void gcdsfum(UM p1,UM p2,UM pr)
1.1 noro 144: {
1.20 ! noro 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: }
1.1 noro 166: }
1.11 noro 167:
168: void mulsfum(UM p1,UM p2,UM pr)
1.1 noro 169: {
1.20 ! noro 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++ )
! 184: if ( *pc1 )
! 185: *pcr = _ADDSF(_MULSF(*pc1,mul),*pcr);
! 186: DEG(pr) = d1 + d2;
1.1 noro 187: }
188:
1.11 noro 189: void mulssfum(UM p,int n,UM pr)
1.1 noro 190: {
1.20 ! noro 191: int *sp,*dp;
! 192: int i;
1.1 noro 193:
1.20 ! noro 194: for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
! 195: i >= 0; i--, dp--, sp-- )
! 196: *dp = _MULSF(*sp,n);
1.1 noro 197: }
1.20 ! noro 198:
1.11 noro 199: void kmulsfum(UM n1,UM n2,UM nr)
1.5 noro 200: {
1.20 ! noro 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(d1+d2+1);
! 224: carry = W_UMALLOC(d1+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));
1.5 noro 244: }
245:
1.11 noro 246: void kmulsfummain(UM n1,UM n2,UM nr)
1.5 noro 247: {
1.20 ! noro 248: int d1,d2,h,len,d;
! 249: UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
1.5 noro 250:
1.20 ! noro 251: d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
! 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);
! 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);
1.5 noro 279: }
280:
1.11 noro 281: int divsfum(UM p1,UM p2,UM pq)
1.1 noro 282: {
1.20 ! noro 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-- )
! 296: *pc1 = _MULSF(*pc1,inv);
! 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-- )
! 303: *pct = _ADDSF(_MULSF(*pc1,tmp),*pct);
! 304: }
! 305: if ( inv != _onesf() ) {
! 306: for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
! 307: *pc1 = _MULSF(*pc1,inv);
! 308: for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )
! 309: *pc1 = _MULSF(*pc1,hd);
! 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;
1.1 noro 315: }
316:
1.11 noro 317: void diffsfum(UM f,UM fd)
1.1 noro 318: {
1.20 ! noro 319: int *dp,*sp;
! 320: int i;
1.1 noro 321:
1.20 ! noro 322: for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
! 323: i >= 1; i--, dp--, sp-- ) {
! 324: *dp = _MULSF(*sp,_itosf(i%current_gfs_p));
! 325: }
! 326: degum(fd,DEG(f) - 1);
1.1 noro 327: }
328:
1.11 noro 329: void monicsfum(UM f)
1.1 noro 330: {
1.20 ! noro 331: int *sp;
! 332: int i,inv;
1.1 noro 333:
1.20 ! noro 334: i = DEG(f); sp = COEF(f)+i;
! 335: inv = _invsf(*sp);
! 336: for ( ; i >= 0; i--, sp-- )
! 337: *sp = _MULSF(*sp,inv);
1.10 noro 338: }
339:
1.11 noro 340: int compsfum(UM a,UM b)
1.10 noro 341: {
1.20 ! noro 342: int i,da,db;
1.10 noro 343:
1.20 ! noro 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: {
1.20 ! noro 370: UM mul,t,s;
! 371: int i,j,d1,d2,dy;
1.4 noro 372:
1.20 ! noro 373: dy = MIN(DEG(f1),DEG(f2));
1.4 noro 374:
1.20 ! noro 375: d1 = degbm(f1);
! 376: d2 = degbm(f2);
! 377: t = W_UMALLOC(d1+d2);
! 378: s = W_UMALLOC(d1+d2);
! 379:
! 380: for ( i = 0; i <= dy; i++ ) {
! 381: mul = COEF(f2)[i];
! 382: if ( DEG(mul) >= 0 )
! 383: for ( j = 0; i+j <= dy; j++ ) {
! 384: if ( COEF(f1)[j] ) {
! 385: kmulsfum(COEF(f1)[j],mul,t);
! 386: addsfum(t,COEF(fr)[i+j],s);
! 387: cpyum(s,COEF(fr)[i+j]);
! 388: }
! 389: }
! 390: }
! 391: DEG(fr) = dy;
1.4 noro 392: }
393:
1.11 noro 394: int degbm(BM f)
1.6 noro 395: {
1.20 ! noro 396: int d,i,dy;
1.6 noro 397:
1.20 ! noro 398: dy = DEG(f);
! 399: d = DEG(COEF(f)[0]);
! 400: for ( i = 1; i <= dy; i++ )
! 401: d = MAX(DEG(COEF(f)[i]),d);
! 402: return d;
1.6 noro 403: }
404:
1.4 noro 405: /* g += f */
406:
1.11 noro 407: void addtosfbm(BM f,BM g)
1.4 noro 408: {
1.20 ! noro 409: int i,d1,d2,dy;
! 410: UM t;
1.4 noro 411:
1.20 ! 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);
! 418: t = W_UMALLOC(MAX(d1,d2));
! 419: for ( i = 0; i <= dy; i++ ) {
! 420: addsfum(COEF(f)[i],COEF(g)[i],t);
! 421: cpyum(t,COEF(g)[i]);
! 422: }
1.4 noro 423: }
424:
1.11 noro 425: void eucsfum(UM f1,UM f2,UM a,UM b)
1.2 noro 426: {
1.20 ! noro 427: UM g1,g2,a1,a2,a3,wm,q,tum;
! 428: int d,dr;
1.2 noro 429:
1.20 ! noro 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);
! 435: DEG(a1) = 0; COEF(a1)[0] = _onesf(); DEG(a2) = -1;
! 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: }
! 444: if ( COEF(g1)[0] != _onesf() )
! 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: {
1.20 ! noro 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: }
1.3 noro 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: {
1.20 ! noro 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: }
1.3 noro 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.20 ! noro 522: int i,j,d,dy;
! 523: UM pwr,ya,w,t,s;
! 524: UM *c;
! 525:
! 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);
! 549: }
! 550: }
1.4 noro 551: }
552:
1.11 noro 553: void clearbm(int n,BM f)
1.4 noro 554: {
1.20 ! noro 555: int i,dy;
! 556: UM *c;
1.4 noro 557:
1.20 ! noro 558: dy = DEG(f);
! 559: for ( c = COEF(f), i = 0; i <= dy; i++ )
! 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 {
1.20 ! noro 566: int pos,len;
! 567: int *r;
! 568: int *hist;
1.12 noro 569: };
570:
571: static NODE insert_lb(NODE g,struct lb *a)
572: {
1.20 ! noro 573: NODE prev,cur,n;
1.12 noro 574:
1.20 ! noro 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;
1.12 noro 593: }
594:
595: static void lnf(int *r,int *h,int n,int len,NODE g)
596: {
1.20 ! noro 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: }
1.12 noro 627: }
628:
629: void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
630: {
1.20 ! noro 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++ ) {
! 695: iftogfs(prev[k],&s);
! 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] ) {
! 731: iftogfs(hist[j],&s);
! 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: }
1.12 noro 749: }
750:
751: static void create_bmono(P c,V x,int i,V y,int j,P *mono)
752: {
1.20 ! noro 753: P t,s;
1.12 noro 754:
1.20 ! noro 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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