Annotation of OpenXM/src/kan96xx/Kan/poly4.c, Revision 1.16
1.16 ! ohara 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly4.c,v 1.15 2005/06/16 05:07:23 takayama Exp $ */
1.1 maekawa 2: #include <stdio.h>
1.16 ! ohara 3: #include <string.h>
1.1 maekawa 4: #include "datatype.h"
5: #include "stackm.h"
6: #include "extern.h"
7: #include "extern2.h"
8: #include "matrix.h"
9: static void shell(int v[],int n);
10: static int degreeOfPrincipalPart(POLY f);
11: static int degreeOfInitW(POLY f,int w[]);
1.10 takayama 12: static int degreeOfInitWS(POLY f,int w[],int s[]);
1.13 takayama 13: static int dDegree(POLY f);
14: static POLY dHomogenize(POLY f);
1.1 maekawa 15:
16: static void shell(v,n)
1.3 takayama 17: int v[];
18: int n;
1.1 maekawa 19: {
20: int gap,i,j,temp;
21:
22: for (gap = n/2; gap > 0; gap /= 2) {
23: for (i = gap; i<n; i++) {
24: for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
1.3 takayama 25: temp = v[j];
26: v[j] = v[j+gap];
27: v[j+gap] = temp;
1.1 maekawa 28: }
29: }
30: }
31: }
32:
33:
34: struct matrixOfPOLY *parts(f,v)
1.3 takayama 35: POLY f;
36: POLY v; /* v must be a single variable, e.g. x */
1.1 maekawa 37: {
38: struct matrixOfPOLY *evPoly;
39: int vi = 0; /* index of v */
40: int vx = 1; /* x --> 1, D--> 0 */
41: int n,evSize,i,k,e;
42: int *ev;
43: struct object *evList;
44: struct object *list;
1.15 takayama 45: struct object ob = OINIT;
1.1 maekawa 46: POLY ans;
47: POLY h;
48: extern struct ring *CurrentRingp;
49: POLY ft;
50:
51:
52: if (f ISZERO || v ISZERO) {
53: evPoly = newMatrixOfPOLY(2,1);
54: getMatrixOfPOLY(evPoly,0,0) = ZERO;
55: getMatrixOfPOLY(evPoly,1,0) = ZERO;
56: return(evPoly);
57: }
58: n = v->m->ringp->n;
59: /* get the index of the variable v */
60: for (i=0; i<n; i++) {
61: if (v->m->e[i].x) {
62: vx = 1; vi = i; break;
63: }else if (v->m->e[i].D) {
64: vx = 0; vi = i; break;
65: }
66: }
67: ft = f;
68: /* get the vector of exponents */
69: evList = NULLLIST;
70: while (ft != POLYNULL) {
71: if (vx) {
72: e = ft->m->e[vi].x;
73: }else{
74: e = ft->m->e[vi].D;
75: }
76: ft = ft->next;
77: ob = KpoInteger(e);
78: if (!memberQ(evList,ob)) {
79: list = newList(&ob);
80: evList = vJoin(evList,list);
81: }
82: }
83: /*printf("evList = "); printObjectList(evList);*/
84: evSize = klength(evList);
85: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
86: if (ev == (int *)NULL) errorPoly("No more memory.");
87: for (i=0; i<evSize; i++) {
88: ev[i] = KopInteger(car(evList));
89: evList = cdr(evList);
90: }
91: /* sort ev */
92: shell(ev,evSize);
93:
94: /* get coefficients */
95: evPoly = newMatrixOfPOLY(2,evSize);
96: for (i=0; i<evSize; i++) {
97: ans = ZERO;
98: getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
99: ft = f;
100: while (ft != POLYNULL) {
101: if (vx) {
1.3 takayama 102: if (ft->m->e[vi].x == ev[i]) {
103: h = newCell(ft->coeffp,monomialCopy(ft->m));
104: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
105: ans = ppAdd(ans,h);
106: }
1.1 maekawa 107: }else{
1.3 takayama 108: if (ft->m->e[vi].D == ev[i]) {
109: h = newCell(ft->coeffp,monomialCopy(ft->m));
110: dset0(h,vi);
111: ans = ppAdd(ans,h);
112: }
1.1 maekawa 113: }
114: ft = ft->next;
115: }
116: getMatrixOfPOLY(evPoly,1,i) = ans;
117: }
118: return(evPoly);
119: }
1.3 takayama 120:
1.1 maekawa 121: struct object parts2(f,v)
1.3 takayama 122: POLY f;
123: POLY v; /* v must be a single variable, e.g. x */
1.1 maekawa 124: {
125: struct matrixOfPOLY *evPoly;
126: int vi = 0; /* index of v */
127: int vx = 1; /* x --> 1, D--> 0 */
128: int n,evSize,i,k,e;
129: int *ev;
130: struct object *evList;
131: struct object *list;
1.15 takayama 132: struct object ob = OINIT;
1.1 maekawa 133: POLY ans;
134: POLY h;
135: POLY ft;
1.15 takayama 136: struct object ob1 = OINIT;
137: struct object ob2 = OINIT;
138: struct object rob = OINIT;
1.1 maekawa 139:
140:
141: if (f ISZERO || v ISZERO) {
142: evPoly = newMatrixOfPOLY(2,1);
143: getMatrixOfPOLY(evPoly,0,0) = ZERO;
144: getMatrixOfPOLY(evPoly,1,0) = ZERO;
145: rob = newObjectArray(2);
146: ob1 = newObjectArray(1);
147: ob2 = newObjectArray(1);
148: putoa(ob1,0,KpoInteger(0));
149: putoa(ob2,0,KpoPOLY(POLYNULL));
150: putoa(rob,0,ob1); putoa(rob,1,ob2);
151: return(rob);
152: }
153: n = v->m->ringp->n;
154: /* get the index of the variable v */
155: for (i=0; i<n; i++) {
156: if (v->m->e[i].x) {
157: vx = 1; vi = i; break;
158: }else if (v->m->e[i].D) {
159: vx = 0; vi = i; break;
160: }
161: }
162: ft = f;
163: /* get the vector of exponents */
164: evList = NULLLIST;
165: while (ft != POLYNULL) {
166: if (vx) {
167: e = ft->m->e[vi].x;
168: }else{
169: e = ft->m->e[vi].D;
170: }
171: ft = ft->next;
172: ob = KpoInteger(e);
173: if (!memberQ(evList,ob)) {
174: list = newList(&ob);
175: evList = vJoin(evList,list);
176: }
177: }
178: /*printf("evList = "); printObjectList(evList);*/
179: evSize = klength(evList);
180: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
181: if (ev == (int *)NULL) errorPoly("No more memory.");
182: for (i=0; i<evSize; i++) {
183: ev[i] = KopInteger(car(evList));
184: evList = cdr(evList);
185: }
186: /* sort ev */
187: shell(ev,evSize);
188:
189: /* get coefficients */
190: evPoly = newMatrixOfPOLY(2,evSize);
191: for (i=0; i<evSize; i++) {
192: ans = ZERO;
193: /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
194: getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
195: ft = f;
196: while (ft != POLYNULL) {
197: if (vx) {
1.3 takayama 198: if (ft->m->e[vi].x == ev[i]) {
199: h = newCell(ft->coeffp,monomialCopy(ft->m));
200: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
201: ans = ppAdd(ans,h);
202: }
1.1 maekawa 203: }else{
1.3 takayama 204: if (ft->m->e[vi].D == ev[i]) {
205: h = newCell(ft->coeffp,monomialCopy(ft->m));
206: dset0(h,vi);
207: ans = ppAdd(ans,h);
208: }
1.1 maekawa 209: }
210: ft = ft->next;
211: }
212: getMatrixOfPOLY(evPoly,1,i) = ans;
213: }
214: rob = newObjectArray(2);
215: ob1 = newObjectArray(evSize);
216: ob2 = newObjectArray(evSize);
217: for (i=0; i<evSize; i++) {
218: putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
219: putoa(ob1,i,KpoInteger(ev[i]));
220: }
221: putoa(rob,0,ob1); putoa(rob,1,ob2);
222: return(rob);
223: }
1.3 takayama 224:
1.1 maekawa 225: int pDegreeWrtV(f,v)
1.3 takayama 226: POLY f;
227: POLY v;
1.1 maekawa 228: {
229: int vx = 1;
230: int vi = 0;
231: int i,n;
232: int ans;
233: if (f ISZERO || v ISZERO) return(0);
234: n = f->m->ringp->n;
235: for (i=0; i<n; i++) {
236: if (v->m->e[i].x) {
237: vx = 1; vi = i;
238: break;
239: }else if (v->m->e[i].D) {
240: vx = 0; vi = i;
241: break;
242: }
243: }
244: if (vx) {
245: ans = f->m->e[vi].x;
246: }else{
247: ans = f->m->e[vi].D;
248: }
249: f = f->next;
250: while (f != POLYNULL) {
251: if (vx) {
252: if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
253: }else{
254: if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
255: }
256: f = f->next;
257: }
258: return(ans);
259: }
260:
261: int containVectorVariable(POLY f)
262: {
263: MONOMIAL tf;
264: static int nn,mm,ll,cc,n,m,l,c;
265: static struct ring *cr = (struct ring *)NULL;
266: int i;
267:
268: if (f ISZERO) return(0);
269: tf = f->m;
270: if (tf->ringp != cr) {
271: n = tf->ringp->n;
272: m = tf->ringp->m;
273: l = tf->ringp->l;
274: c = tf->ringp->c;
275: nn = tf->ringp->nn;
276: mm = tf->ringp->mm;
277: ll = tf->ringp->ll;
278: cc = tf->ringp->cc;
279: cr = tf->ringp;
280: }
281:
282: while (f != POLYNULL) {
283: tf = f->m;
284: for (i=cc; i<c; i++) {
285: if ( tf->e[i].x ) return(1);
286: if ( tf->e[i].D ) return(1);
287: }
288: for (i=ll; i<l; i++) {
289: if (tf->e[i].x) return(1);
290: if (tf->e[i].D) return(1);
291: }
292: for (i=mm; i<m; i++) {
293: if (tf->e[i].x) return(1);
294: if (tf->e[i].D) return(1);
295: }
296: for (i=nn; i<n; i++) {
297: if (tf->e[i].x) return(1);
298: if (tf->e[i].D) return(1);
299: }
300: f = f->next;
301: }
302: return(0);
303:
304: }
305:
306: POLY homogenize(f)
1.3 takayama 307: POLY f;
308: /* homogenize by using (*grade)(f) */
1.1 maekawa 309: {
310: POLY t;
311: int maxg;
312: int flag,d;
1.13 takayama 313: extern int Homogenize;
1.1 maekawa 314:
315: if (f == ZERO) return(f);
1.13 takayama 316: if (Homogenize == 3) { /* double homogenization Dx x = x Dx + h H */
317: return dHomogenize(f);
318: }
1.1 maekawa 319: t = f; maxg = (*grade)(f); flag = 0;
320: while (t != POLYNULL) {
321: d = (*grade)(t);
322: if (d != maxg) flag = 1;
323: if (d > maxg) {
324: maxg = d;
325: }
326: t = t->next;
327: }
328: if (flag == 0) return(f);
329:
330: f = pmCopy(f); /* You can rewrite the monomial parts */
331: t = f;
332: while (t != POLYNULL) {
333: d = (*grade)(t);
334: if (d != maxg) {
335: t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
336: }
337: t = t->next;
338: }
339: return(f);
340: }
341:
342: int isHomogenized(f)
1.3 takayama 343: POLY f;
1.1 maekawa 344: {
345: POLY t;
346: extern int Homogenize_vec;
347: int maxg;
348: if (!Homogenize_vec) return(isHomogenized_vec(f));
349: if (f == ZERO) return(1);
1.4 takayama 350: if (f->m->ringp->weightedHomogenization) {
351: return 1; /* BUG: do not chech in case of one-zero homogenization */
352: }
1.1 maekawa 353: maxg = (*grade)(f);
354: t = f;
355: while (t != POLYNULL) {
356: if (maxg != (*grade)(t)) return(0);
357: t = t->next;
358: }
359: return(1);
360: }
361:
362: int isHomogenized_vec(f)
1.3 takayama 363: POLY f;
1.1 maekawa 364: {
1.3 takayama 365: /* This is not efficient version. *grade should be grade_module1v(). */
1.1 maekawa 366: POLY t;
367: int ggg;
368: if (f == ZERO) return(1);
1.4 takayama 369: if (f->m->ringp->weightedHomogenization) {
370: return 1; /* BUG: do not chech in case of one-zero homogenization */
371: }
1.1 maekawa 372: while (f != POLYNULL) {
373: t = f;
374: ggg = (*grade)(f);
375: while (t != POLYNULL) {
376: if ((*isSameComponent)(f,t)) {
1.3 takayama 377: if (ggg != (*grade)(t)) return(0);
1.1 maekawa 378: }
379: t = t->next;
380: }
381: f = f->next;
382: }
383: return(1);
384: }
385:
1.13 takayama 386: static POLY dHomogenize(f)
387: POLY f;
388: {
389: POLY t;
390: int maxg, maxdg;
391: int flag,d,dd,neg;
392:
393: if (f == ZERO) return(f);
1.14 takayama 394:
395: t = f;
396: maxg = (*grade)(f);
397: while (t != POLYNULL) {
398: dd = (*grade)(t);
399: if (maxg < dd) maxg = dd;
400: t = t->next;
401: }
402: /* fprintf(stderr,"maxg=%d\n",maxg); */
403:
404: t = f;
405: maxdg = dDegree(f);
406: while (t != POLYNULL) {
407: dd = dDegree(t);
408: if (maxdg < dd) maxdg = dd;
409: t = t->next;
410: }
411: /* fprintf(stderr,"maxdg=%d\n",maxdg); */
412:
413: t = f;
414: flag = 0;
1.13 takayama 415: while (t != POLYNULL) {
416: d = (*grade)(t);
417: if (d != maxg) flag = 1;
418: if (d > maxg) {
419: maxg = d;
420: }
421: d = dDegree(f);
422: if (d > maxdg) {
423: maxdg = d;
424: }
425: t = t->next;
426: }
427: if (flag == 0) return(f);
428:
429: t = f; neg = 0;
430: while (t != POLYNULL) {
431: d = (*grade)(t);
432: dd = dDegree(t);
433: if (maxg-d-(maxdg-dd) < neg) {
434: neg = maxg-d-(maxdg-dd);
435: }
436: t = t->next;
437: }
438: neg = -neg;
439:
440: f = pmCopy(f); /* You can rewrite the monomial parts */
441: t = f;
442: while (t != POLYNULL) {
443: d = (*grade)(t);
444: dd = dDegree(t);
445: t->m->e[0].D += maxdg-dd; /* h */
446: t->m->e[0].x += maxg-d-(maxdg-dd)+neg; /* Multiply H */
447: /* example Dx^2+Dx+x */
448: t = t->next;
449: }
450: return(f);
451: }
1.1 maekawa 452:
453: static int degreeOfPrincipalPart(f)
1.3 takayama 454: POLY f;
1.1 maekawa 455: {
456: int n,i,dd;
457: if (f ISZERO) return(0);
458: n = f->m->ringp->n; dd = 0;
459: /* D[0] is homogenization var */
460: for (i=1; i<n; i++) {
1.13 takayama 461: dd += f->m->e[i].D;
462: }
463: return(dd);
464: }
465:
466: static int dDegree(f)
467: POLY f;
468: {
469: int nn,i,dd,m;
470: if (f ISZERO) return(0);
471: nn = f->m->ringp->nn; dd = 0;
472: m = f->m->ringp->m;
473: for (i=m; i<nn; i++) {
1.1 maekawa 474: dd += f->m->e[i].D;
475: }
476: return(dd);
477: }
478:
479: POLY POLYToPrincipalPart(f)
1.3 takayama 480: POLY f;
1.1 maekawa 481: {
482: POLY node;
483: struct listPoly nod;
484: POLY h;
485: POLY g;
486: int maxd = -20000; /* very big negative number */
487: int dd;
488: node = &nod; node->next = POLYNULL; h = node;
489:
490: g = pCopy(f); /* shallow copy */
491: while (!(f ISZERO)) {
492: dd = degreeOfPrincipalPart(f);
493: if (dd > maxd) maxd = dd;
494: f = f->next;
495: }
496: while (!(g ISZERO)) {
497: dd = degreeOfPrincipalPart(g);
498: if (dd == maxd) {
499: h->next = g;
500: h = h->next;
501: }
502: g = g->next;
503: }
504: h->next = POLYNULL;
505: return(node->next);
506: }
507:
508: static int degreeOfInitW(f,w)
1.3 takayama 509: POLY f;
510: int w[];
1.1 maekawa 511: {
512: int n,i,dd;
513: if (f ISZERO) {
514: errorPoly("degreeOfInitW(0,w) ");
515: }
516: n = f->m->ringp->n; dd = 0;
517: for (i=0; i<n; i++) {
518: dd += (f->m->e[i].D)*w[n+i];
519: dd += (f->m->e[i].x)*w[i];
520: }
521: return(dd);
522: }
523:
524: POLY POLYToInitW(f,w)
1.3 takayama 525: POLY f;
526: int w[]; /* weight vector */
1.1 maekawa 527: {
528: POLY h;
529: POLY g;
530: int maxd;
531: int dd;
1.11 takayama 532: h = POLYNULL;
1.1 maekawa 533:
1.11 takayama 534: /*printf("1:%s\n",POLYToString(f,'*',1));*/
1.1 maekawa 535: if (f ISZERO) return(f);
536: maxd = degreeOfInitW(f,w);
1.11 takayama 537: g = f;
1.1 maekawa 538: while (!(f ISZERO)) {
539: dd = degreeOfInitW(f,w);
540: if (dd > maxd) maxd = dd;
541: f = f->next;
542: }
543: while (!(g ISZERO)) {
544: dd = degreeOfInitW(g,w);
545: if (dd == maxd) {
1.11 takayama 546: h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.1 maekawa 547: }
548: g = g->next;
549: }
1.11 takayama 550: /*printf("2:%s\n",POLYToString(h,'*',1));*/
551: return(h);
1.10 takayama 552: }
553:
554: static int degreeOfInitWS(f,w,s)
555: POLY f;
556: int w[];
557: int s[];
558: {
559: int n,i,dd;
560: if (f ISZERO) {
561: errorPoly("degreeOfInitWS(0,w) ");
562: }
1.11 takayama 563: if (s == (int *) NULL) return degreeOfInitW(f,w);
1.10 takayama 564: n = f->m->ringp->n; dd = 0;
565: for (i=0; i<n-1; i++) {
566: dd += (f->m->e[i].D)*w[n+i];
567: dd += (f->m->e[i].x)*w[i];
568: }
569: dd += s[(f->m->e[n-1].x)];
570: return(dd);
571: }
572:
573: POLY POLYToInitWS(f,w,s)
574: POLY f;
575: int w[]; /* weight vector */
576: int s[]; /* shift vector */
577: {
578: POLY h;
579: POLY g;
580: int maxd;
581: int dd;
1.11 takayama 582: h = POLYNULL;
1.10 takayama 583:
1.11 takayama 584: /*printf("1s:%s\n",POLYToString(f,'*',1));*/
1.10 takayama 585: if (f ISZERO) return(f);
586: maxd = degreeOfInitWS(f,w,s);
1.11 takayama 587: g = f;
1.10 takayama 588: while (!(f ISZERO)) {
589: dd = degreeOfInitWS(f,w,s);
590: if (dd > maxd) maxd = dd;
591: f = f->next;
592: }
593: while (!(g ISZERO)) {
594: dd = degreeOfInitWS(g,w,s);
595: if (dd == maxd) {
1.11 takayama 596: h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.10 takayama 597: }
598: g = g->next;
599: }
1.11 takayama 600: /*printf("2s:%s\n",POLYToString(h,'*',1));*/
601: return(h);
1.10 takayama 602: }
603:
604: int ordWsAll(f,w,s)
605: POLY f;
606: int w[]; /* weight vector */
607: int s[]; /* shift vector */
608: {
609: int maxd;
610: int dd;
611:
612: if (f ISZERO) errorPoly("ordWsAll(0,w,s) ");
613: maxd = degreeOfInitWS(f,w,s);
614: while (!(f ISZERO)) {
615: dd = degreeOfInitWS(f,w,s);
616: if (dd > maxd) maxd = dd;
617: f = f->next;
618: }
619: return maxd;
1.1 maekawa 620: }
621:
622:
623: /*
624: 1.The substitution "ringp->multiplication = ...." is allowed only in
625: KsetUpRing(), so the check in KswitchFunction is not necessary.
626: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
627: 3.If Schreyer = 1, then the system always generates a new ring.
628: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
629: 5.When mmLarger == tower (in tower.sm1, tower-sugar.sm1), we do
630: as follows with our own risk.
631: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
632: */
633: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
634: {
635: struct ring *rrr;
636: int i,j,k;
637: int a=0;
638: for (k=0; k<rp; k++) {
639: rrr = rstack[k];
640: if (rrr->p != newRingp->p) { a=1; goto bbb ; }
641: if (rrr->n != newRingp->n) { a=2; goto bbb ; }
642: if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
643: if (rrr->m != newRingp->m) { a=4; goto bbb ; }
644: if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
645: if (rrr->l != newRingp->l) { a=6; goto bbb ; }
646: if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
647: if (rrr->c != newRingp->c) { a=8; goto bbb ; }
648: if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
649: for (i=0; i<rrr->n; i++) {
650: if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
651: if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
652: }
653: if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
654: for (i=0; i<rrr->orderMatrixSize; i++) {
655: for (j=0; j<2*(rrr->n); j++) {
1.3 takayama 656: if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
657: { a=13; goto bbb ; }
1.1 maekawa 658: }
659: }
660: if (rrr->next != newRingp->next) { a=14; goto bbb ; }
661: if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
662: /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
663: if (newRingp->schreyer == 1) { a=16; goto bbb; }
1.4 takayama 664: if (rrr->weightedHomogenization != newRingp->weightedHomogenization) { a=16; goto bbb; }
1.7 takayama 665: if (rrr->degreeShiftSize != newRingp->degreeShiftSize) {
666: a = 17; goto bbb;
667: }
668: if (rrr->degreeShiftN != newRingp->degreeShiftN) {
669: a = 17; goto bbb;
670: }
671: for (i=0; i < rrr->degreeShiftSize; i++) {
672: for (j=0; j< rrr->degreeShiftN; j++) {
673: if (rrr->degreeShift[i*(rrr->degreeShiftN)+j] !=
674: newRingp->degreeShift[i*(rrr->degreeShiftN)+j]) {
675: a = 17; goto bbb;
676: }
677: }
678: }
679:
1.1 maekawa 680: /* The following fields are ignored.
681: void *gbListTower;
682: int *outputOrder;
683: char *name;
684: */
685: /* All tests are passed. */
686: return(k);
687: bbb: ;
688: /* for debugging. */
689: /* fprintf(stderr," reason=%d, ",a); */
690: }
691: return(-1);
692: }
1.6 takayama 693:
694: /* s->1 */
695: POLY goDeHomogenizeS(POLY f) {
1.9 takayama 696: POLY lRule[1];
697: POLY rRule[1];
698: struct ring *rp;
699: POLY ans;
700: /* printf("1:[%s]\n",POLYToString(f,'*',1)); */
701: if (f == POLYNULL) return f;
702: rp = f->m->ringp;
703: if (rp->next == NULL) {
704: lRule[0] = cxx(1,0,1,rp);
705: rRule[0] = cxx(1,0,0,rp);
706: ans=replace(f,lRule,rRule,1);
707: }else{
708: struct coeff *cp;
709: POLY t;
710: POLY nc;
711: ans = POLYNULL;
712: while (f != POLYNULL) {
713: cp = f->coeffp;
714: if (cp->tag == POLY_COEFF) {
715: t = goDeHomogenizeS((cp->val).f);
1.11 takayama 716: nc = newCell(polyToCoeff(t,f->m->ringp),monomialCopy(f->m));
1.9 takayama 717: ans = ppAddv(ans,nc);
718: f = f->next;
719: }else{
720: ans = f; break;
721: }
722: }
723: }
724: /* printf("2:[%s]\n",POLYToString(ans,'*',1)); */
725: return ans;
726: }
727:
728: POLY goDeHomogenizeS_buggy(POLY f) {
1.6 takayama 729: POLY node;
730: POLY lastf;
731: struct listPoly nod;
732: POLY h;
733: POLY tf;
734: int gt,first;
735:
1.9 takayama 736: printf("1:[%s]\n",POLYToString(f,'*',1));
1.6 takayama 737: if (f == POLYNULL) return(POLYNULL);
738: node = &nod;
739: node->next = POLYNULL;
740: lastf = POLYNULL;
741: first = 1;
742: while (f != POLYNULL) {
743: tf = newCell(f->coeffp,monomialCopy(f->m));
744: tf->m->e[0].x = 0; /* H, s variable in the G-O paper. */
745: if (first) {
746: node->next = tf;
747: lastf = tf;
748: first = 0;
749: }else{
750: gt = (*mmLarger)(lastf,tf);
751: if (gt == 1) {
752: lastf->next = tf;
753: lastf = tf;
754: }else{
755: h = node->next;
756: h = ppAddv(h,tf);
757: node->next = h;
758: lastf = h;
759: while (lastf->next != POLYNULL) {
760: lastf = lastf->next;
761: }
762: }
763: }
764: f = f->next;
765: }
1.9 takayama 766: printf("2:[%s]\n",POLYToString(node->next,'*',1));
1.6 takayama 767: return (node->next);
768: }
769:
1.5 takayama 770: /* Granger-Oaku's homogenization for the ecart tangent cone.
771: Note: 2003.07.10.
772: ds[] is the degree shift.
773: ei ( element index ). If it is < 0, then e[n-1]->x will be used,
774: else ei is used.
1.6 takayama 775: if onlyS is set to 1, then input is assumed to be (u,v)-h-homogeneous.
1.5 takayama 776: */
1.6 takayama 777: POLY goHomogenize(POLY f,int u[],int v[],int ds[],int dssize,int ei,int onlyS)
1.5 takayama 778: {
779: POLY node;
780: POLY lastf;
781: struct listPoly nod;
782: POLY h;
783: POLY tf;
784: int gt,first,m,mp,t,tp,dsIdx,message;
1.7 takayama 785: struct ring *rp;
1.5 takayama 786:
787: message = 1;
788: if (f == POLYNULL) return(POLYNULL);
1.7 takayama 789: rp = f->m->ringp;
1.12 takayama 790: /*
1.7 takayama 791: if ((rp->degreeShiftSize == 0) && (dssize > 0)) {
792: warningPoly("You are trying to homogenize a polynomial with degree shift. However, the polynomial belongs to the ring without degreeShift option. It may cause a trouble in comparison in free module.\n");
793: }
1.12 takayama 794: */
1.5 takayama 795: node = &nod;
796: node->next = POLYNULL;
797: lastf = POLYNULL;
798: first = 1;
799: while (f != POLYNULL) {
800: if (first) {
801: t = m = dGrade1(f);
802: tp = mp = uvGrade1(f,u,v,ds,dssize,ei);
803: }else{
804: t = dGrade1(f);
805: tp = uvGrade1(f,u,v,ds,dssize,ei);
806: if (t > m) m = t;
807: if (tp < mp) mp = tp;
808: }
809:
810: tf = newCell(f->coeffp,monomialCopy(f->m));
811: /* Automatic dehomogenize. Not += */
812: if (message && ((tf->m->e[0].D != 0) || (tf->m->e[0].x != 0))) {
813: /*go-debug fprintf(stderr,"Automatic dehomogenize and homogenize.\n"); */
814: message = 0;
815: }
1.6 takayama 816: if (!onlyS) {
817: tf->m->e[0].D = -t; /* h */
818: }
1.5 takayama 819: tf->m->e[0].x = tp; /* H, s variable in the G-O paper. */
820: /*go-debug printf("t(h)=%d, tp(uv+ds)=%d\n",t,tp); */
821: if (first) {
822: node->next = tf;
823: lastf = tf;
824: first = 0;
825: }else{
826: gt = (*mmLarger)(lastf,tf);
827: if (gt == 1) {
828: lastf->next = tf;
829: lastf = tf;
830: }else{
831: /*go-debug printf("?\n"); */
832: h = node->next;
833: h = ppAddv(h,tf);
834: node->next = h;
835: lastf = h;
836: while (lastf->next != POLYNULL) {
837: lastf = lastf->next;
838: }
839: }
840: }
841: f = f->next;
842: }
843: h = node->next;
844: /*go-debug printf("m=%d, mp=%d\n",m,mp); */
845: while (h != POLYNULL) {
1.12 takayama 846: /*go-debug printf("Old: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x); */
1.6 takayama 847: if (!onlyS) h->m->e[0].D += m; /* h */
1.5 takayama 848: h->m->e[0].x += -mp; /* H, s*/
849: /*go-debug printf("New: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x); */
850: h = h->next;
851: }
852: return (node->next);
853: }
854:
855: /* u[] = -1, v[] = 1 */
1.6 takayama 856: POLY goHomogenize11(POLY f,int ds[],int dssize,int ei,int onlyS)
1.5 takayama 857: {
858: int r;
859: int i,t,n,m,nn;
860: MONOMIAL tf;
861: static int *u;
862: static int *v;
863: static struct ring *cr = (struct ring *)NULL;
864:
865: if (f == POLYNULL) return POLYNULL;
866:
867: tf = f->m;
868: if (tf->ringp != cr) {
869: n = tf->ringp->n;
870: m = tf->ringp->m;
871: nn = tf->ringp->nn;
872: cr = tf->ringp;
873: u = (int *)sGC_malloc(sizeof(int)*n);
874: v = (int *)sGC_malloc(sizeof(int)*n);
875: for (i=0; i<n; i++) u[i]=v[i]=0;
876: for (i=m; i<nn; i++) {
877: u[i] = -1; v[i] = 1;
878: }
879: }
1.6 takayama 880: return(goHomogenize(f,u,v,ds,dssize,ei,onlyS));
1.5 takayama 881: }
882:
1.6 takayama 883: POLY goHomogenize_dsIdx(POLY f,int u[],int v[],int dsIdx,int ei,int onlyS)
1.5 takayama 884: {
885: if (f == POLYNULL) return POLYNULL;
886: }
1.6 takayama 887: POLY goHomogenize11_dsIdx(POLY f,int ds[],int dsIdx,int ei,int onlyS)
1.5 takayama 888: {
889: if (f == POLYNULL) return POLYNULL;
1.8 takayama 890: }
891:
892: /* cf. KsetUpRing() in kanExport0.c */
893: struct ring *newRingOverFp(struct ring *rp,int p) {
894: struct ring *newRingp;
895: char *ringName = NULL;
896: char pstr[64];
897: sprintf(pstr,"_%d",p);
898: ringName = (char *)sGC_malloc(128);
899: newRingp = (struct ring *)sGC_malloc(sizeof(struct ring));
900: if (newRingp == NULL) errorPoly("No more memory.\n");
901: strcpy(ringName,rp->name);
902: strcat(ringName,pstr);
903: *newRingp = *rp;
904: newRingp->p = p;
905: newRingp->name = ringName;
906: return newRingp;
907: }
908:
909: /*
910: P = 3001;
911: L = [ ];
912: while (P<10000) {
913: L=cons(P,L);
914: P = pari(nextprime,P+1);
915: }
916: print(L);
917: */
918: #define N799 799
919: static int nextPrime(void) {
920: static int pt = 0;
921: static int tb[N799] =
922: {3001,3011,3019,3023,3037,3041,3049,3061,3067,3079,3083,3089,3109,3119,3121,3137,3163,3167,3169,3181,3187,3191,3203,3209,3217,3221,3229,3251,3253,3257,3259,3271,3299,3301,3307,3313,3319,3323,3329,3331,3343,3347,3359,3361,3371,3373,3389,3391,3407,3413,3433,3449,3457,3461,3463,3467,3469,3491,3499,3511,3517,3527,3529,3533,3539,3541,3547,3557,3559,3571,3581,3583,3593,3607,3613,3617,3623,3631,3637,3643,3659,3671,3673,3677,3691,3697,3701,3709,3719,3727,3733,3739,3761,3767,3769,3779,3793,3797,3803,3821,3823,3833,3847,3851,3853,3863,3877,3881,3889,3907,3911,3917,3919,3923,3929,3931,3943,3947,3967,3989,
923: 4001,4003,4007,4013,4019,4021,4027,4049,4051,4057,4073,4079,4091,4093,4099,4111,4127,4129,4133,4139,4153,4157,4159,4177,4201,4211,4217,4219,4229,4231,4241,4243,4253,4259,4261,4271,4273,4283,4289,4297,4327,4337,4339,4349,4357,4363,4373,4391,4397,4409,4421,4423,4441,4447,4451,4457,4463,4481,4483,4493,4507,4513,4517,4519,4523,4547,4549,4561,4567,4583,4591,4597,4603,4621,4637,4639,4643,4649,4651,4657,4663,4673,4679,4691,4703,4721,4723,4729,4733,4751,4759,4783,4787,4789,4793,4799,4801,4813,4817,4831,4861,4871,4877,4889,4903,4909,4919,4931,4933,4937,4943,4951,4957,4967,4969,4973,4987,4993,4999,
924: 5003,5009,5011,5021,5023,5039,5051,5059,5077,5081,5087,5099,5101,5107,5113,5119,5147,5153,5167,5171,5179,5189,5197,5209,5227,5231,5233,5237,5261,5273,5279,5281,5297,5303,5309,5323,5333,5347,5351,5381,5387,5393,5399,5407,5413,5417,5419,5431,5437,5441,5443,5449,5471,5477,5479,5483,5501,5503,5507,5519,5521,5527,5531,5557,5563,5569,5573,5581,5591,5623,5639,5641,5647,5651,5653,5657,5659,5669,5683,5689,5693,5701,5711,5717,5737,5741,5743,5749,5779,5783,5791,5801,5807,5813,5821,5827,5839,5843,5849,5851,5857,5861,5867,5869,5879,5881,5897,5903,5923,5927,5939,5953,5981,5987,
925: 6007,6011,6029,6037,6043,6047,6053,6067,6073,6079,6089,6091,6101,6113,6121,6131,6133,6143,6151,6163,6173,6197,6199,6203,6211,6217,6221,6229,6247,6257,6263,6269,6271,6277,6287,6299,6301,6311,6317,6323,6329,6337,6343,6353,6359,6361,6367,6373,6379,6389,6397,6421,6427,6449,6451,6469,6473,6481,6491,6521,6529,6547,6551,6553,6563,6569,6571,6577,6581,6599,6607,6619,6637,6653,6659,6661,6673,6679,6689,6691,6701,6703,6709,6719,6733,6737,6761,6763,6779,6781,6791,6793,6803,6823,6827,6829,6833,6841,6857,6863,6869,6871,6883,6899,6907,6911,6917,6947,6949,6959,6961,6967,6971,6977,6983,6991,6997,
926: 7001,7013,7019,7027,7039,7043,7057,7069,7079,7103,7109,7121,7127,7129,7151,7159,7177,7187,7193,7207,7211,7213,7219,7229,7237,7243,7247,7253,7283,7297,7307,7309,7321,7331,7333,7349,7351,7369,7393,7411,7417,7433,7451,7457,7459,7477,7481,7487,7489,7499,7507,7517,7523,7529,7537,7541,7547,7549,7559,7561,7573,7577,7583,7589,7591,7603,7607,7621,7639,7643,7649,7669,7673,7681,7687,7691,7699,7703,7717,7723,7727,7741,7753,7757,7759,7789,7793,7817,7823,7829,7841,7853,7867,7873,7877,7879,7883,7901,7907,7919,7927,7933,7937,7949,7951,7963,7993,
927: 8009,8011,8017,8039,8053,8059,8069,8081,8087,8089,8093,8101,8111,8117,8123,8147,8161,8167,8171,8179,8191,8209,8219,8221,8231,8233,8237,8243,8263,8269,8273,8287,8291,8293,8297,8311,8317,8329,8353,8363,8369,8377,8387,8389,8419,8423,8429,8431,8443,8447,8461,8467,8501,8513,8521,8527,8537,8539,8543,8563,8573,8581,8597,8599,8609,8623,8627,8629,8641,8647,8663,8669,8677,8681,8689,8693,8699,8707,8713,8719,8731,8737,8741,8747,8753,8761,8779,8783,8803,8807,8819,8821,8831,8837,8839,8849,8861,8863,8867,8887,8893,8923,8929,8933,8941,8951,8963,8969,8971,8999,
928: 9001,9007,9011,9013,9029,9041,9043,9049,9059,9067,9091,9103,9109,9127,9133,9137,9151,9157,9161,9173,9181,9187,9199,9203,9209,9221,9227,9239,9241,9257,9277,9281,9283,9293,9311,9319,9323,9337,9341,9343,9349,9371,9377,9391,9397,9403,9413,9419,9421,9431,9433,9437,9439,9461,9463,9467,9473,9479,9491,9497,9511,9521,9533,9539,9547,9551,9587,9601,9613,9619,9623,9629,9631,9643,9649,9661,9677,9679,9689,9697,9719,9721,9733,9739,9743,9749,9767,9769,9781,9787,9791,9803,9811,9817,9829,9833,9839,9851,9857,9859,9871,9883,9887,9901,9907,9923,9929,9931,9941,9949,9967,9973};
929:
930: if (pt <N799) {
931: return tb[pt++];
932: }else{
933: pt = 0;
934: return tb[pt++];
935: }
936: }
937:
938: int getPrime(int p) {
939: int i;
940: if (p <= 2) return nextPrime();
941: for (i=2; i<p; i++) {
942: if (p % i == 0) {
943: return nextPrime();
944: }
945: }
946: return p;
1.5 takayama 947: }
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