Annotation of OpenXM/src/kan96xx/Kan/poly3a.c, Revision 1.1
1.1 ! takayama 1: /* $OpenXM$ */
! 2: #include <stdio.h>
! 3: #include "datatype.h"
! 4: #include "extern2.h"
! 5:
! 6: extern int Homogenize;
! 7:
! 8: #define DIFFERENTIAL 0
! 9: #define DIFFERENCE 1
! 10: #define I(i,j) (i*N0+j)
! 11: #define LSIZE 5000
! 12: #define MULTTABLE_MAX 10
! 13: struct multTable {
! 14: int V[N0]; /* index of variables */
! 15: int Vtype[N0];
! 16: struct coeff ** CList; /* binomial coeff */
! 17: int *EList; /* for homogenization */
! 18: int *DList;
! 19: int *MList;
! 20: int *Mark;
! 21: int Lsize;
! 22: int Plist;
! 23: int Maxv;
! 24: POLY *RList;
! 25: POLY *RListRoot;
! 26: struct coeff *Tc;
! 27: int myplace;
! 28: };
! 29: static struct multTable * MultTable[MULTTABLE_MAX];
! 30: static int MultTableMark[MULTTABLE_MAX];
! 31: static struct multTable *newMultTable(struct coeff *cf);
! 32: static void freeMultTable(struct multTable *mp);
! 33: static struct multTable *initMT(int place);
! 34: static void makeTable_differential_difference(struct coeff *c, struct exps e[],
! 35: struct ring *ringp, struct multTable *mt);
! 36: static void monomialMult_differential_difference(struct exps e[],
! 37: POLY f, struct multTable *mt);
! 38: static void outputMultTable(struct multTable *mt);
! 39:
! 40:
! 41: POLY muMult_differential_difference(POLY f,POLY g);
! 42:
! 43: /*
! 44: These local variables should be set and restored from and to multTable.
! 45: */
! 46: #define mt_declare_locals() \
! 47: int *V; int *Vtype; int Plist; int Maxv; \
! 48: struct coeff ** CList; \
! 49: int *EList; int *DList; int *MList; int *Mark; int Lsize; \
! 50: POLY *RList; POLY *RListRoot; struct coeff *Tc;
! 51:
! 52: #define mt_set_mt(mt) \
! 53: {mt->EList = EList; \
! 54: mt->DList = DList; \
! 55: mt->MList = MList; \
! 56: mt->Mark = Mark; \
! 57: mt->Lsize = Lsize; \
! 58: mt->RList = RList; \
! 59: mt->RListRoot = RListRoot; \
! 60: mt->Tc = Tc; \
! 61: mt->Plist = Plist; \
! 62: mt->Maxv = Maxv; \
! 63: mt->CList = CList; \
! 64: }
! 65:
! 66: #define mt_get_mt(mt) \
! 67: {EList = mt->EList; \
! 68: DList = mt->DList; \
! 69: MList = mt->MList; \
! 70: Mark = mt->Mark; \
! 71: Lsize = mt->Lsize; \
! 72: RList = mt->RList; \
! 73: RListRoot = mt->RListRoot; \
! 74: Tc = mt->Tc; \
! 75: Plist = mt->Plist; \
! 76: Maxv = mt->Maxv; \
! 77: V = mt->V; \
! 78: Vtype = mt->Vtype; \
! 79: CList = mt->CList; \
! 80: }
! 81:
! 82:
! 83:
! 84: static struct multTable *initMT(int place) {
! 85: int i;
! 86: mt_declare_locals()
! 87: struct multTable *mt;
! 88:
! 89: Lsize = LSIZE;
! 90: CList = (struct coeff **)sGC_malloc(sizeof(struct coeff *)*Lsize);
! 91: EList = (int *)sGC_malloc(sizeof(int)*Lsize);
! 92: Mark = (int *)sGC_malloc(sizeof(int)*Lsize);
! 93: DList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
! 94: MList = (int *)sGC_malloc(sizeof(int)*Lsize*N0);
! 95: RList = (POLY *)sGC_malloc(sizeof(POLY)*Lsize);
! 96: RListRoot = (POLY *)sGC_malloc(sizeof(POLY)*Lsize);
! 97: for (i=0; i<Lsize; i++)
! 98: RListRoot[i] = newCell((struct coeff *)NULL,(MONOMIAL)NULL);
! 99: Tc = NULL;
! 100:
! 101: mt = (struct multTable *)sGC_malloc(sizeof(struct multTable));
! 102: mt_set_mt(mt);
! 103: mt->myplace = place;
! 104: return mt;
! 105: }
! 106:
! 107: static struct multTable *newMultTable(struct coeff *cf) {
! 108: static int p = -1;
! 109: int i;
! 110: if (p+1 >= MULTTABLE_MAX-1) {
! 111: errorPoly("newMultTable(): MULTTABLE_MAX is not large enough.\n");
! 112: }
! 113: if (p<0) {
! 114: p = 0;
! 115: MultTable[p] = initMT(p); MultTableMark[p] = 1;
! 116: MultTable[p]->Tc = coeffCopy(cf);
! 117: p++; return( MultTable[p-1] );
! 118: }
! 119: for (i=0; i<p; i++) {
! 120: if (! MultTableMark[i]) {
! 121: MultTable[i]->Tc = coeffCopy(cf);
! 122: return( MultTable[i] );
! 123: }
! 124: }
! 125: MultTable[p] = initMT(p); MultTableMark[p] = 1;
! 126: MultTable[p]->Tc = coeffCopy(cf);
! 127: p++; return( MultTable[p-1] );
! 128: }
! 129: static void freeMultTable(struct multTable *mp) {
! 130: MultTableMark[mp->myplace] = 0;
! 131: }
! 132:
! 133:
! 134: void makeTable_differential_difference(c,e,ringp,mt)
! 135: struct coeff *c; /* read only */
! 136: struct exps e[];
! 137: struct ring *ringp;
! 138: struct multTable *mt;
! 139: {
! 140: int i,j,k,p,q,deg,m,n,l,t;
! 141: mt_declare_locals()
! 142:
! 143: mt_get_mt(mt);
! 144: m = ringp->m; n = ringp->n; l = ringp->l;
! 145: /* initialize */
! 146: Maxv = 0; Plist = 1; deg = 0;
! 147:
! 148: for (i=l; i<m; i++) {
! 149: if (e[i].D != 0) {
! 150: V[Maxv] = i;
! 151: Vtype[Maxv] = DIFFERENCE;
! 152: DList[I(0,Maxv)] = e[i].D;
! 153: MList[I(0,Maxv)] = 0;
! 154: deg += e[i].D;
! 155: Maxv++;
! 156: }
! 157: }
! 158: for (i=m; i<n; i++) {
! 159: if (e[i].D != 0) {
! 160: V[Maxv] = i;
! 161: Vtype[Maxv] = DIFFERENTIAL;
! 162: DList[I(0,Maxv)] = e[i].D;
! 163: MList[I(0,Maxv)] = 0;
! 164: deg += e[i].D;
! 165: Maxv++;
! 166: }
! 167: }
! 168: CList[0] = coeffCopy(c);
! 169: Mark[0] = 0;
! 170:
! 171: for (i=0; i<Maxv; i++) {
! 172: k = Plist;
! 173: /* Copy j-th row to k-th row and modify it. */
! 174: for (j=0; j<Plist; j++) {
! 175: for (q=1; q<=DList[I(j,i)]; q++) {
! 176: for (p=0; p<Maxv; p++) { /* copy */
! 177: DList[I(k,p)] = DList[I(j,p)];
! 178: MList[I(k,p)] = MList[I(j,p)];
! 179: }
! 180: /* modify */
! 181: DList[I(k,i)] -= q;
! 182: MList[I(k,i)] += q;
! 183:
! 184: CiiComb(Tc,DList[I(j,i)],q);
! 185: /* Tc->val.bigp is read only. */
! 186: CList[k] = coeffCopy(Tc);
! 187: Cmult(CList[k],CList[k],CList[j]);
! 188: /*CList[k] = normalize(CList[j]*BiiComb(DList[I(j,i)],q));*/
! 189:
! 190: Mark[k] = 0;
! 191: k++;
! 192: if (k>= Lsize) {
! 193: errorPoly("makeTable_differential_difference(): Lsize is not large enough.\n");
! 194: }
! 195: }
! 196: }
! 197: Plist = k;
! 198: }
! 199:
! 200: if (Homogenize == 1) {
! 201: for (i=0; i<Plist; i++) {
! 202: t = deg;
! 203: for (j=0; j<Maxv; j++) {
! 204: t -= MList[I(i,j)];
! 205: }
! 206: for (j=0; j<Maxv; j++) {
! 207: if (Vtype[j] == DIFFERENTIAL) {
! 208: t += DList[I(i,j)];
! 209: }
! 210: }
! 211: EList[i] = t;
! 212: }
! 213: }else if (Homogenize == 2) { /* Dx x = x Dx + h */
! 214: for (i=0; i<Plist; i++) {
! 215: t = deg;
! 216: for (j=0; j<Maxv; j++) {
! 217: t -= MList[I(i,j)];
! 218: }
! 219: EList[i] = t;
! 220: }
! 221: }
! 222:
! 223:
! 224: mt_set_mt(mt);
! 225:
! 226: }
! 227:
! 228:
! 229: void monomialMult_differential_difference(e,f,mt)
! 230: struct exps e[];
! 231: POLY f;
! 232: /* (e) * f = [Plist] monomials */
! 233: struct multTable *mt;
! 234: {
! 235:
! 236: int n,k,c,l,q,i,m;
! 237: struct coeff *a;
! 238: struct monomial tmp;
! 239: struct ring *ringp;
! 240: POLY mm;
! 241: mt_declare_locals()
! 242:
! 243: mt_get_mt(mt);
! 244:
! 245: tmp.ringp = ringp = f->m->ringp;
! 246: n = ringp->n; c = ringp->c; l = ringp->l; m = ringp->m;
! 247: for (k=Plist-1; k>=0; k--) {
! 248: /* coeff */
! 249: a = coeffCopy(CList[k]);
! 250: Cmult(a,a,f->coeffp);
! 251: if (isZero(a)) goto no;
! 252: /* initialize tmp */
! 253: for (i=0; i<n; i++) {
! 254: tmp.e[i] = f->m->e[i];
! 255: }
! 256: if (Homogenize) {
! 257: tmp.e[0].D += EList[k]; /* homogenization.
! 258: e[0].D will be added later. */
! 259: }
! 260:
! 261: /* from m to n: Differential and difference variables. */
! 262: for (i=0; i<Maxv; i++) {
! 263: if (Vtype[i] == DIFFERENTIAL) {
! 264: CiiPoch(Tc,tmp.e[V[i]].x,DList[I(k,i)]);
! 265: /* Tc->val.bigp is read only */
! 266: Cmult(a,a,Tc);
! 267: /*printf("k=%d V[i]=%d a=%s\n",k,V[i],coeffToString(a));*/
! 268: if (isZero(a)) goto no;
! 269: tmp.e[V[i]].D += MList[I(k,i)]; /* monomial add for D. */
! 270: tmp.e[V[i]].x -= DList[I(k,i)]; /* differentiate */
! 271: }else{
! 272: CiiPower(Tc,tmp.e[V[i]].x,DList[I(k,i)]);
! 273: /* Tc->val.bigp is read only */
! 274: Cmult(a,a,Tc);
! 275: /*printf("k=%d V[i]=%d a=%s\n",k,V[i],coeffToString(a));*/
! 276: if (isZero(a)) goto no;
! 277: tmp.e[V[i]].D += MList[I(k,i)]; /* monomial add for D. */
! 278: }
! 279: }
! 280:
! 281:
! 282: /* q-variables. Compute q before updating tmp.e[i]. */
! 283: if (l-c > 0) {
! 284: q =0;
! 285: for (i=c; i<l; i++) {
! 286: q += (e[i].D)*(tmp.e[i].x); /* Don't repeat these things. */
! 287: tmp.e[i].D += e[i].D;
! 288: }
! 289: /*printf("l=%d, q=%d\n",l,q);*/
! 290: if (ringp->next == (struct ring *)NULL) {
! 291: tmp.e[0].x += q;
! 292: }else{
! 293: Cmult(a,a,polyToCoeff(cxx(1,0,q,ringp->next),ringp));
! 294: /* x[0]^q */
! 295: }
! 296: }
! 297:
! 298: /* Update tmp.e[i].x */
! 299: for (i=0; i<n; i++) {
! 300: tmp.e[i].x += e[i].x;
! 301: }
! 302:
! 303:
! 304: /* commutative variables */
! 305: for (i=0; i<c; i++) {
! 306: tmp.e[i].D += e[i].D;
! 307: }
! 308:
! 309: mm = newCell(a,monomialCopy(&tmp));
! 310: RList[k]->next = mm;
! 311: RList[k] = RList[k]->next;
! 312: no: ;
! 313: }
! 314: }
! 315:
! 316: /* Note also that mpMult_diff assumes coefficients and Dx commutes each other*/
! 317: POLY mpMult_diff(POLY f,POLY g)
! 318: {
! 319: int k;
! 320: POLY r,temp;
! 321: mt_declare_locals()
! 322: struct multTable *mt;
! 323:
! 324: /* printf("mpMult_diff(%s,%s)\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
! 325:
! 326: if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
! 327: checkRing(f,g);
! 328:
! 329:
! 330: if (isConstant(f)) return(cpMult(f->coeffp,g));
! 331:
! 332: mt = newMultTable(f->coeffp);
! 333: makeTable_differential_difference(f->coeffp,f->m->e,f->m->ringp,mt);
! 334: mt_get_mt(mt);
! 335: /* outputMultTable(mt); */
! 336:
! 337: for (k=0; k<Plist; k++) {
! 338: RList[k] = RListRoot[k];
! 339: RList[k]->next = POLYNULL;
! 340: }
! 341:
! 342: while (g != POLYNULL) {
! 343: monomialMult_differential_difference(f->m->e,g,mt);
! 344: g = g->next;
! 345: }
! 346: r = POLYNULL;
! 347: for (k=0; k<Plist; k++) {
! 348: temp = RListRoot[k]->next;
! 349: r = ppAddv(r,temp);
! 350: }
! 351:
! 352: freeMultTable(mt);
! 353: return(r);
! 354: }
! 355:
! 356:
! 357: static void outputMultTable(struct multTable *mt) {
! 358: int i,j;
! 359: mt_declare_locals()
! 360:
! 361: mt_get_mt(mt);
! 362: printf("Maxv = %d Plist=%d",Maxv,Plist);
! 363: printf("\n------- V ----------------\n");
! 364: for (i=0; i<Maxv; i++) printf("%5d",V[i]);
! 365: printf("\n------- Vtype ----------------\n");
! 366: for (i=0; i<Maxv; i++) printf("%5d",Vtype[i]);
! 367: printf("\n------ DList --------------\n");
! 368: for (i=0; i<Plist; i++) {
! 369: for (j=0; j<Maxv; j++) {
! 370: printf("%5d",DList[I(i,j)]);
! 371: }
! 372: putchar('\n');
! 373: }
! 374: printf("\n--------- MList ------------\n");
! 375: for (i=0; i<Plist; i++) {
! 376: for (j=0; j<Maxv; j++) {
! 377: printf("%5d",MList[I(i,j)]);
! 378: }
! 379: printf(" | e=%5d M=%1d c=%s\n",EList[i],Mark[i],coeffToString(CList[i]));
! 380: }
! 381: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to poly3a.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan