Annotation of OpenXM/src/kan96xx/Kan/poly4.c, Revision 1.2
1.2 ! takayama 1: /* $OpenXM$ */
1.1 maekawa 2: #include <stdio.h>
3: #include "datatype.h"
4: #include "stackm.h"
5: #include "extern.h"
6: #include "extern2.h"
7: #include "matrix.h"
8: static void shell(int v[],int n);
9: static int degreeOfPrincipalPart(POLY f);
10: static int degreeOfInitW(POLY f,int w[]);
11:
12:
13: static void shell(v,n)
14: int v[];
15: int n;
16: {
17: int gap,i,j,temp;
18:
19: for (gap = n/2; gap > 0; gap /= 2) {
20: for (i = gap; i<n; i++) {
21: for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
22: temp = v[j];
23: v[j] = v[j+gap];
24: v[j+gap] = temp;
25: }
26: }
27: }
28: }
29:
30:
31: struct matrixOfPOLY *parts(f,v)
32: POLY f;
33: POLY v; /* v must be a single variable, e.g. x */
34: {
35: struct matrixOfPOLY *evPoly;
36: int vi = 0; /* index of v */
37: int vx = 1; /* x --> 1, D--> 0 */
38: int n,evSize,i,k,e;
39: int *ev;
40: struct object *evList;
41: struct object *list;
42: struct object ob;
43: POLY ans;
44: POLY h;
45: extern struct ring *CurrentRingp;
46: POLY ft;
47:
48:
49: if (f ISZERO || v ISZERO) {
50: evPoly = newMatrixOfPOLY(2,1);
51: getMatrixOfPOLY(evPoly,0,0) = ZERO;
52: getMatrixOfPOLY(evPoly,1,0) = ZERO;
53: return(evPoly);
54: }
55: n = v->m->ringp->n;
56: /* get the index of the variable v */
57: for (i=0; i<n; i++) {
58: if (v->m->e[i].x) {
59: vx = 1; vi = i; break;
60: }else if (v->m->e[i].D) {
61: vx = 0; vi = i; break;
62: }
63: }
64: ft = f;
65: /* get the vector of exponents */
66: evList = NULLLIST;
67: while (ft != POLYNULL) {
68: if (vx) {
69: e = ft->m->e[vi].x;
70: }else{
71: e = ft->m->e[vi].D;
72: }
73: ft = ft->next;
74: ob = KpoInteger(e);
75: if (!memberQ(evList,ob)) {
76: list = newList(&ob);
77: evList = vJoin(evList,list);
78: }
79: }
80: /*printf("evList = "); printObjectList(evList);*/
81: evSize = klength(evList);
82: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
83: if (ev == (int *)NULL) errorPoly("No more memory.");
84: for (i=0; i<evSize; i++) {
85: ev[i] = KopInteger(car(evList));
86: evList = cdr(evList);
87: }
88: /* sort ev */
89: shell(ev,evSize);
90:
91: /* get coefficients */
92: evPoly = newMatrixOfPOLY(2,evSize);
93: for (i=0; i<evSize; i++) {
94: ans = ZERO;
95: getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
96: ft = f;
97: while (ft != POLYNULL) {
98: if (vx) {
99: if (ft->m->e[vi].x == ev[i]) {
100: h = newCell(ft->coeffp,monomialCopy(ft->m));
101: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
102: ans = ppAdd(ans,h);
103: }
104: }else{
105: if (ft->m->e[vi].D == ev[i]) {
106: h = newCell(ft->coeffp,monomialCopy(ft->m));
107: dset0(h,vi);
108: ans = ppAdd(ans,h);
109: }
110: }
111: ft = ft->next;
112: }
113: getMatrixOfPOLY(evPoly,1,i) = ans;
114: }
115: return(evPoly);
116: }
117:
118: struct object parts2(f,v)
119: POLY f;
120: POLY v; /* v must be a single variable, e.g. x */
121: {
122: struct matrixOfPOLY *evPoly;
123: int vi = 0; /* index of v */
124: int vx = 1; /* x --> 1, D--> 0 */
125: int n,evSize,i,k,e;
126: int *ev;
127: struct object *evList;
128: struct object *list;
129: struct object ob;
130: POLY ans;
131: POLY h;
132: POLY ft;
133: struct object ob1,ob2,rob;
134:
135:
136: if (f ISZERO || v ISZERO) {
137: evPoly = newMatrixOfPOLY(2,1);
138: getMatrixOfPOLY(evPoly,0,0) = ZERO;
139: getMatrixOfPOLY(evPoly,1,0) = ZERO;
140: rob = newObjectArray(2);
141: ob1 = newObjectArray(1);
142: ob2 = newObjectArray(1);
143: putoa(ob1,0,KpoInteger(0));
144: putoa(ob2,0,KpoPOLY(POLYNULL));
145: putoa(rob,0,ob1); putoa(rob,1,ob2);
146: return(rob);
147: }
148: n = v->m->ringp->n;
149: /* get the index of the variable v */
150: for (i=0; i<n; i++) {
151: if (v->m->e[i].x) {
152: vx = 1; vi = i; break;
153: }else if (v->m->e[i].D) {
154: vx = 0; vi = i; break;
155: }
156: }
157: ft = f;
158: /* get the vector of exponents */
159: evList = NULLLIST;
160: while (ft != POLYNULL) {
161: if (vx) {
162: e = ft->m->e[vi].x;
163: }else{
164: e = ft->m->e[vi].D;
165: }
166: ft = ft->next;
167: ob = KpoInteger(e);
168: if (!memberQ(evList,ob)) {
169: list = newList(&ob);
170: evList = vJoin(evList,list);
171: }
172: }
173: /*printf("evList = "); printObjectList(evList);*/
174: evSize = klength(evList);
175: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
176: if (ev == (int *)NULL) errorPoly("No more memory.");
177: for (i=0; i<evSize; i++) {
178: ev[i] = KopInteger(car(evList));
179: evList = cdr(evList);
180: }
181: /* sort ev */
182: shell(ev,evSize);
183:
184: /* get coefficients */
185: evPoly = newMatrixOfPOLY(2,evSize);
186: for (i=0; i<evSize; i++) {
187: ans = ZERO;
188: /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
189: getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
190: ft = f;
191: while (ft != POLYNULL) {
192: if (vx) {
193: if (ft->m->e[vi].x == ev[i]) {
194: h = newCell(ft->coeffp,monomialCopy(ft->m));
195: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
196: ans = ppAdd(ans,h);
197: }
198: }else{
199: if (ft->m->e[vi].D == ev[i]) {
200: h = newCell(ft->coeffp,monomialCopy(ft->m));
201: dset0(h,vi);
202: ans = ppAdd(ans,h);
203: }
204: }
205: ft = ft->next;
206: }
207: getMatrixOfPOLY(evPoly,1,i) = ans;
208: }
209: rob = newObjectArray(2);
210: ob1 = newObjectArray(evSize);
211: ob2 = newObjectArray(evSize);
212: for (i=0; i<evSize; i++) {
213: putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
214: putoa(ob1,i,KpoInteger(ev[i]));
215: }
216: putoa(rob,0,ob1); putoa(rob,1,ob2);
217: return(rob);
218: }
219:
220: int pDegreeWrtV(f,v)
221: POLY f;
222: POLY v;
223: {
224: int vx = 1;
225: int vi = 0;
226: int i,n;
227: int ans;
228: if (f ISZERO || v ISZERO) return(0);
229: n = f->m->ringp->n;
230: for (i=0; i<n; i++) {
231: if (v->m->e[i].x) {
232: vx = 1; vi = i;
233: break;
234: }else if (v->m->e[i].D) {
235: vx = 0; vi = i;
236: break;
237: }
238: }
239: if (vx) {
240: ans = f->m->e[vi].x;
241: }else{
242: ans = f->m->e[vi].D;
243: }
244: f = f->next;
245: while (f != POLYNULL) {
246: if (vx) {
247: if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
248: }else{
249: if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
250: }
251: f = f->next;
252: }
253: return(ans);
254: }
255:
256: int containVectorVariable(POLY f)
257: {
258: MONOMIAL tf;
259: static int nn,mm,ll,cc,n,m,l,c;
260: static struct ring *cr = (struct ring *)NULL;
261: int i;
262:
263: if (f ISZERO) return(0);
264: tf = f->m;
265: if (tf->ringp != cr) {
266: n = tf->ringp->n;
267: m = tf->ringp->m;
268: l = tf->ringp->l;
269: c = tf->ringp->c;
270: nn = tf->ringp->nn;
271: mm = tf->ringp->mm;
272: ll = tf->ringp->ll;
273: cc = tf->ringp->cc;
274: cr = tf->ringp;
275: }
276:
277: while (f != POLYNULL) {
278: tf = f->m;
279: for (i=cc; i<c; i++) {
280: if ( tf->e[i].x ) return(1);
281: if ( tf->e[i].D ) return(1);
282: }
283: for (i=ll; i<l; i++) {
284: if (tf->e[i].x) return(1);
285: if (tf->e[i].D) return(1);
286: }
287: for (i=mm; i<m; i++) {
288: if (tf->e[i].x) return(1);
289: if (tf->e[i].D) return(1);
290: }
291: for (i=nn; i<n; i++) {
292: if (tf->e[i].x) return(1);
293: if (tf->e[i].D) return(1);
294: }
295: f = f->next;
296: }
297: return(0);
298:
299: }
300:
301: POLY homogenize(f)
302: POLY f;
303: /* homogenize by using (*grade)(f) */
304: {
305: POLY t;
306: int maxg;
307: int flag,d;
308:
309: if (f == ZERO) return(f);
310: t = f; maxg = (*grade)(f); flag = 0;
311: while (t != POLYNULL) {
312: d = (*grade)(t);
313: if (d != maxg) flag = 1;
314: if (d > maxg) {
315: maxg = d;
316: }
317: t = t->next;
318: }
319: if (flag == 0) return(f);
320:
321: f = pmCopy(f); /* You can rewrite the monomial parts */
322: t = f;
323: while (t != POLYNULL) {
324: d = (*grade)(t);
325: if (d != maxg) {
326: t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
327: }
328: t = t->next;
329: }
330: return(f);
331: }
332:
333: int isHomogenized(f)
334: POLY f;
335: {
336: POLY t;
337: extern int Homogenize_vec;
338: int maxg;
339: if (!Homogenize_vec) return(isHomogenized_vec(f));
340: if (f == ZERO) return(1);
341: maxg = (*grade)(f);
342: t = f;
343: while (t != POLYNULL) {
344: if (maxg != (*grade)(t)) return(0);
345: t = t->next;
346: }
347: return(1);
348: }
349:
350: int isHomogenized_vec(f)
351: POLY f;
352: {
353: /* This is not efficient version. *grade should be grade_module1v(). */
354: POLY t;
355: int ggg;
356: if (f == ZERO) return(1);
357: while (f != POLYNULL) {
358: t = f;
359: ggg = (*grade)(f);
360: while (t != POLYNULL) {
361: if ((*isSameComponent)(f,t)) {
362: if (ggg != (*grade)(t)) return(0);
363: }
364: t = t->next;
365: }
366: f = f->next;
367: }
368: return(1);
369: }
370:
371:
372: static int degreeOfPrincipalPart(f)
373: POLY f;
374: {
375: int n,i,dd;
376: if (f ISZERO) return(0);
377: n = f->m->ringp->n; dd = 0;
378: /* D[0] is homogenization var */
379: for (i=1; i<n; i++) {
380: dd += f->m->e[i].D;
381: }
382: return(dd);
383: }
384:
385: POLY POLYToPrincipalPart(f)
386: POLY f;
387: {
388: POLY node;
389: struct listPoly nod;
390: POLY h;
391: POLY g;
392: int maxd = -20000; /* very big negative number */
393: int dd;
394: node = &nod; node->next = POLYNULL; h = node;
395:
396: g = pCopy(f); /* shallow copy */
397: while (!(f ISZERO)) {
398: dd = degreeOfPrincipalPart(f);
399: if (dd > maxd) maxd = dd;
400: f = f->next;
401: }
402: while (!(g ISZERO)) {
403: dd = degreeOfPrincipalPart(g);
404: if (dd == maxd) {
405: h->next = g;
406: h = h->next;
407: }
408: g = g->next;
409: }
410: h->next = POLYNULL;
411: return(node->next);
412: }
413:
414: static int degreeOfInitW(f,w)
415: POLY f;
416: int w[];
417: {
418: int n,i,dd;
419: if (f ISZERO) {
420: errorPoly("degreeOfInitW(0,w) ");
421: }
422: n = f->m->ringp->n; dd = 0;
423: for (i=0; i<n; i++) {
424: dd += (f->m->e[i].D)*w[n+i];
425: dd += (f->m->e[i].x)*w[i];
426: }
427: return(dd);
428: }
429:
430: POLY POLYToInitW(f,w)
431: POLY f;
432: int w[]; /* weight vector */
433: {
434: POLY node;
435: struct listPoly nod;
436: POLY h;
437: POLY g;
438: int maxd;
439: int dd;
440: node = &nod; node->next = POLYNULL; h = node;
441:
442: if (f ISZERO) return(f);
443: maxd = degreeOfInitW(f,w);
444: g = pCopy(f); /* shallow copy */
445: while (!(f ISZERO)) {
446: dd = degreeOfInitW(f,w);
447: if (dd > maxd) maxd = dd;
448: f = f->next;
449: }
450: while (!(g ISZERO)) {
451: dd = degreeOfInitW(g,w);
452: if (dd == maxd) {
453: h->next = g;
454: h = h->next;
455: }
456: g = g->next;
457: }
458: h->next = POLYNULL;
459: return(node->next);
460: }
461:
462:
463: /*
464: 1.The substitution "ringp->multiplication = ...." is allowed only in
465: KsetUpRing(), so the check in KswitchFunction is not necessary.
466: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
467: 3.If Schreyer = 1, then the system always generates a new ring.
468: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
469: 5.When mmLarger == tower (in tower.sm1, tower-sugar.sm1), we do
470: as follows with our own risk.
471: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
472: */
473: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
474: {
475: struct ring *rrr;
476: int i,j,k;
477: int a=0;
478: for (k=0; k<rp; k++) {
479: rrr = rstack[k];
480: if (rrr->p != newRingp->p) { a=1; goto bbb ; }
481: if (rrr->n != newRingp->n) { a=2; goto bbb ; }
482: if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
483: if (rrr->m != newRingp->m) { a=4; goto bbb ; }
484: if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
485: if (rrr->l != newRingp->l) { a=6; goto bbb ; }
486: if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
487: if (rrr->c != newRingp->c) { a=8; goto bbb ; }
488: if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
489: for (i=0; i<rrr->n; i++) {
490: if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
491: if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
492: }
493: if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
494: for (i=0; i<rrr->orderMatrixSize; i++) {
495: for (j=0; j<2*(rrr->n); j++) {
496: if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
497: { a=13; goto bbb ; }
498: }
499: }
500: if (rrr->next != newRingp->next) { a=14; goto bbb ; }
501: if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
502: /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
503: if (newRingp->schreyer == 1) { a=16; goto bbb; }
504: /* The following fields are ignored.
505: void *gbListTower;
506: int *outputOrder;
507: char *name;
508: */
509: /* All tests are passed. */
510: return(k);
511: bbb: ;
512: /* for debugging. */
513: /* fprintf(stderr," reason=%d, ",a); */
514: }
515: return(-1);
516: }
517:
518:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to poly4.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan