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