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