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