Annotation of OpenXM/src/kan96xx/Kan/poly4.c, Revision 1.14
1.14 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly4.c,v 1.13 2004/06/12 07:29:46 takayama Exp $ */
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[]);
1.10 takayama 11: static int degreeOfInitWS(POLY f,int w[],int s[]);
1.13 takayama 12: static int dDegree(POLY f);
13: static POLY dHomogenize(POLY f);
1.1 maekawa 14:
15: static void shell(v,n)
1.3 takayama 16: int v[];
17: int n;
1.1 maekawa 18: {
19: int gap,i,j,temp;
20:
21: for (gap = n/2; gap > 0; gap /= 2) {
22: for (i = gap; i<n; i++) {
23: for (j=i-gap ; j>=0 && v[j]<v[j+gap]; j -= gap) {
1.3 takayama 24: temp = v[j];
25: v[j] = v[j+gap];
26: v[j+gap] = temp;
1.1 maekawa 27: }
28: }
29: }
30: }
31:
32:
33: struct matrixOfPOLY *parts(f,v)
1.3 takayama 34: POLY f;
35: POLY v; /* v must be a single variable, e.g. x */
1.1 maekawa 36: {
37: struct matrixOfPOLY *evPoly;
38: int vi = 0; /* index of v */
39: int vx = 1; /* x --> 1, D--> 0 */
40: int n,evSize,i,k,e;
41: int *ev;
42: struct object *evList;
43: struct object *list;
44: struct object ob;
45: POLY ans;
46: POLY h;
47: extern struct ring *CurrentRingp;
48: POLY ft;
49:
50:
51: if (f ISZERO || v ISZERO) {
52: evPoly = newMatrixOfPOLY(2,1);
53: getMatrixOfPOLY(evPoly,0,0) = ZERO;
54: getMatrixOfPOLY(evPoly,1,0) = ZERO;
55: return(evPoly);
56: }
57: n = v->m->ringp->n;
58: /* get the index of the variable v */
59: for (i=0; i<n; i++) {
60: if (v->m->e[i].x) {
61: vx = 1; vi = i; break;
62: }else if (v->m->e[i].D) {
63: vx = 0; vi = i; break;
64: }
65: }
66: ft = f;
67: /* get the vector of exponents */
68: evList = NULLLIST;
69: while (ft != POLYNULL) {
70: if (vx) {
71: e = ft->m->e[vi].x;
72: }else{
73: e = ft->m->e[vi].D;
74: }
75: ft = ft->next;
76: ob = KpoInteger(e);
77: if (!memberQ(evList,ob)) {
78: list = newList(&ob);
79: evList = vJoin(evList,list);
80: }
81: }
82: /*printf("evList = "); printObjectList(evList);*/
83: evSize = klength(evList);
84: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
85: if (ev == (int *)NULL) errorPoly("No more memory.");
86: for (i=0; i<evSize; i++) {
87: ev[i] = KopInteger(car(evList));
88: evList = cdr(evList);
89: }
90: /* sort ev */
91: shell(ev,evSize);
92:
93: /* get coefficients */
94: evPoly = newMatrixOfPOLY(2,evSize);
95: for (i=0; i<evSize; i++) {
96: ans = ZERO;
97: getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp);
98: ft = f;
99: while (ft != POLYNULL) {
100: if (vx) {
1.3 takayama 101: if (ft->m->e[vi].x == ev[i]) {
102: h = newCell(ft->coeffp,monomialCopy(ft->m));
103: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
104: ans = ppAdd(ans,h);
105: }
1.1 maekawa 106: }else{
1.3 takayama 107: if (ft->m->e[vi].D == ev[i]) {
108: h = newCell(ft->coeffp,monomialCopy(ft->m));
109: dset0(h,vi);
110: ans = ppAdd(ans,h);
111: }
1.1 maekawa 112: }
113: ft = ft->next;
114: }
115: getMatrixOfPOLY(evPoly,1,i) = ans;
116: }
117: return(evPoly);
118: }
1.3 takayama 119:
1.1 maekawa 120: struct object parts2(f,v)
1.3 takayama 121: POLY f;
122: POLY v; /* v must be a single variable, e.g. x */
1.1 maekawa 123: {
124: struct matrixOfPOLY *evPoly;
125: int vi = 0; /* index of v */
126: int vx = 1; /* x --> 1, D--> 0 */
127: int n,evSize,i,k,e;
128: int *ev;
129: struct object *evList;
130: struct object *list;
131: struct object ob;
132: POLY ans;
133: POLY h;
134: POLY ft;
135: struct object ob1,ob2,rob;
136:
137:
138: if (f ISZERO || v ISZERO) {
139: evPoly = newMatrixOfPOLY(2,1);
140: getMatrixOfPOLY(evPoly,0,0) = ZERO;
141: getMatrixOfPOLY(evPoly,1,0) = ZERO;
142: rob = newObjectArray(2);
143: ob1 = newObjectArray(1);
144: ob2 = newObjectArray(1);
145: putoa(ob1,0,KpoInteger(0));
146: putoa(ob2,0,KpoPOLY(POLYNULL));
147: putoa(rob,0,ob1); putoa(rob,1,ob2);
148: return(rob);
149: }
150: n = v->m->ringp->n;
151: /* get the index of the variable v */
152: for (i=0; i<n; i++) {
153: if (v->m->e[i].x) {
154: vx = 1; vi = i; break;
155: }else if (v->m->e[i].D) {
156: vx = 0; vi = i; break;
157: }
158: }
159: ft = f;
160: /* get the vector of exponents */
161: evList = NULLLIST;
162: while (ft != POLYNULL) {
163: if (vx) {
164: e = ft->m->e[vi].x;
165: }else{
166: e = ft->m->e[vi].D;
167: }
168: ft = ft->next;
169: ob = KpoInteger(e);
170: if (!memberQ(evList,ob)) {
171: list = newList(&ob);
172: evList = vJoin(evList,list);
173: }
174: }
175: /*printf("evList = "); printObjectList(evList);*/
176: evSize = klength(evList);
177: ev = (int *)sGC_malloc(sizeof(int)*(evSize+1));
178: if (ev == (int *)NULL) errorPoly("No more memory.");
179: for (i=0; i<evSize; i++) {
180: ev[i] = KopInteger(car(evList));
181: evList = cdr(evList);
182: }
183: /* sort ev */
184: shell(ev,evSize);
185:
186: /* get coefficients */
187: evPoly = newMatrixOfPOLY(2,evSize);
188: for (i=0; i<evSize; i++) {
189: ans = ZERO;
190: /* getMatrixOfPOLY(evPoly,0,i) = cxx(ev[i],0,0,CurrentRingp); */
191: getMatrixOfPOLY(evPoly,0,i) = POLYNULL;
192: ft = f;
193: while (ft != POLYNULL) {
194: if (vx) {
1.3 takayama 195: if (ft->m->e[vi].x == ev[i]) {
196: h = newCell(ft->coeffp,monomialCopy(ft->m));
197: xset0(h,vi); /* touch monomial part, so you need to copy it above. */
198: ans = ppAdd(ans,h);
199: }
1.1 maekawa 200: }else{
1.3 takayama 201: if (ft->m->e[vi].D == ev[i]) {
202: h = newCell(ft->coeffp,monomialCopy(ft->m));
203: dset0(h,vi);
204: ans = ppAdd(ans,h);
205: }
1.1 maekawa 206: }
207: ft = ft->next;
208: }
209: getMatrixOfPOLY(evPoly,1,i) = ans;
210: }
211: rob = newObjectArray(2);
212: ob1 = newObjectArray(evSize);
213: ob2 = newObjectArray(evSize);
214: for (i=0; i<evSize; i++) {
215: putoa(ob2,i,KpoPOLY(getMatrixOfPOLY(evPoly,1,i)));
216: putoa(ob1,i,KpoInteger(ev[i]));
217: }
218: putoa(rob,0,ob1); putoa(rob,1,ob2);
219: return(rob);
220: }
1.3 takayama 221:
1.1 maekawa 222: int pDegreeWrtV(f,v)
1.3 takayama 223: POLY f;
224: POLY v;
1.1 maekawa 225: {
226: int vx = 1;
227: int vi = 0;
228: int i,n;
229: int ans;
230: if (f ISZERO || v ISZERO) return(0);
231: n = f->m->ringp->n;
232: for (i=0; i<n; i++) {
233: if (v->m->e[i].x) {
234: vx = 1; vi = i;
235: break;
236: }else if (v->m->e[i].D) {
237: vx = 0; vi = i;
238: break;
239: }
240: }
241: if (vx) {
242: ans = f->m->e[vi].x;
243: }else{
244: ans = f->m->e[vi].D;
245: }
246: f = f->next;
247: while (f != POLYNULL) {
248: if (vx) {
249: if (f->m->e[vi].x > ans) ans = f->m->e[vi].x;
250: }else{
251: if (f->m->e[vi].D > ans) ans = f->m->e[vi].D;
252: }
253: f = f->next;
254: }
255: return(ans);
256: }
257:
258: int containVectorVariable(POLY f)
259: {
260: MONOMIAL tf;
261: static int nn,mm,ll,cc,n,m,l,c;
262: static struct ring *cr = (struct ring *)NULL;
263: int i;
264:
265: if (f ISZERO) return(0);
266: tf = f->m;
267: if (tf->ringp != cr) {
268: n = tf->ringp->n;
269: m = tf->ringp->m;
270: l = tf->ringp->l;
271: c = tf->ringp->c;
272: nn = tf->ringp->nn;
273: mm = tf->ringp->mm;
274: ll = tf->ringp->ll;
275: cc = tf->ringp->cc;
276: cr = tf->ringp;
277: }
278:
279: while (f != POLYNULL) {
280: tf = f->m;
281: for (i=cc; i<c; i++) {
282: if ( tf->e[i].x ) return(1);
283: if ( tf->e[i].D ) return(1);
284: }
285: for (i=ll; i<l; i++) {
286: if (tf->e[i].x) return(1);
287: if (tf->e[i].D) return(1);
288: }
289: for (i=mm; i<m; i++) {
290: if (tf->e[i].x) return(1);
291: if (tf->e[i].D) return(1);
292: }
293: for (i=nn; i<n; i++) {
294: if (tf->e[i].x) return(1);
295: if (tf->e[i].D) return(1);
296: }
297: f = f->next;
298: }
299: return(0);
300:
301: }
302:
303: POLY homogenize(f)
1.3 takayama 304: POLY f;
305: /* homogenize by using (*grade)(f) */
1.1 maekawa 306: {
307: POLY t;
308: int maxg;
309: int flag,d;
1.13 takayama 310: extern int Homogenize;
1.1 maekawa 311:
312: if (f == ZERO) return(f);
1.13 takayama 313: if (Homogenize == 3) { /* double homogenization Dx x = x Dx + h H */
314: return dHomogenize(f);
315: }
1.1 maekawa 316: t = f; maxg = (*grade)(f); flag = 0;
317: while (t != POLYNULL) {
318: d = (*grade)(t);
319: if (d != maxg) flag = 1;
320: if (d > maxg) {
321: maxg = d;
322: }
323: t = t->next;
324: }
325: if (flag == 0) return(f);
326:
327: f = pmCopy(f); /* You can rewrite the monomial parts */
328: t = f;
329: while (t != POLYNULL) {
330: d = (*grade)(t);
331: if (d != maxg) {
332: t->m->e[0].D += maxg-d; /* Multiply h^(maxg-d) */
333: }
334: t = t->next;
335: }
336: return(f);
337: }
338:
339: int isHomogenized(f)
1.3 takayama 340: POLY f;
1.1 maekawa 341: {
342: POLY t;
343: extern int Homogenize_vec;
344: int maxg;
345: if (!Homogenize_vec) return(isHomogenized_vec(f));
346: if (f == ZERO) return(1);
1.4 takayama 347: if (f->m->ringp->weightedHomogenization) {
348: return 1; /* BUG: do not chech in case of one-zero homogenization */
349: }
1.1 maekawa 350: maxg = (*grade)(f);
351: t = f;
352: while (t != POLYNULL) {
353: if (maxg != (*grade)(t)) return(0);
354: t = t->next;
355: }
356: return(1);
357: }
358:
359: int isHomogenized_vec(f)
1.3 takayama 360: POLY f;
1.1 maekawa 361: {
1.3 takayama 362: /* This is not efficient version. *grade should be grade_module1v(). */
1.1 maekawa 363: POLY t;
364: int ggg;
365: if (f == ZERO) return(1);
1.4 takayama 366: if (f->m->ringp->weightedHomogenization) {
367: return 1; /* BUG: do not chech in case of one-zero homogenization */
368: }
1.1 maekawa 369: while (f != POLYNULL) {
370: t = f;
371: ggg = (*grade)(f);
372: while (t != POLYNULL) {
373: if ((*isSameComponent)(f,t)) {
1.3 takayama 374: if (ggg != (*grade)(t)) return(0);
1.1 maekawa 375: }
376: t = t->next;
377: }
378: f = f->next;
379: }
380: return(1);
381: }
382:
1.13 takayama 383: static POLY dHomogenize(f)
384: POLY f;
385: {
386: POLY t;
387: int maxg, maxdg;
388: int flag,d,dd,neg;
389:
390: if (f == ZERO) return(f);
1.14 ! takayama 391:
! 392: t = f;
! 393: maxg = (*grade)(f);
! 394: while (t != POLYNULL) {
! 395: dd = (*grade)(t);
! 396: if (maxg < dd) maxg = dd;
! 397: t = t->next;
! 398: }
! 399: /* fprintf(stderr,"maxg=%d\n",maxg); */
! 400:
! 401: t = f;
! 402: maxdg = dDegree(f);
! 403: while (t != POLYNULL) {
! 404: dd = dDegree(t);
! 405: if (maxdg < dd) maxdg = dd;
! 406: t = t->next;
! 407: }
! 408: /* fprintf(stderr,"maxdg=%d\n",maxdg); */
! 409:
! 410: t = f;
! 411: flag = 0;
1.13 takayama 412: while (t != POLYNULL) {
413: d = (*grade)(t);
414: if (d != maxg) flag = 1;
415: if (d > maxg) {
416: maxg = d;
417: }
418: d = dDegree(f);
419: if (d > maxdg) {
420: maxdg = d;
421: }
422: t = t->next;
423: }
424: if (flag == 0) return(f);
425:
426: t = f; neg = 0;
427: while (t != POLYNULL) {
428: d = (*grade)(t);
429: dd = dDegree(t);
430: if (maxg-d-(maxdg-dd) < neg) {
431: neg = maxg-d-(maxdg-dd);
432: }
433: t = t->next;
434: }
435: neg = -neg;
436:
437: f = pmCopy(f); /* You can rewrite the monomial parts */
438: t = f;
439: while (t != POLYNULL) {
440: d = (*grade)(t);
441: dd = dDegree(t);
442: t->m->e[0].D += maxdg-dd; /* h */
443: t->m->e[0].x += maxg-d-(maxdg-dd)+neg; /* Multiply H */
444: /* example Dx^2+Dx+x */
445: t = t->next;
446: }
447: return(f);
448: }
1.1 maekawa 449:
450: static int degreeOfPrincipalPart(f)
1.3 takayama 451: POLY f;
1.1 maekawa 452: {
453: int n,i,dd;
454: if (f ISZERO) return(0);
455: n = f->m->ringp->n; dd = 0;
456: /* D[0] is homogenization var */
457: for (i=1; i<n; i++) {
1.13 takayama 458: dd += f->m->e[i].D;
459: }
460: return(dd);
461: }
462:
463: static int dDegree(f)
464: POLY f;
465: {
466: int nn,i,dd,m;
467: if (f ISZERO) return(0);
468: nn = f->m->ringp->nn; dd = 0;
469: m = f->m->ringp->m;
470: for (i=m; i<nn; i++) {
1.1 maekawa 471: dd += f->m->e[i].D;
472: }
473: return(dd);
474: }
475:
476: POLY POLYToPrincipalPart(f)
1.3 takayama 477: POLY f;
1.1 maekawa 478: {
479: POLY node;
480: struct listPoly nod;
481: POLY h;
482: POLY g;
483: int maxd = -20000; /* very big negative number */
484: int dd;
485: node = &nod; node->next = POLYNULL; h = node;
486:
487: g = pCopy(f); /* shallow copy */
488: while (!(f ISZERO)) {
489: dd = degreeOfPrincipalPart(f);
490: if (dd > maxd) maxd = dd;
491: f = f->next;
492: }
493: while (!(g ISZERO)) {
494: dd = degreeOfPrincipalPart(g);
495: if (dd == maxd) {
496: h->next = g;
497: h = h->next;
498: }
499: g = g->next;
500: }
501: h->next = POLYNULL;
502: return(node->next);
503: }
504:
505: static int degreeOfInitW(f,w)
1.3 takayama 506: POLY f;
507: int w[];
1.1 maekawa 508: {
509: int n,i,dd;
510: if (f ISZERO) {
511: errorPoly("degreeOfInitW(0,w) ");
512: }
513: n = f->m->ringp->n; dd = 0;
514: for (i=0; i<n; i++) {
515: dd += (f->m->e[i].D)*w[n+i];
516: dd += (f->m->e[i].x)*w[i];
517: }
518: return(dd);
519: }
520:
521: POLY POLYToInitW(f,w)
1.3 takayama 522: POLY f;
523: int w[]; /* weight vector */
1.1 maekawa 524: {
525: POLY h;
526: POLY g;
527: int maxd;
528: int dd;
1.11 takayama 529: h = POLYNULL;
1.1 maekawa 530:
1.11 takayama 531: /*printf("1:%s\n",POLYToString(f,'*',1));*/
1.1 maekawa 532: if (f ISZERO) return(f);
533: maxd = degreeOfInitW(f,w);
1.11 takayama 534: g = f;
1.1 maekawa 535: while (!(f ISZERO)) {
536: dd = degreeOfInitW(f,w);
537: if (dd > maxd) maxd = dd;
538: f = f->next;
539: }
540: while (!(g ISZERO)) {
541: dd = degreeOfInitW(g,w);
542: if (dd == maxd) {
1.11 takayama 543: h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.1 maekawa 544: }
545: g = g->next;
546: }
1.11 takayama 547: /*printf("2:%s\n",POLYToString(h,'*',1));*/
548: return(h);
1.10 takayama 549: }
550:
551: static int degreeOfInitWS(f,w,s)
552: POLY f;
553: int w[];
554: int s[];
555: {
556: int n,i,dd;
557: if (f ISZERO) {
558: errorPoly("degreeOfInitWS(0,w) ");
559: }
1.11 takayama 560: if (s == (int *) NULL) return degreeOfInitW(f,w);
1.10 takayama 561: n = f->m->ringp->n; dd = 0;
562: for (i=0; i<n-1; i++) {
563: dd += (f->m->e[i].D)*w[n+i];
564: dd += (f->m->e[i].x)*w[i];
565: }
566: dd += s[(f->m->e[n-1].x)];
567: return(dd);
568: }
569:
570: POLY POLYToInitWS(f,w,s)
571: POLY f;
572: int w[]; /* weight vector */
573: int s[]; /* shift vector */
574: {
575: POLY h;
576: POLY g;
577: int maxd;
578: int dd;
1.11 takayama 579: h = POLYNULL;
1.10 takayama 580:
1.11 takayama 581: /*printf("1s:%s\n",POLYToString(f,'*',1));*/
1.10 takayama 582: if (f ISZERO) return(f);
583: maxd = degreeOfInitWS(f,w,s);
1.11 takayama 584: g = f;
1.10 takayama 585: while (!(f ISZERO)) {
586: dd = degreeOfInitWS(f,w,s);
587: if (dd > maxd) maxd = dd;
588: f = f->next;
589: }
590: while (!(g ISZERO)) {
591: dd = degreeOfInitWS(g,w,s);
592: if (dd == maxd) {
1.11 takayama 593: h = ppAdd(h,newCell(g->coeffp,g->m)); /* it might be slow. */
1.10 takayama 594: }
595: g = g->next;
596: }
1.11 takayama 597: /*printf("2s:%s\n",POLYToString(h,'*',1));*/
598: return(h);
1.10 takayama 599: }
600:
601: int ordWsAll(f,w,s)
602: POLY f;
603: int w[]; /* weight vector */
604: int s[]; /* shift vector */
605: {
606: int maxd;
607: int dd;
608:
609: if (f ISZERO) errorPoly("ordWsAll(0,w,s) ");
610: maxd = degreeOfInitWS(f,w,s);
611: while (!(f ISZERO)) {
612: dd = degreeOfInitWS(f,w,s);
613: if (dd > maxd) maxd = dd;
614: f = f->next;
615: }
616: return maxd;
1.1 maekawa 617: }
618:
619:
620: /*
621: 1.The substitution "ringp->multiplication = ...." is allowed only in
622: KsetUpRing(), so the check in KswitchFunction is not necessary.
623: 2.mmLarger != matrix and AvoidTheSameRing==1, then error
624: 3.If Schreyer = 1, then the system always generates a new ring.
625: 4.The execution of set_order_by_matrix is not allowed when Avoid... == 1.
626: 5.When mmLarger == tower (in tower.sm1, tower-sugar.sm1), we do
627: as follows with our own risk.
628: [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv
629: */
630: int isTheSameRing(struct ring *rstack[],int rp, struct ring *newRingp)
631: {
632: struct ring *rrr;
633: int i,j,k;
634: int a=0;
635: for (k=0; k<rp; k++) {
636: rrr = rstack[k];
637: if (rrr->p != newRingp->p) { a=1; goto bbb ; }
638: if (rrr->n != newRingp->n) { a=2; goto bbb ; }
639: if (rrr->nn != newRingp->nn) { a=3; goto bbb ; }
640: if (rrr->m != newRingp->m) { a=4; goto bbb ; }
641: if (rrr->mm != newRingp->mm) { a=5; goto bbb ; }
642: if (rrr->l != newRingp->l) { a=6; goto bbb ; }
643: if (rrr->ll != newRingp->ll) { a=7; goto bbb ; }
644: if (rrr->c != newRingp->c) { a=8; goto bbb ; }
645: if (rrr->cc != newRingp->cc) { a=9; goto bbb ; }
646: for (i=0; i<rrr->n; i++) {
647: if (strcmp(rrr->x[i],newRingp->x[i])!=0) { a=10; goto bbb ; }
648: if (strcmp(rrr->D[i],newRingp->D[i])!=0) { a=11; goto bbb ; }
649: }
650: if (rrr->orderMatrixSize != newRingp->orderMatrixSize) { a=12; goto bbb ; }
651: for (i=0; i<rrr->orderMatrixSize; i++) {
652: for (j=0; j<2*(rrr->n); j++) {
1.3 takayama 653: if (rrr->order[i*2*(rrr->n)+j] != newRingp->order[i*2*(rrr->n)+j])
654: { a=13; goto bbb ; }
1.1 maekawa 655: }
656: }
657: if (rrr->next != newRingp->next) { a=14; goto bbb ; }
658: if (rrr->multiplication != newRingp->multiplication) { a=15; goto bbb ; }
659: /* if (rrr->schreyer != newRingp->schreyer) { a=16; goto bbb ; }*/
660: if (newRingp->schreyer == 1) { a=16; goto bbb; }
1.4 takayama 661: if (rrr->weightedHomogenization != newRingp->weightedHomogenization) { a=16; goto bbb; }
1.7 takayama 662: if (rrr->degreeShiftSize != newRingp->degreeShiftSize) {
663: a = 17; goto bbb;
664: }
665: if (rrr->degreeShiftN != newRingp->degreeShiftN) {
666: a = 17; goto bbb;
667: }
668: for (i=0; i < rrr->degreeShiftSize; i++) {
669: for (j=0; j< rrr->degreeShiftN; j++) {
670: if (rrr->degreeShift[i*(rrr->degreeShiftN)+j] !=
671: newRingp->degreeShift[i*(rrr->degreeShiftN)+j]) {
672: a = 17; goto bbb;
673: }
674: }
675: }
676:
1.1 maekawa 677: /* The following fields are ignored.
678: void *gbListTower;
679: int *outputOrder;
680: char *name;
681: */
682: /* All tests are passed. */
683: return(k);
684: bbb: ;
685: /* for debugging. */
686: /* fprintf(stderr," reason=%d, ",a); */
687: }
688: return(-1);
689: }
1.6 takayama 690:
691: /* s->1 */
692: POLY goDeHomogenizeS(POLY f) {
1.9 takayama 693: POLY lRule[1];
694: POLY rRule[1];
695: struct ring *rp;
696: POLY ans;
697: /* printf("1:[%s]\n",POLYToString(f,'*',1)); */
698: if (f == POLYNULL) return f;
699: rp = f->m->ringp;
700: if (rp->next == NULL) {
701: lRule[0] = cxx(1,0,1,rp);
702: rRule[0] = cxx(1,0,0,rp);
703: ans=replace(f,lRule,rRule,1);
704: }else{
705: struct coeff *cp;
706: POLY t;
707: POLY nc;
708: ans = POLYNULL;
709: while (f != POLYNULL) {
710: cp = f->coeffp;
711: if (cp->tag == POLY_COEFF) {
712: t = goDeHomogenizeS((cp->val).f);
1.11 takayama 713: nc = newCell(polyToCoeff(t,f->m->ringp),monomialCopy(f->m));
1.9 takayama 714: ans = ppAddv(ans,nc);
715: f = f->next;
716: }else{
717: ans = f; break;
718: }
719: }
720: }
721: /* printf("2:[%s]\n",POLYToString(ans,'*',1)); */
722: return ans;
723: }
724:
725: POLY goDeHomogenizeS_buggy(POLY f) {
1.6 takayama 726: POLY node;
727: POLY lastf;
728: struct listPoly nod;
729: POLY h;
730: POLY tf;
731: int gt,first;
732:
1.9 takayama 733: printf("1:[%s]\n",POLYToString(f,'*',1));
1.6 takayama 734: if (f == POLYNULL) return(POLYNULL);
735: node = &nod;
736: node->next = POLYNULL;
737: lastf = POLYNULL;
738: first = 1;
739: while (f != POLYNULL) {
740: tf = newCell(f->coeffp,monomialCopy(f->m));
741: tf->m->e[0].x = 0; /* H, s variable in the G-O paper. */
742: if (first) {
743: node->next = tf;
744: lastf = tf;
745: first = 0;
746: }else{
747: gt = (*mmLarger)(lastf,tf);
748: if (gt == 1) {
749: lastf->next = tf;
750: lastf = tf;
751: }else{
752: h = node->next;
753: h = ppAddv(h,tf);
754: node->next = h;
755: lastf = h;
756: while (lastf->next != POLYNULL) {
757: lastf = lastf->next;
758: }
759: }
760: }
761: f = f->next;
762: }
1.9 takayama 763: printf("2:[%s]\n",POLYToString(node->next,'*',1));
1.6 takayama 764: return (node->next);
765: }
766:
1.5 takayama 767: /* Granger-Oaku's homogenization for the ecart tangent cone.
768: Note: 2003.07.10.
769: ds[] is the degree shift.
770: ei ( element index ). If it is < 0, then e[n-1]->x will be used,
771: else ei is used.
1.6 takayama 772: if onlyS is set to 1, then input is assumed to be (u,v)-h-homogeneous.
1.5 takayama 773: */
1.6 takayama 774: POLY goHomogenize(POLY f,int u[],int v[],int ds[],int dssize,int ei,int onlyS)
1.5 takayama 775: {
776: POLY node;
777: POLY lastf;
778: struct listPoly nod;
779: POLY h;
780: POLY tf;
781: int gt,first,m,mp,t,tp,dsIdx,message;
1.7 takayama 782: struct ring *rp;
1.5 takayama 783:
784: message = 1;
785: if (f == POLYNULL) return(POLYNULL);
1.7 takayama 786: rp = f->m->ringp;
1.12 takayama 787: /*
1.7 takayama 788: if ((rp->degreeShiftSize == 0) && (dssize > 0)) {
789: 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");
790: }
1.12 takayama 791: */
1.5 takayama 792: node = &nod;
793: node->next = POLYNULL;
794: lastf = POLYNULL;
795: first = 1;
796: while (f != POLYNULL) {
797: if (first) {
798: t = m = dGrade1(f);
799: tp = mp = uvGrade1(f,u,v,ds,dssize,ei);
800: }else{
801: t = dGrade1(f);
802: tp = uvGrade1(f,u,v,ds,dssize,ei);
803: if (t > m) m = t;
804: if (tp < mp) mp = tp;
805: }
806:
807: tf = newCell(f->coeffp,monomialCopy(f->m));
808: /* Automatic dehomogenize. Not += */
809: if (message && ((tf->m->e[0].D != 0) || (tf->m->e[0].x != 0))) {
810: /*go-debug fprintf(stderr,"Automatic dehomogenize and homogenize.\n"); */
811: message = 0;
812: }
1.6 takayama 813: if (!onlyS) {
814: tf->m->e[0].D = -t; /* h */
815: }
1.5 takayama 816: tf->m->e[0].x = tp; /* H, s variable in the G-O paper. */
817: /*go-debug printf("t(h)=%d, tp(uv+ds)=%d\n",t,tp); */
818: if (first) {
819: node->next = tf;
820: lastf = tf;
821: first = 0;
822: }else{
823: gt = (*mmLarger)(lastf,tf);
824: if (gt == 1) {
825: lastf->next = tf;
826: lastf = tf;
827: }else{
828: /*go-debug printf("?\n"); */
829: h = node->next;
830: h = ppAddv(h,tf);
831: node->next = h;
832: lastf = h;
833: while (lastf->next != POLYNULL) {
834: lastf = lastf->next;
835: }
836: }
837: }
838: f = f->next;
839: }
840: h = node->next;
841: /*go-debug printf("m=%d, mp=%d\n",m,mp); */
842: while (h != POLYNULL) {
1.12 takayama 843: /*go-debug printf("Old: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x); */
1.6 takayama 844: if (!onlyS) h->m->e[0].D += m; /* h */
1.5 takayama 845: h->m->e[0].x += -mp; /* H, s*/
846: /*go-debug printf("New: h=%d, s=%d\n",h->m->e[0].D,h->m->e[0].x); */
847: h = h->next;
848: }
849: return (node->next);
850: }
851:
852: /* u[] = -1, v[] = 1 */
1.6 takayama 853: POLY goHomogenize11(POLY f,int ds[],int dssize,int ei,int onlyS)
1.5 takayama 854: {
855: int r;
856: int i,t,n,m,nn;
857: MONOMIAL tf;
858: static int *u;
859: static int *v;
860: static struct ring *cr = (struct ring *)NULL;
861:
862: if (f == POLYNULL) return POLYNULL;
863:
864: tf = f->m;
865: if (tf->ringp != cr) {
866: n = tf->ringp->n;
867: m = tf->ringp->m;
868: nn = tf->ringp->nn;
869: cr = tf->ringp;
870: u = (int *)sGC_malloc(sizeof(int)*n);
871: v = (int *)sGC_malloc(sizeof(int)*n);
872: for (i=0; i<n; i++) u[i]=v[i]=0;
873: for (i=m; i<nn; i++) {
874: u[i] = -1; v[i] = 1;
875: }
876: }
1.6 takayama 877: return(goHomogenize(f,u,v,ds,dssize,ei,onlyS));
1.5 takayama 878: }
879:
1.6 takayama 880: POLY goHomogenize_dsIdx(POLY f,int u[],int v[],int dsIdx,int ei,int onlyS)
1.5 takayama 881: {
882: if (f == POLYNULL) return POLYNULL;
883: }
1.6 takayama 884: POLY goHomogenize11_dsIdx(POLY f,int ds[],int dsIdx,int ei,int onlyS)
1.5 takayama 885: {
886: if (f == POLYNULL) return POLYNULL;
1.8 takayama 887: }
888:
889: /* cf. KsetUpRing() in kanExport0.c */
890: struct ring *newRingOverFp(struct ring *rp,int p) {
891: struct ring *newRingp;
892: char *ringName = NULL;
893: char pstr[64];
894: sprintf(pstr,"_%d",p);
895: ringName = (char *)sGC_malloc(128);
896: newRingp = (struct ring *)sGC_malloc(sizeof(struct ring));
897: if (newRingp == NULL) errorPoly("No more memory.\n");
898: strcpy(ringName,rp->name);
899: strcat(ringName,pstr);
900: *newRingp = *rp;
901: newRingp->p = p;
902: newRingp->name = ringName;
903: return newRingp;
904: }
905:
906: /*
907: P = 3001;
908: L = [ ];
909: while (P<10000) {
910: L=cons(P,L);
911: P = pari(nextprime,P+1);
912: }
913: print(L);
914: */
915: #define N799 799
916: static int nextPrime(void) {
917: static int pt = 0;
918: static int tb[N799] =
919: {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,
920: 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,
921: 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,
922: 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,
923: 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,
924: 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,
925: 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};
926:
927: if (pt <N799) {
928: return tb[pt++];
929: }else{
930: pt = 0;
931: return tb[pt++];
932: }
933: }
934:
935: int getPrime(int p) {
936: int i;
937: if (p <= 2) return nextPrime();
938: for (i=2; i<p; i++) {
939: if (p % i == 0) {
940: return nextPrime();
941: }
942: }
943: return p;
1.5 takayama 944: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>