[BACK]Return to poly3.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Kan

Annotation of OpenXM/src/kan96xx/Kan/poly3.c, Revision 1.7

1.7     ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly3.c,v 1.6 2002/09/08 10:49:50 takayama Exp $ */
1.1       maekawa     2: #include <stdio.h>
                      3: #include "datatype.h"
                      4: #include "extern2.h"
                      5:
                      6: int Homogenize = 1;
                      7:
                      8: #define I(i,j) (i*N0+j)
                      9: #define LSIZE 5000
                     10: static int V[N0];  /* index of variables */
                     11: static struct coeff ** CList;
                     12: static int *EList;
                     13: static int *DList;
                     14: static int *MList;
                     15: static int *Mark;
                     16: static int Lsize;
                     17: static int Plist;
                     18: static int Maxv;
                     19: static POLY *RList;
                     20: static POLY *RListRoot;
                     21: static struct coeff *Tc = (struct coeff *)NULL;
1.3       takayama   22: /* It is initialized in mpMult_diff() */
1.1       maekawa    23:
                     24: void initT(void) {
                     25:   int i;
                     26:   Lsize = LSIZE;
                     27:   CList = (struct coeff **)sGC_malloc(sizeof(struct coeff *)*Lsize);
                     28:   EList = (int *)sGC_malloc(sizeof(int)*Lsize);
                     29:   Mark = (int *)sGC_malloc(sizeof(int)*Lsize);
                     30:   /* The following line causes the warning 'needed to allocate blacklisted..'
1.3       takayama   31:      DList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
                     32:      MList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
1.1       maekawa    33:   */
1.7     ! takayama   34:   DList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
        !            35:   MList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
1.1       maekawa    36:
                     37:   RList = (POLY *)sGC_malloc(sizeof(POLY)*Lsize);
                     38:   RListRoot = (POLY *)sGC_malloc(sizeof(POLY)*Lsize);
                     39:   for (i=0; i<Lsize; i++)
                     40:     RListRoot[i] = newCell((struct coeff *)NULL,(MONOMIAL)NULL);
                     41: }
                     42:
                     43: void makeTable(c,e,ringp)
1.3       takayama   44:      struct coeff *c; /* read only */
                     45:      struct exps e[];
                     46:      struct ring *ringp;
1.1       maekawa    47: {
                     48:   int i,j,k,p,q,deg,m,n;
                     49:   m = ringp->m; n = ringp->n;
                     50:   /* initialize */
                     51:   Maxv = 0; Plist = 1; deg = 0;
                     52:   for (i=m; i<n; i++) {
                     53:     if (e[i].D != 0) {
                     54:       V[Maxv] = i;
                     55:       DList[I(0,Maxv)] = e[i].D;
                     56:       MList[I(0,Maxv)] = 0;
                     57:       deg += e[i].D;
                     58:       Maxv++;
                     59:     }
                     60:   }
                     61:   CList[0] = coeffCopy(c);
                     62:   EList[0] = deg*2; Mark[0] = 0;
                     63:
                     64:   for (i=0; i<Maxv; i++) {
                     65:     k = Plist;
                     66:     /* Copy j-th row to k-th row and modify it. */
                     67:     for (j=0; j<Plist; j++) {
                     68:       for (q=1; q<=DList[I(j,i)]; q++) {
1.3       takayama   69:         for (p=0; p<Maxv; p++) { /* copy */
                     70:           DList[I(k,p)] = DList[I(j,p)];
                     71:           MList[I(k,p)] = MList[I(j,p)];
                     72:         }
                     73:         /* modify */
                     74:         DList[I(k,i)] -= q;
                     75:         MList[I(k,i)] += q;
                     76:
                     77:         CiiComb(Tc,DList[I(j,i)],q);
                     78:         /* Tc->val.bigp is read only. */
                     79:         CList[k] = coeffCopy(Tc);
                     80:         Cmult(CList[k],CList[k],CList[j]);
                     81:         /*CList[k] = normalize(CList[j]*BiiComb(DList[I(j,i)],q));*/
                     82:
                     83:         EList[k] = EList[j]-2*q;
                     84:         Mark[k] = 0;
                     85:         k++;
                     86:         if (k>= Lsize) {
                     87:           errorPoly("makeTable(): Lsize is not large enough.\n");
                     88:         }
1.1       maekawa    89:       }
                     90:     }
                     91:     Plist = k;
                     92:   }
                     93: }
                     94:
                     95: void monomialMult_diff(e,f)
1.3       takayama   96:      struct exps e[];
                     97:      POLY f;
                     98:      /* (e) * f = [Plist] monomials  */
1.1       maekawa    99: {
                    100:
1.6       takayama  101:   int n,k,c,l,q,i,m, weightedHomogenization;
1.1       maekawa   102:   struct coeff *a;
                    103:   struct monomial tmp;
                    104:   struct ring *ringp;
                    105:   POLY mm;
                    106:
                    107:   tmp.ringp = ringp = f->m->ringp;
                    108:   n = ringp->n; c = ringp->c; l = ringp->l; m = ringp->m;
1.6       takayama  109:   weightedHomogenization = ringp->weightedHomogenization;
1.1       maekawa   110:   for (k=Plist-1; k>=0; k--) {
                    111:     /* coeff */
                    112:     a = coeffCopy(CList[k]);
                    113:     Cmult(a,a,f->coeffp);
                    114:     if (isZero(a)) goto no;
                    115:     /* initialize tmp */
                    116:     for (i=0; i<n; i++) {
                    117:       tmp.e[i] = f->m->e[i];
                    118:     }
1.6       takayama  119:     if ((!weightedHomogenization) && Homogenize) {
1.1       maekawa   120:       tmp.e[0].D += EList[k]; /* homogenization.
1.3       takayama  121:                                  e[0].D will be added later. */
1.6       takayama  122:     }else if (weightedHomogenization && Homogenize) {
                    123:       tmp.e[0].D += EList[k]/2 ; /* homogenization.  Weight is (1,0) (special).
                    124:                                                                  */
                    125:        }
1.1       maekawa   126:
                    127:     /* from m to n:  Differential variables. */
                    128:     for (i=0; i<Maxv; i++) {
                    129:       CiiPoch(Tc,tmp.e[V[i]].x,DList[I(k,i)]);
                    130:       /* Tc->val.bigp is read only */
                    131:       Cmult(a,a,Tc);
                    132:       /*printf("k=%d V[i]=%d a=%s\n",k,V[i],coeffToString(a));*/
                    133:       if (isZero(a)) goto no;
                    134:       tmp.e[V[i]].D += MList[I(k,i)]; /* monomial add */
                    135:       tmp.e[V[i]].x -= DList[I(k,i)]; /* differentiate */
                    136:     }
                    137:
                    138:
                    139:     /* difference variables are commutative */
                    140:
                    141:     /* q-variables. Compute q before updating tmp.e[i]. */
                    142:     if (l-c > 0) {
                    143:       q =0;
                    144:       for (i=c; i<l; i++) {
1.3       takayama  145:         q += (e[i].D)*(tmp.e[i].x);  /* Don't repeat these things. */
                    146:         tmp.e[i].D += e[i].D;
1.1       maekawa   147:       }
                    148:       /*printf("l=%d, q=%d\n",l,q);*/
                    149:       if (ringp->next == (struct ring *)NULL) {
1.3       takayama  150:         tmp.e[0].x += q;
1.1       maekawa   151:       }else{
1.3       takayama  152:         Cmult(a,a,polyToCoeff(cxx(1,0,q,ringp->next),ringp));
                    153:         /* x[0]^q */
1.1       maekawa   154:       }
                    155:     }
                    156:
                    157:     /* Update tmp.e[i].x */
                    158:     for (i=0; i<n; i++) {
                    159:       tmp.e[i].x += e[i].x;
                    160:     }
                    161:
                    162:
                    163:     /* commutative variables */
                    164:     for (i=0; i<c; i++) {
                    165:       tmp.e[i].D += e[i].D;
                    166:     }
                    167:     /* Difference variables */
                    168:     for (i=l; i<m; i++) {
                    169:       tmp.e[i].D += e[i].D;
                    170:     }
                    171:     /***** x_i ----> x_i + e[i].D h  Substitution and homogenization */
                    172:     /*** They will be done in mpMult_diff() */
                    173:
                    174:     mm = newCell(a,monomialCopy(&tmp));
                    175:     RList[k]->next = mm;
                    176:     RList[k] = RList[k]->next;
                    177:   no: ;
                    178:   }
                    179: }
                    180:
                    181: /* Note that you cannot call mpMult_diff recursively. */
                    182: /* Note also that mpMult_diff assumes coefficients and Dx commutes each other*/
1.5       takayama  183: POLY mpMult_diff(POLY f,POLY g)
1.1       maekawa   184: {
                    185:   int k;
                    186:   POLY r,temp;
                    187:
                    188:   /* printf("mpMult_diff(%s,%s)\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
                    189:
                    190:   if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
                    191:   checkRing(f,g);
                    192:   Tc = coeffCopy(f->coeffp);
                    193:
                    194:   if (isConstant(f)) return(cpMult(f->coeffp,g));
                    195:
                    196:   makeTable(f->coeffp,f->m->e,f->m->ringp);
                    197:   /* outputTable(); */
                    198:   for (k=0; k<Plist; k++) {
                    199:     RList[k] = RListRoot[k];
                    200:     RList[k]->next = POLYNULL;
                    201:   }
                    202:
                    203:   while (g != POLYNULL) {
                    204:     monomialMult_diff(f->m->e,g);
                    205:     g = g->next;
                    206:   }
                    207:   r = POLYNULL;
                    208:   for (k=0; k<Plist; k++) {
                    209:     temp = RListRoot[k]->next;
                    210:     r = ppAddv(r,temp);
                    211:   }
                    212:
                    213:
                    214:   /***** x_i ----> x_i + e[i].D h  Substitution and homogenization
                    215:          for difference variables */
                    216:   /*** They are implemented in _difference, but slow.*/
                    217:
                    218:   return(r);
                    219: }
                    220:
1.5       takayama  221: POLY mpMult_difference_org(POLY f,POLY g)
1.1       maekawa   222: {
                    223:   POLY r;
                    224:   int m,l;
                    225:   POLY lRule[N0];
                    226:   POLY rRule[N0];
                    227:   int size;
                    228:   int i;
                    229:
                    230:   if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
                    231:   checkRing(f,g);
                    232:   m = f->m->ringp->m;
                    233:   l = f->m->ringp->l;
                    234:
                    235:   r = mpMult_diff(f,g);
                    236:
                    237:   /***** x_i ----> x_i + e[i].D h  Substitution and homogenization
                    238:          for difference variables */
                    239:   size = 0;
                    240:   for (i=l; i<m; i++) {
                    241:     if (f->m->e[i].D) {
                    242:       lRule[size] = cxx(1,i,1,f->m->ringp);
                    243:       if (Homogenize) {
1.3       takayama  244:         rRule[size] = ppAdd(cxx(1,i,1,f->m->ringp),cdd(f->m->e[i].D,0,1,f->m->ringp));
                    245:         /* x_i               + e[i].D  h */
1.1       maekawa   246:       }else{
1.3       takayama  247:         rRule[size] = ppAdd(cxx(1,i,1,f->m->ringp),cdd(f->m->e[i].D,0,0,f->m->ringp));
                    248:         /* x_i               + e[i].D   */
                    249:       }
1.1       maekawa   250:       size++;
                    251:     }
                    252:   }
                    253:
                    254:   /* It's a dirty trick. */
                    255:   r = replace_poly(r,lRule,rRule,size);
                    256:
                    257:   return(r);
                    258: }
                    259:
                    260: outputTable() {
                    261:   int i,j;
                    262:   printf("Maxv = %d Plist=%d\n",Maxv,Plist);
                    263:   for (i=0; i<Maxv; i++) printf("%5d",V[i]);
                    264:   printf("\n------ DList --------------\n");
                    265:   for (i=0; i<Plist; i++) {
                    266:     for (j=0; j<Maxv; j++) {
                    267:       printf("%5d",DList[I(i,j)]);
                    268:     }
                    269:     putchar('\n');
                    270:   }
                    271:   printf("\n--------- MList ------------\n");
                    272:   for (i=0; i<Plist; i++) {
                    273:     for (j=0; j<Maxv; j++) {
                    274:       printf("%5d",MList[I(i,j)]);
                    275:     }
                    276:     printf(" |  e=%5d M=%1d c=%s\n",EList[i],Mark[i],coeffToString(CList[i]));
                    277:   }
                    278: }

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