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