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