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