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