Annotation of OpenXM/src/kan96xx/Kan/poly3a.c, Revision 1.2
1.2 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly3a.c,v 1.1 2002/02/04 07:58:28 takayama Exp $ */
1.1 takayama 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*/
1.2 ! takayama 317: POLY mpMult_difference(POLY f,POLY g)
1.1 takayama 318: {
319: int k;
320: POLY r,temp;
321: mt_declare_locals()
322: struct multTable *mt;
323:
1.2 ! takayama 324: /* printf("mpMult_difference(%s,%s)\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
1.1 takayama 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