[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.5

1.5     ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly3.c,v 1.4 2002/02/04 07:58:28 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:   */
                     34:   DList = (int *)malloc(sizeof(int)*Lsize*N0);
                     35:   MList = (int *)malloc(sizeof(int)*Lsize*N0);
                     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:
                    101:   int n,k,c,l,q,i,m;
                    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;
                    109:   for (k=Plist-1; k>=0; k--) {
                    110:     /* coeff */
                    111:     a = coeffCopy(CList[k]);
                    112:     Cmult(a,a,f->coeffp);
                    113:     if (isZero(a)) goto no;
                    114:     /* initialize tmp */
                    115:     for (i=0; i<n; i++) {
                    116:       tmp.e[i] = f->m->e[i];
                    117:     }
                    118:     if (Homogenize) {
                    119:       tmp.e[0].D += EList[k]; /* homogenization.
1.3       takayama  120:                                  e[0].D will be added later. */
1.1       maekawa   121:     }
                    122:
                    123:     /* from m to n:  Differential variables. */
                    124:     for (i=0; i<Maxv; i++) {
                    125:       CiiPoch(Tc,tmp.e[V[i]].x,DList[I(k,i)]);
                    126:       /* Tc->val.bigp is read only */
                    127:       Cmult(a,a,Tc);
                    128:       /*printf("k=%d V[i]=%d a=%s\n",k,V[i],coeffToString(a));*/
                    129:       if (isZero(a)) goto no;
                    130:       tmp.e[V[i]].D += MList[I(k,i)]; /* monomial add */
                    131:       tmp.e[V[i]].x -= DList[I(k,i)]; /* differentiate */
                    132:     }
                    133:
                    134:
                    135:     /* difference variables are commutative */
                    136:
                    137:     /* q-variables. Compute q before updating tmp.e[i]. */
                    138:     if (l-c > 0) {
                    139:       q =0;
                    140:       for (i=c; i<l; i++) {
1.3       takayama  141:         q += (e[i].D)*(tmp.e[i].x);  /* Don't repeat these things. */
                    142:         tmp.e[i].D += e[i].D;
1.1       maekawa   143:       }
                    144:       /*printf("l=%d, q=%d\n",l,q);*/
                    145:       if (ringp->next == (struct ring *)NULL) {
1.3       takayama  146:         tmp.e[0].x += q;
1.1       maekawa   147:       }else{
1.3       takayama  148:         Cmult(a,a,polyToCoeff(cxx(1,0,q,ringp->next),ringp));
                    149:         /* x[0]^q */
1.1       maekawa   150:       }
                    151:     }
                    152:
                    153:     /* Update tmp.e[i].x */
                    154:     for (i=0; i<n; i++) {
                    155:       tmp.e[i].x += e[i].x;
                    156:     }
                    157:
                    158:
                    159:     /* commutative variables */
                    160:     for (i=0; i<c; i++) {
                    161:       tmp.e[i].D += e[i].D;
                    162:     }
                    163:     /* Difference variables */
                    164:     for (i=l; i<m; i++) {
                    165:       tmp.e[i].D += e[i].D;
                    166:     }
                    167:     /***** x_i ----> x_i + e[i].D h  Substitution and homogenization */
                    168:     /*** They will be done in mpMult_diff() */
                    169:
                    170:     mm = newCell(a,monomialCopy(&tmp));
                    171:     RList[k]->next = mm;
                    172:     RList[k] = RList[k]->next;
                    173:   no: ;
                    174:   }
                    175: }
                    176:
                    177: /* Note that you cannot call mpMult_diff recursively. */
                    178: /* Note also that mpMult_diff assumes coefficients and Dx commutes each other*/
1.5     ! takayama  179: POLY mpMult_diff(POLY f,POLY g)
1.1       maekawa   180: {
                    181:   int k;
                    182:   POLY r,temp;
                    183:
                    184:   /* printf("mpMult_diff(%s,%s)\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
                    185:
                    186:   if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
                    187:   checkRing(f,g);
                    188:   Tc = coeffCopy(f->coeffp);
                    189:
                    190:   if (isConstant(f)) return(cpMult(f->coeffp,g));
                    191:
                    192:   makeTable(f->coeffp,f->m->e,f->m->ringp);
                    193:   /* outputTable(); */
                    194:   for (k=0; k<Plist; k++) {
                    195:     RList[k] = RListRoot[k];
                    196:     RList[k]->next = POLYNULL;
                    197:   }
                    198:
                    199:   while (g != POLYNULL) {
                    200:     monomialMult_diff(f->m->e,g);
                    201:     g = g->next;
                    202:   }
                    203:   r = POLYNULL;
                    204:   for (k=0; k<Plist; k++) {
                    205:     temp = RListRoot[k]->next;
                    206:     r = ppAddv(r,temp);
                    207:   }
                    208:
                    209:
                    210:   /***** x_i ----> x_i + e[i].D h  Substitution and homogenization
                    211:          for difference variables */
                    212:   /*** They are implemented in _difference, but slow.*/
                    213:
                    214:   return(r);
                    215: }
                    216:
1.5     ! takayama  217: POLY mpMult_difference_org(POLY f,POLY g)
1.1       maekawa   218: {
                    219:   POLY r;
                    220:   int m,l;
                    221:   POLY lRule[N0];
                    222:   POLY rRule[N0];
                    223:   int size;
                    224:   int i;
                    225:
                    226:   if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
                    227:   checkRing(f,g);
                    228:   m = f->m->ringp->m;
                    229:   l = f->m->ringp->l;
                    230:
                    231:   r = mpMult_diff(f,g);
                    232:
                    233:   /***** x_i ----> x_i + e[i].D h  Substitution and homogenization
                    234:          for difference variables */
                    235:   size = 0;
                    236:   for (i=l; i<m; i++) {
                    237:     if (f->m->e[i].D) {
                    238:       lRule[size] = cxx(1,i,1,f->m->ringp);
                    239:       if (Homogenize) {
1.3       takayama  240:         rRule[size] = ppAdd(cxx(1,i,1,f->m->ringp),cdd(f->m->e[i].D,0,1,f->m->ringp));
                    241:         /* x_i               + e[i].D  h */
1.1       maekawa   242:       }else{
1.3       takayama  243:         rRule[size] = ppAdd(cxx(1,i,1,f->m->ringp),cdd(f->m->e[i].D,0,0,f->m->ringp));
                    244:         /* x_i               + e[i].D   */
                    245:       }
1.1       maekawa   246:       size++;
                    247:     }
                    248:   }
                    249:
                    250:   /* It's a dirty trick. */
                    251:   r = replace_poly(r,lRule,rRule,size);
                    252:
                    253:   return(r);
                    254: }
                    255:
                    256: outputTable() {
                    257:   int i,j;
                    258:   printf("Maxv = %d Plist=%d\n",Maxv,Plist);
                    259:   for (i=0; i<Maxv; i++) printf("%5d",V[i]);
                    260:   printf("\n------ DList --------------\n");
                    261:   for (i=0; i<Plist; i++) {
                    262:     for (j=0; j<Maxv; j++) {
                    263:       printf("%5d",DList[I(i,j)]);
                    264:     }
                    265:     putchar('\n');
                    266:   }
                    267:   printf("\n--------- MList ------------\n");
                    268:   for (i=0; i<Plist; i++) {
                    269:     for (j=0; j<Maxv; j++) {
                    270:       printf("%5d",MList[I(i,j)]);
                    271:     }
                    272:     printf(" |  e=%5d M=%1d c=%s\n",EList[i],Mark[i],coeffToString(CList[i]));
                    273:   }
                    274: }

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