[BACK]Return to Mgfs.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2018 / engine

Annotation of OpenXM_contrib2/asir2018/engine/Mgfs.c, Revision 1.1

1.1     ! noro        1: /* $OpenXM$ */
        !             2:
        !             3: #include "ca.h"
        !             4: #include "inline.h"
        !             5:
        !             6: extern int current_gfs_p;
        !             7: extern int up_kara_mag, current_gfs_q1;
        !             8: extern int *current_gfs_plus1;
        !             9:
        !            10: #if defined(__GNUC__)
        !            11: #define INLINE static inline
        !            12: #elif defined(VISUAL) || defined(__MINGW32__)
        !            13: #define INLINE __inline
        !            14: #else
        !            15: #define INLINE
        !            16: #endif
        !            17:
        !            18: INLINE int _ADDSF(int a,int b)
        !            19: {
        !            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:   }
        !            61: }
        !            62:
        !            63: INLINE int _MULSF(int a,int b)
        !            64: {
        !            65:   int c;
        !            66:
        !            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:   }
        !            79: }
        !            80:
        !            81: void addsfum(UM p1,UM p2,UM pr)
        !            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++ )
        !            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;
        !           106: }
        !           107:
        !           108: void subsfum(UM p1,UM p2,UM pr)
        !           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:
        !           143: void gcdsfum(UM p1,UM p2,UM pr)
        !           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: }
        !           167:
        !           168: void mulsfum(UM p1,UM p2,UM pr)
        !           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++ ) != 0 )
        !           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;
        !           187: }
        !           188:
        !           189: void mulssfum(UM p,int n,UM pr)
        !           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-- )
        !           196:     *dp = _MULSF(*sp,n);
        !           197: }
        !           198:
        !           199: void kmulsfum(UM n1,UM n2,UM nr)
        !           200: {
        !           201:   UM n,t,s,m,carry;
        !           202:   int d,d1,d2,len,i,l;
        !           203:   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 = (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));
        !           244: }
        !           245:
        !           246: void kmulsfummain(UM n1,UM n2,UM nr)
        !           247: {
        !           248:   int d1,d2,h,len,d;
        !           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:   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);
        !           279: }
        !           280:
        !           281: int divsfum(UM p1,UM p2,UM pq)
        !           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-- )
        !           296:       *pc1 = _MULSF(*pc1,inv);
        !           297:   } else
        !           298:     inv = _onesf();
        !           299:   for ( i = dd, ct = c1+d1; i >= 0; i-- )
        !           300:     if ( ( tmp = *ct-- ) != 0 ) {
        !           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;
        !           315: }
        !           316:
        !           317: void diffsfum(UM f,UM fd)
        !           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-- ) {
        !           324:     *dp = _MULSF(*sp,_itosf(i%current_gfs_p));
        !           325:   }
        !           326:   degum(fd,DEG(f) - 1);
        !           327: }
        !           328:
        !           329: void monicsfum(UM f)
        !           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-- )
        !           337:     *sp = _MULSF(*sp,inv);
        !           338: }
        !           339:
        !           340: int compsfum(UM a,UM b)
        !           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:   }
        !           361: }
        !           362:
        !           363: void addsfarray(int,int *,int *);
        !           364: void mulsfarray_trunc(int,int *,int *,int *);
        !           365:
        !           366: /* f1 = f1->c[0]+f1->c[1]*y+...,  f2 = f2->c[0]+f2->c[1]*y+... mod y^n */
        !           367:
        !           368: void mulsfbm(BM f1,BM f2,BM fr)
        !           369: {
        !           370:   UM mul,t,s;
        !           371:   int i,j,d1,d2,dy;
        !           372:
        !           373:   dy = MIN(DEG(f1),DEG(f2));
        !           374:
        !           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;
        !           392: }
        !           393:
        !           394: int degbm(BM f)
        !           395: {
        !           396:   int d,i,dy;
        !           397:
        !           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;
        !           403: }
        !           404:
        !           405: /* g += f */
        !           406:
        !           407: void addtosfbm(BM f,BM g)
        !           408: {
        !           409:   int i,d1,d2,dy;
        !           410:   UM t;
        !           411:
        !           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:   }
        !           423: }
        !           424:
        !           425: void eucsfum(UM f1,UM f2,UM a,UM b)
        !           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);
        !           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);
        !           455: }
        !           456:
        !           457: void shiftsfum(UM,int,UM);
        !           458:
        !           459: void shiftsflum(int n,LUM f,int ev)
        !           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:
        !           483: void shiftsfum(UM f,int a,UM g)
        !           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:
        !           518: /* f(y) -> f(y+a) */
        !           519:
        !           520: void shiftsfbm(BM f,int a)
        !           521: {
        !           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:   }
        !           551: }
        !           552:
        !           553: void clearbm(int n,BM f)
        !           554: {
        !           555:   int i,dy;
        !           556:   UM *c;
        !           557:
        !           558:   dy = DEG(f);
        !           559:   for ( c = COEF(f), i = 0; i <= dy; i++ )
        !           560:     clearum(c[i],n);
        !           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] ) != 0 ) {
        !           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++ ) {
        !           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:   }
        !           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;
        !           778: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>