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