Annotation of OpenXM/src/kan96xx/Kan/order.c, Revision 1.14
1.14 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/order.c,v 1.13 2004/09/13 11:24:11 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:
8: /* The format of order.
9: Example: graded lexicographic order
10: x_{N-1} x_{N-2} ... x_0 D_{N-1} .... D_{0}
11: 1 1 1 1 1
12: 1 0 0 0 0
13: 0 1 0 0 0
14: ..............................................
15:
16: (ringp->order)[i][j] should be (ringp->order)[i*2*N+j].
17: All order matrix is generated by functions in smacro.sm1
18: */
19:
20: static void warningOrder(char *s);
21: static void errorOrder(char *s);
22:
23: void setOrderByMatrix(order,n,c,l,omsize)
1.4 takayama 24: int order[];
25: int n,c,l,omsize;
1.1 maekawa 26: {
27: int i,j;
28: int *Order;
29: extern struct ring *CurrentRingp;
30:
31: switch_mmLarger("default");
1.4 takayama 32: /* q-case */
1.1 maekawa 33: if ( l-c > 0) {
34: switch_mmLarger("qmatrix");
35: }
36:
37: Order = (int *)sGC_malloc(sizeof(int)*(2*n)*(omsize));
38: if (Order == (int *)NULL) errorOrder("No memory.");
39: CurrentRingp->order = Order;
40: CurrentRingp->orderMatrixSize = omsize;
41: for (i=0; i<omsize; i++) {
42: for (j=0; j<2*n; j++) {
43: Order[i*2*n+j] = order[i*2*n+j];
44: }
45: }
46: }
47:
48: void showRing(level,ringp)
1.4 takayama 49: int level;
50: struct ring *ringp;
1.1 maekawa 51: {
52: int i,j;
53: FILE *fp;
54: char tmp[100];
55: int N,M,L,C,NN,MM,LL,CC;
56: char **TransX,**TransD;
57: int *Order;
58: int P;
59: char *mtype;
60: extern char *F_isSameComponent;
1.5 takayama 61: POLY f;
1.6 takayama 62: POLY fx;
63: POLY fd;
64: POLY rf;
1.1 maekawa 65: fp = stdout;
66:
67: N=ringp->n; M = ringp->m; L = ringp->l; C = ringp->c;
68: NN=ringp->nn; MM = ringp->mm; LL = ringp->ll; CC = ringp->cc;
69: TransX = ringp->x; TransD = ringp->D;
70: Order = ringp->order;
71: P = ringp->p;
72:
73:
74: fprintf(fp,"\n---------- the current ring ---- name: %s------\n",ringp->name);
75: fprintf(fp,"Characteristic is %d. ",P);
76: fprintf(fp,"N0=%d N=%d NN=%d M=%d MM=%d L=%d LL=%d C=%d CC=%d omsize=%d\n",N0,N,NN,M,MM,L,LL,C,CC,ringp->orderMatrixSize);
77: fprintf(fp,"\n");
78:
79: /* print identifier names */
80: if (N-M >0) {
81: fprintf(fp,"Differential variables: ");
82: for (i=M; i<N; i++) fprintf(fp," %4s ",TransX[i]);
83: for (i=M; i<N; i++) fprintf(fp," %4s ",TransD[i]);
84: fprintf(fp,"\n");
85: fprintf(fp,"where ");
86: for (i=M; i<N; i++) {
1.6 takayama 87: fx = cxx(1,i,1,ringp); fd = cdd(1,i,1,ringp);
88: rf = ppSub(ppMult(fd,fx),ppMult(fx,fd));
89: fprintf(fp," %s %s - %s %s = %s, ",TransD[i],TransX[i],
90: TransX[i],TransD[i],POLYToString(rf,'*',0));
1.1 maekawa 91: }
92: fprintf(fp,"\n\n");
93: }
94: if (M-L >0) {
95: fprintf(fp,"Difference variables: ");
96: for (i=L; i<M; i++) fprintf(fp," %4s ",TransX[i]);
97: for (i=L; i<M; i++) fprintf(fp," %4s ",TransD[i]);
98: fprintf(fp,"\n");
99: fprintf(fp,"where ");
100: for (i=L; i<M; i++) {
1.5 takayama 101: fprintf(fp," %s %s - %s %s = ",TransD[i],TransX[i],
102: TransX[i],TransD[i]);
103: f=ppSub(ppMult(cdd(1,i,1,ringp),cxx(1,i,1,ringp)),
104: ppMult(cxx(1,i,1,ringp),cdd(1,i,1,ringp)));
105: fprintf(fp," %s, ",POLYToString(f,'*',0));
1.1 maekawa 106: }
107: fprintf(fp,"\n\n");
108: }
109: if (L-C >0) {
110: fprintf(fp,"q-Difference variables: ");
111: for (i=C; i<L; i++) fprintf(fp," %4s ",TransX[i]);
112: for (i=C; i<L; i++) fprintf(fp," %4s ",TransD[i]);
113: fprintf(fp,"\n");
114: fprintf(fp,"where ");
115: for (i=C; i<L; i++) {
116: fprintf(fp," %s %s = %s %s %s, ",TransD[i],TransX[i],
1.4 takayama 117: TransX[0],
118: TransX[i],TransD[i]);
1.1 maekawa 119: }
120: fprintf(fp,"\n\n");
121: }
122: if (C>0) {
123: fprintf(fp,"Commutative variables: ");
124: for (i=0; i<C; i++) fprintf(fp," %4s ",TransX[i]);
125: for (i=0; i<C; i++) fprintf(fp," %4s ",TransD[i]);
126: fprintf(fp,"\n\n");
127: }
128:
129: if (strcmp(F_isSameComponent,"x") == 0) {
130: fprintf(fp,"Integral or summation or graduation variables are : ");
131: for (i=CC; i<C; i++) fprintf(fp," %4s ",TransX[i]);
132: for (i=LL; i<L; i++) fprintf(fp," %4s ",TransX[i]);
133: for (i=MM; i<M; i++) fprintf(fp," %4s ",TransX[i]);
134: for (i=NN; i<N; i++) fprintf(fp," %4s ",TransX[i]);
135: fprintf(fp,"\n");
136: }else if (strcmp(F_isSameComponent,"xd") == 0) {
137: fprintf(fp,"Graduation variables are : ");
138: for (i=CC; i<C; i++) fprintf(fp," %4s ",TransX[i]);
139: for (i=LL; i<L; i++) fprintf(fp," %4s ",TransX[i]);
140: for (i=MM; i<M; i++) fprintf(fp," %4s ",TransX[i]);
141: for (i=NN; i<N; i++) fprintf(fp," %4s ",TransX[i]);
142: for (i=CC; i<C; i++) fprintf(fp," %4s ",TransD[i]);
143: for (i=LL; i<L; i++) fprintf(fp," %4s ",TransD[i]);
144: for (i=MM; i<M; i++) fprintf(fp," %4s ",TransD[i]);
145: for (i=NN; i<N; i++) fprintf(fp," %4s ",TransD[i]);
146: fprintf(fp,"\n");
147: }else {
148: fprintf(fp,"Unknown graduation variable specification.\n\n");
149: }
150: fprintf(fp,"The homogenization variable is : ");
151: fprintf(fp," %4s ",TransD[0]);
152: fprintf(fp,"\n");
153:
154:
155:
156: fprintf(fp,"-------------------------------------------\n");
157: fprintf(fp,"Output order : ");
158: for (i=0; i<2*N; i++) {
159: if (ringp->outputOrder[i] < N) {
160: fprintf(fp,"%s ",TransX[ringp->outputOrder[i]]);
161: }else{
162: fprintf(fp,"%s ",TransD[(ringp->outputOrder[i])-N]);
163: }
164: }
165: fprintf(fp,"\n");
166:
167: if (ringp->multiplication == mpMult_poly) {
168: mtype = "poly";
169: }else if (ringp->multiplication == mpMult_diff) {
170: mtype = "diff";
171: }else if (ringp->multiplication == mpMult_difference) {
172: mtype = "difference";
173: }else {
174: mtype = "unknown";
175: }
176: fprintf(fp,"Multiplication function --%s(%xH).\n",
1.4 takayama 177: mtype,(unsigned int) ringp->multiplication);
1.1 maekawa 178: if (ringp->schreyer) {
179: fprintf(fp,"schreyer=1, gbListTower=");
180: printObjectList((struct object *)(ringp->gbListTower));
181: fprintf(fp,"\n");
182: }
1.7 takayama 183: if (ringp->degreeShiftSize) {
1.8 takayama 184: fprintf(fp,"degreeShift vector (N=%d,Size=%d)= \n[\n",ringp->degreeShiftN,ringp->degreeShiftSize);
1.7 takayama 185: {
1.8 takayama 186: int i,j;
187: for (i=0; i<ringp->degreeShiftN; i++) {
188: fprintf(fp," [");
189: for (j=0; j< ringp->degreeShiftSize; j++) {
190: fprintf(fp," %d ",ringp->degreeShift[i*(ringp->degreeShiftSize)+j]);
191: }
192: fprintf(fp,"]\n");
1.7 takayama 193: }
194: }
195: fprintf(fp,"]\n");
196: }
197: fprintf(fp,"--- weight vectors ---\n");
1.1 maekawa 198: if (level) printOrder(ringp);
1.13 takayama 199:
200: if (ringp->partialEcart) {
201: fprintf(fp,"--- partialEcartGlobalVarX ---\n");
202: for (i=0; i<ringp->partialEcart; i++) {
203: fprintf(fp," %4s ",TransX[ringp->partialEcartGlobalVarX[i]]);
204: }
205: fprintf(fp,"\n");
206: }
1.1 maekawa 207:
208: if (ringp->next != (struct ring *)NULL) {
209: fprintf(fp,"\n\n-------- The next ring is .... --------------\n");
210: showRing(level,ringp->next);
211: }
212: }
213:
214: /***************************************************************
215: functions related to order
216: ******************************************************************/
217: #define xtoi(k) ((N-1)-(k))
218: #define dtoi(k) ((2*N-1)-(k))
219: #define itox(k) ((N-1)-(k))
220: #define itod(k) ((2*N-1)-(k))
221: #define isX(i) (i<N? 1: 0)
222: #define isD(i) (i<N? 0: 1)
223: /****************************************************
224: i : 0 1 N-1 N 2N-1
225: x :x_{N-1} x_{N-2} x_0
226: d : D_{N-1} D_{0}
227: if (isX(i)) x_{itox(i)}
228: if (isD(i)) D_{itod(i)}
229: ******************************************************/
230: /* xtoi(0):N-1 xtoi(1):N-2 ....
231: dtoi(0):2N-1 dtoi(1):2N-2 ...
232: itod(N):N-1 dtoi(N-1):N ...
233: */
234:
235: void printOrder(ringp)
1.4 takayama 236: struct ring *ringp;
1.1 maekawa 237: {
238: int i,j;
239: FILE *fp;
240: char tmp[100];
241: int N,M,L,C,NN,MM,LL,CC;
242: char **TransX,**TransD;
243: int *Order;
244: int P;
245: int omsize;
246: extern char *F_isSameComponent;
247:
248: N=ringp->n; M = ringp->m; L = ringp->l; C = ringp->c;
249: NN=ringp->nn; MM = ringp->mm; LL = ringp->ll; CC = ringp->cc;
250: TransX = ringp->x; TransD = ringp->D;
251: Order = ringp->order;
252: P = ringp->p;
253: omsize = ringp->orderMatrixSize;
254:
255: fp = stdout;
256:
257:
258: for (i=0; i<2*N; i++) printf("%4d",i);
259: fprintf(fp,"\n");
260:
261: /* print variables names */
262: for (i=0; i<N; i++) {
263: sprintf(tmp,"x%d",N-1-i);
264: fprintf(fp,"%4s",tmp);
265: }
266: for (i=0; i<N; i++) {
267: sprintf(tmp,"D%d",N-1-i);
268: fprintf(fp,"%4s",tmp);
269: }
270: fprintf(fp,"\n");
271:
272: /* print identifier names */
273: for (i=0; i<N; i++) fprintf(fp,"%4s",TransX[itox(i)]);
274: for (i=N; i<2*N; i++) fprintf(fp,"%4s",TransD[itod(i)]);
275: fprintf(fp,"\n");
276:
277: /* print D: differential DE: differential, should be eliminated
1.4 takayama 278: E: difference
279: Q: q-difference
280: C: commutative
1.1 maekawa 281: */
282: if (strcmp(F_isSameComponent,"x")== 0 || strcmp(F_isSameComponent,"xd")==0) {
283: for (i=0; i<N; i++) {
284: if ((NN<=itox(i)) && (itox(i)<N)) fprintf(fp,"%4s","DE");
285: if ((M<=itox(i)) && (itox(i)<NN)) fprintf(fp,"%4s","D");
286: if ((MM<=itox(i)) && (itox(i)<M)) fprintf(fp,"%4s","EE");
287: if ((L<=itox(i)) && (itox(i)<MM)) fprintf(fp,"%4s","E");
288: if ((LL<=itox(i)) && (itox(i)<L)) fprintf(fp,"%4s","QE");
289: if ((C<=itox(i)) && (itox(i)<LL)) fprintf(fp,"%4s","Q");
290: if ((CC<=itox(i)) && (itox(i)<C)) fprintf(fp,"%4s","CE");
291: if ((0<=itox(i)) && (itox(i)<CC)) fprintf(fp,"%4s","C");
292: }
293: }
294: if (strcmp(F_isSameComponent,"x")==0) {
295: for (i=N; i<2*N; i++) {
296: if ((M<=itod(i)) && (itod(i)<N)) fprintf(fp,"%4s","D");
297: if ((L<=itod(i)) && (itod(i)<M)) fprintf(fp,"%4s","E");
298: if ((C<=itod(i)) && (itod(i)<L)) fprintf(fp,"%4s","Q");
299: if ((0<=itod(i)) && (itod(i)<C)) fprintf(fp,"%4s","C");
300: }
301: }else if (strcmp(F_isSameComponent,"xd")==0) {
302: for (i=N; i<2*N; i++) {
303: if ((NN<=itod(i)) && (itod(i)<N)) fprintf(fp,"%4s","DE");
304: if ((M<=itod(i)) && (itod(i)<NN)) fprintf(fp,"%4s","D");
305: if ((MM<=itod(i)) && (itod(i)<M)) fprintf(fp,"%4s","EE");
306: if ((L<=itod(i)) && (itod(i)<MM)) fprintf(fp,"%4s","E");
307: if ((LL<=itod(i)) && (itod(i)<L)) fprintf(fp,"%4s","QE");
308: if ((C<=itod(i)) && (itod(i)<LL)) fprintf(fp,"%4s","Q");
309: if ((CC<=itod(i)) && (itod(i)<C)) fprintf(fp,"%4s","CE");
310: if ((0<=itod(i)) && (itod(i)<CC)) fprintf(fp,"%4s","C");
311: }
312: } else {
313: fprintf(fp,"Unknown graduation variable type.\n");
314: }
315: fprintf(fp,"\n");
316:
317: for (i=0; i< omsize; i++) {
318: for (j=0; j<2*N; j++) {
319: fprintf(fp,"%4d", Order[i*2*N+j]);
320: }
321: fprintf(fp,"\n");
322: }
323: fprintf(fp,"\n");
324:
325: }
326:
327: struct object oGetOrderMatrix(struct ring *ringp)
328: {
1.14 ! takayama 329: struct object rob = OINIT;
! 330: struct object ob2 = OINIT;
1.1 maekawa 331: int n,i,j,m;
332: int *om;
333: n = ringp->n;
334: m = ringp->orderMatrixSize;
335: om = ringp->order;
336: if (m<=0) m = 1;
337: rob = newObjectArray(m);
338: for (i=0; i<m; i++) {
339: ob2 = newObjectArray(2*n);
340: for (j=0; j<2*n; j++) {
341: putoa(ob2,j,KpoInteger(om[2*n*i+j]));
342: }
343: putoa(rob,i,ob2);
344: }
345: return(rob);
346: }
347:
348:
349: int mmLarger_matrix(ff,gg)
1.4 takayama 350: POLY ff; POLY gg;
1.1 maekawa 351: {
352: int exp[2*N0]; /* exponents */
353: int i,k;
354: int sum,flag;
355: int *Order;
356: int N;
357: MONOMIAL f,g;
358: struct ring *rp;
359: int in2;
360: int *from, *to;
361: int omsize;
1.7 takayama 362: int dssize;
1.8 takayama 363: int dsn;
1.7 takayama 364: int *degreeShiftVector;
1.1 maekawa 365:
366: if (ff == POLYNULL ) {
367: if (gg == POLYNULL) return( 2 );
368: else return( 0 );
369: }
370: if (gg == POLYNULL) {
371: if (ff == POLYNULL) return( 2 );
372: else return( 1 );
373: }
374: f = ff->m; g=gg->m;
375:
376: rp = f->ringp;
377: Order = rp->order;
378: N = rp->n;
379: from = rp->from;
380: to = rp->to;
381: omsize = rp->orderMatrixSize;
1.7 takayama 382: if (dssize = rp->degreeShiftSize) {
383: degreeShiftVector = rp->degreeShift; /* Note. 2003.06.26 */
1.8 takayama 384: dsn = rp->degreeShiftN;
1.7 takayama 385: }
1.1 maekawa 386:
387: flag = 1;
388: for (i=N-1,k=0; i>=0; i--,k++) {
389: exp[k] = (f->e[i].x) - (g->e[i].x);
390: exp[k+N] = (f->e[i].D) - (g->e[i].D);
391: if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
392: }
393: if (flag==1) return(2);
394: /* exp > 0 <---> f>g
395: exp = 0 <---> f=g
396: exp < 0 <---> f<g
397: */
398: for (i=0; i< omsize; i++) {
399: sum = 0; in2 = i*2*N;
400: /* for (k=0; k<2*N; k++) sum += exp[k]*Order[in2+k]; */
401: for (k=from[i]; k<to[i]; k++) sum += exp[k]*Order[in2+k];
1.8 takayama 402: if (dssize && ( i < dsn)) { /* Note, 2003.06.26 */
1.7 takayama 403: if ((f->e[N-1].x < dssize) && (f->e[N-1].x >= 0) &&
404: (g->e[N-1].x < dssize) && (g->e[N-1].x >= 0)) {
1.8 takayama 405: sum += degreeShiftVector[i*dssize+ (f->e[N-1].x)]
406: -degreeShiftVector[i*dssize+ (g->e[N-1].x)];
1.7 takayama 407: }else{
1.9 takayama 408: /*warningOrder("Size mismatch in the degree shift vector. It is ignored.");*/
1.7 takayama 409: }
410: }
1.1 maekawa 411: if (sum > 0) return(1);
412: if (sum < 0) return(0);
413: }
414: return(2);
415: }
416:
417: /* This should be used in case of q */
418: int mmLarger_qmatrix(ff,gg)
1.4 takayama 419: POLY ff; POLY gg;
1.1 maekawa 420: {
421: int exp[2*N0]; /* exponents */
422: int i,k;
423: int sum,flag;
424: int *Order;
425: int N;
426: MONOMIAL f,g;
427: int omsize;
428:
429: if (ff == POLYNULL ) {
430: if (gg == POLYNULL) return( 2 );
431: else return( 0 );
432: }
433: if (gg == POLYNULL) {
434: if (ff == POLYNULL) return( 2 );
435: else return( 1 );
436: }
437: f = ff->m; g = gg->m;
438: Order = f->ringp->order;
439: N = f->ringp->n;
440: omsize = f->ringp->orderMatrixSize;
441:
442: flag = 1;
443: for (i=N-1,k=0; i>=0; i--,k++) {
444: exp[k] = (f->e[i].x) - (g->e[i].x);
445: exp[k+N] = (f->e[i].D) - (g->e[i].D);
446: if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
447: }
448: if (flag==1) return(2);
449: /* exp > 0 <---> f>g
450: exp = 0 <---> f=g
451: exp < 0 <---> f<g
452: */
453: for (i=0; i< omsize; i++) {
454: sum = 0;
455: /* In case of q, you should do as follows */
456: for (k=0; k<N-1; k++) sum += exp[k]*Order[i*2*N+k]; /* skip k= N-1 -->q */
457: for (k=N; k<2*N-1; k++) sum += exp[k]*Order[i*2*N+k]; /* SKip k= 2*N-1 */
458: if (sum > 0) return(1);
459: else if (sum < 0) return(0);
460: }
461: if (exp[N-1] > 0) return(1);
462: else if (exp[N-1] < 0) return(0);
463: else return(2);
464: }
465:
466: /* x(N-1)>x(N-2)>....>D(N-1)>....>D(0) */
467: mmLarger_pureLexicographic(f,g)
1.4 takayama 468: POLY f;
469: POLY g;
1.1 maekawa 470: {
471: int i,r;
472: int n;
473: MONOMIAL fm,gm;
474: /* Note that this function ignores the order matrix of the given
475: ring. */
476: if (f == POLYNULL ) {
477: if (g == POLYNULL) return( 2 );
478: else return( 0 );
479: }
480: if (g == POLYNULL) {
481: if (f == POLYNULL) return( 2 );
482: else return( 1 );
483: }
484:
485:
486: fm = f->m; gm = g->m;
487: n = fm->ringp->n;
488: for (i=n-1; i>=0; i--) {
489: r = (fm->e[i].x) - (gm->e[i].x);
490: if (r > 0) return(1);
491: else if (r < 0) return(0);
492: else ;
493: }
494:
495: for (i=n-1; i>=0; i--) {
496: r = (fm->e[i].D) - (gm->e[i].D);
497: if (r > 0) return(1);
498: else if (r < 0) return(0);
499: else ;
500: }
501:
502: return(2);
503:
504: }
505:
506:
507: void setFromTo(ringp)
1.4 takayama 508: struct ring *ringp;
1.1 maekawa 509: {
510: int n;
511: int i,j,oasize;
512: if (ringp->order == (int *)NULL) errorOrder("setFromTo(); no order matrix.");
513: n = (ringp->n)*2;
514: oasize = ringp->orderMatrixSize;
515: ringp->from = (int *)sGC_malloc(sizeof(int)*oasize);
516: ringp->to = (int *)sGC_malloc(sizeof(int)*oasize);
517: if (ringp->from == (int *)NULL || ringp->to == (int *)NULL) {
518: errorOrder("setFromTo(): No memory.");
519: }
520: for (i=0; i<oasize; i++) {
521: ringp->from[i] = 0; ringp->to[i] = n;
522: for (j=0; j<n; j++) {
523: if (ringp->order[i*n+j] != 0) {
1.4 takayama 524: ringp->from[i] = j;
525: break;
1.1 maekawa 526: }
527: }
528: for (j=n-1; j>=0; j--) {
529: if (ringp->order[i*n+j] != 0) {
1.4 takayama 530: ringp->to[i] = j+1;
531: break;
1.1 maekawa 532: }
533: }
534: }
535: }
536:
537: /* It ignores h and should be used with mmLarger_tower */
538: /* cf. mmLarger_matrix. h always must be checked at last. */
539: static int mmLarger_matrix_schreyer(ff,gg)
1.4 takayama 540: POLY ff; POLY gg;
1.1 maekawa 541: {
542: int exp[2*N0]; /* exponents */
543: int i,k;
544: int sum,flag;
545: int *Order;
546: int N;
547: MONOMIAL f,g;
548: struct ring *rp;
549: int in2;
550: int *from, *to;
551: int omsize;
552:
553: if (ff == POLYNULL ) {
554: if (gg == POLYNULL) return( 2 );
555: else return( 0 );
556: }
557: if (gg == POLYNULL) {
558: if (ff == POLYNULL) return( 2 );
559: else return( 1 );
560: }
561: f = ff->m; g=gg->m;
562:
563: rp = f->ringp;
564: Order = rp->order;
565: N = rp->n;
566: from = rp->from;
567: to = rp->to;
568: omsize = rp->orderMatrixSize;
569:
570: flag = 1;
571: for (i=N-1,k=0; i>0; i--,k++) {
572: exp[k] = (f->e[i].x) - (g->e[i].x);
573: exp[k+N] = (f->e[i].D) - (g->e[i].D);
574: if ((exp[k] != 0) || (exp[k+N] != 0)) flag =0;
575: }
576: exp[N-1] = (f->e[0].x) - (g->e[0].x);
577: exp[2*N-1] = 0; /* f->e[0].D - g->e[0].D. Ignore h! */
578: if ((exp[N-1] != 0) || (exp[2*N-1] != 0)) flag =0;
579:
580: if (flag==1) return(2);
581: /* exp > 0 <---> f>g
582: exp = 0 <---> f=g
583: exp < 0 <---> f<g
584: */
585: for (i=0; i< omsize; i++) {
586: sum = 0; in2 = i*2*N;
587: /* for (k=0; k<2*N; k++) sum += exp[k]*Order[in2+k]; */
588: for (k=from[i]; k<to[i]; k++) sum += exp[k]*Order[in2+k];
589: if (sum > 0) return(1);
590: if (sum < 0) return(0);
591: }
592: return(2);
593: }
594:
595: int mmLarger_tower(POLY f,POLY g) {
596: struct object *gbList;
597: int r;
598: if (f == POLYNULL) {
599: if (g == POLYNULL) return(2);
600: else return(0);
601: }
602: if (g == POLYNULL) {
603: if (f == POLYNULL) return(2);
604: else return(1);
605: }
606: if (!(f->m->ringp->schreyer) || !(g->m->ringp->schreyer))
607: return(mmLarger_matrix(f,g));
1.4 takayama 608: /* modifiable: mmLarger_qmatrix */
1.1 maekawa 609: gbList = (struct object *)(g->m->ringp->gbListTower);
610: if (gbList == NULL) return(mmLarger_matrix(f,g));
1.4 takayama 611: /* modifiable: mmLarger_qmatrix */
1.1 maekawa 612: if (gbList->tag != Slist) {
613: warningOrder("mmLarger_tower(): gbList must be in Slist.\n");
614: return(1);
615: }
616: if (klength(gbList) ==0) return(mmLarger_matrix(f,g));
1.4 takayama 617: /* modifiable: mmLarger_qmatrix */
1.1 maekawa 618:
619: r = mmLarger_tower3(f,g,gbList);
620: /* printf("mmLarger_tower3(%s,%s) --> %d\n",POLYToString(head(f),'*',1),POLYToString(head(g),'*',1),r); */
621: if (r == 2) { /* Now, compare by h */
622: if (f->m->e[0].D > g->m->e[0].D) return(1);
623: else if (f->m->e[0].D < g->m->e[0].D) return(0);
624: else return(2);
625: }else{
626: return(r);
627: }
628: }
629:
630: int mmLarger_tower3(POLY f,POLY g,struct object *gbList)
631: { /* gbList is assumed to be Slist */
632: int n,fv,gv,t,r,nn;
633: POLY fm;
634: POLY gm;
1.14 ! takayama 635: struct object gb = OINIT;
1.1 maekawa 636:
637: if (f == POLYNULL) {
638: if (g == POLYNULL) return(2);
639: else return(0);
640: }
641: if (g == POLYNULL) {
642: if (f == POLYNULL) return(2);
643: else return(1); /* It assumes the zero is the minimum element!! */
644: }
645: n = f->m->ringp->n;
646: nn = f->m->ringp->nn;
647: /* critical and modifiable */ /* m e_u > m e_v <==> m g_u > m g_v */
1.4 takayama 648: /* or equal and u < v */
1.1 maekawa 649: fv = f->m->e[nn].x ; /* extract component (vector) number of f! */
650: gv = g->m->e[nn].x ;
651: if (fv == gv) { /* They have the same component number. */
652: return(mmLarger_matrix_schreyer(f,g));
653: }
654:
655: if (gbList == NULL) return(mmLarger_matrix_schreyer(f,g));
1.4 takayama 656: /* modifiable: mmLarger_qmatrix */
1.1 maekawa 657: if (gbList->tag != Slist) {
658: warningOrder("mmLarger_tower(): gbList must be in Slist.\n");
659: return(1);
660: }
661: if (klength(gbList) ==0) return(mmLarger_matrix(f,g));
1.4 takayama 662: /* modifiable: mmLarger_qmatrix */
1.1 maekawa 663: gb = car(gbList); /* each entry must be monomials */
664: if (gb.tag != Sarray) {
665: warningOrder("mmLarger_tower3(): car(gbList) must be an array.\n");
666: return(1);
667: }
668: t = getoaSize(gb);
669: if (t == 0) return(mmLarger_tower3(f,g,cdr(gbList)));
670:
671: fm = pmCopy(head(f)); fm->m->e[nn].x = 0; /* f is not modified. */
672: gm = pmCopy(head(g)); gm->m->e[nn].x = 0;
673: if (fv >= t || gv >= t) {
674: warningOrder("mmLarger_tower3(): incompatible input and gbList.\n");
675: printf("Length of gb is %d, f is %s, g is %s\n",t,KPOLYToString(f),
1.4 takayama 676: KPOLYToString(g));
1.3 takayama 677: KSexecuteString(" show_ring ");
1.1 maekawa 678: return(1);
679: }
680: /* mpMult_poly is too expensive to call. @@@*/
681: r = mmLarger_tower3(mpMult_poly(fm,KopPOLY(getoa(gb,fv))),
682: mpMult_poly(gm,KopPOLY(getoa(gb,gv))),
683: cdr(gbList));
684: if (r != 2) return(r);
685: else if (fv == gv) return(2);
686: else if (fv > gv) return(0); /* modifiable */
687: else if (fv < gv) return(1); /* modifiable */
688: }
1.11 takayama 689:
690: static struct object auxPruneZeroRow(struct object ob) {
691: int i,m,size;
1.14 ! takayama 692: struct object obt = OINIT;
! 693: struct object rob = OINIT;
1.11 takayama 694: m = getoaSize(ob);
695: size=0;
696: for (i=0; i<m; i++) {
697: obt = getoa(ob,i);
698: if (getoaSize(obt) != 0) size++;
699: }
700: if (size == m) return ob;
701: rob = newObjectArray(size);
702: for (i=0, size=0; i<m; i++) {
703: obt = getoa(ob,i);
704: if (getoaSize(obt) != 0) {
705: putoa(rob,size,obt); size++;
706: }
707: }
708: return rob;
709: }
1.12 takayama 710: static struct object oRingToOXringStructure_long(struct ring *ringp)
1.10 takayama 711: {
1.14 ! takayama 712: struct object rob = OINIT;
! 713: struct object ob2 = OINIT;
! 714: struct object obMat = OINIT;
! 715: struct object obV = OINIT;
! 716: struct object obShift = OINIT;
! 717: struct object obt = OINIT;
1.10 takayama 718: char **TransX; char **TransD;
719: int n,i,j,m,p,nonzero;
720: int *om;
721: n = ringp->n;
722: m = ringp->orderMatrixSize;
723: om = ringp->order;
724: TransX = ringp->x; TransD = ringp->D;
725: if (m<=0) m = 1;
726: /*test: (1). getRing /rr set rr (oxRingStructure) dc */
727: obMat = newObjectArray(m);
728: for (i=0; i<m; i++) {
729: nonzero = 0;
730: for (j=0; j<2*n; j++) {
731: if (om[2*n*i+j] != 0) nonzero++;
732: }
733: ob2 = newObjectArray(nonzero*2);
734: nonzero=0;
735: for (j=0; j<2*n; j++) {
736: /* fprintf(stderr,"%d, ",nonzero); */
737: if (om[2*n*i+j] != 0) {
738: if (j < n) {
739: putoa(ob2,nonzero,KpoString(TransX[n-1-j])); nonzero++;
740: }else{
741: putoa(ob2,nonzero,KpoString(TransD[n-1-(j-n)])); nonzero++;
742: }
743: putoa(ob2,nonzero,KpoUniversalNumber(newUniversalNumber(om[2*n*i+j]))); nonzero++;
744: }
745: }
746: /* printObject(ob2,0,stderr); fprintf(stderr,".\n"); */
747: putoa(obMat,i,ob2);
748: }
1.11 takayama 749: obMat = auxPruneZeroRow(obMat);
1.10 takayama 750: /* printObject(obMat,0,stderr); */
751:
752: obV = newObjectArray(2*n);
753: for (i=0; i<n; i++) putoa(obV,i,KpoString(TransX[n-1-i]));
754: for (i=0; i<n; i++) putoa(obV,i+n,KpoString(TransD[n-1-i]));
755: /* printObject(obV,0,stderr); */
756:
757: if (ringp->degreeShiftSize) {
758: /*test:
759: [(x) ring_of_differential_operators [[(x)]] weight_vector 0
760: [(weightedHomogenization) 1 (degreeShift) [[1 2 1]]] ] define_ring ;
761: (1). getRing /rr set rr (oxRingStructure) dc message
762: */
763: obShift = newObjectArray(ringp->degreeShiftN);
764: for (i=0; i<ringp->degreeShiftN; i++) {
765: obt = newObjectArray(ringp->degreeShiftSize);
766: for (j=0; j< ringp->degreeShiftSize; j++) {
767: putoa(obt,j,KpoUniversalNumber(newUniversalNumber(ringp->degreeShift[i*(ringp->degreeShiftSize)+j])));
768: }
769: putoa(obShift,i,obt);
770: }
771: /* printObject(obShift,0,stderr); */
772: }
773:
774: p = 0;
775: if (ringp->degreeShiftSize) {
776: rob = newObjectArray(3);
777: obt = newObjectArray(2);
778: putoa(obt,0,KpoString("degreeShift"));
779: putoa(obt,1,obShift);
780: putoa(rob,p, obt); p++;
781: }else {
782: rob = newObjectArray(2);
783: }
784:
785: obt = newObjectArray(2);
786: putoa(obt,0,KpoString("v"));
787: putoa(obt,1,obV);
788: putoa(rob,p, obt); p++;
789:
790: obt = newObjectArray(2);
791: putoa(obt,0,KpoString("order"));
792: putoa(obt,1,obMat);
793: putoa(rob,p, obt); p++;
794:
1.12 takayama 795: return(rob);
796: }
797: static int auxEffectiveVar(int idx,int n) {
798: int x;
799: if (idx < n) x=1; else x=0;
800: if (x) {
801: if ((idx >= 1) && (idx < n-1)) return 1;
802: else return 0;
803: }else{
804: if ( 1 <= idx-n ) return 1;
805: else return 0;
806: }
807: }
808: /*test:
809: [(x,y) ring_of_differential_operators [[(Dx) 1 (Dy) 1]]
810: weight_vector 0] define_ring
811: (x). getRing (oxRingStructure) dc ::
812: */
813: static struct object oRingToOXringStructure_short(struct ring *ringp)
814: {
1.14 ! takayama 815: struct object rob = OINIT;
! 816: struct object ob2 = OINIT;
! 817: struct object obMat = OINIT;
! 818: struct object obV = OINIT;
! 819: struct object obShift = OINIT;
! 820: struct object obt = OINIT;
1.12 takayama 821: char **TransX; char **TransD;
822: int n,i,j,m,p,nonzero;
823: int *om;
824: n = ringp->n;
825: m = ringp->orderMatrixSize;
826: om = ringp->order;
827: TransX = ringp->x; TransD = ringp->D;
828: if (m<=0) m = 1;
829: /*test: (1). getRing /rr set rr (oxRingStructure) dc */
830: obMat = newObjectArray(m);
831: for (i=0; i<m; i++) {
832: nonzero = 0;
833: for (j=0; j<2*n; j++) {
834: if ((om[2*n*i+j] != 0) && auxEffectiveVar(j,n)) nonzero++;
835: }
836: ob2 = newObjectArray(nonzero*2);
837: nonzero=0;
838: for (j=0; j<2*n; j++) {
839: /* fprintf(stderr,"%d, ",nonzero); */
840: if ((om[2*n*i+j] != 0) && auxEffectiveVar(j,n)) {
841: if (j < n) {
842: putoa(ob2,nonzero,KpoString(TransX[n-1-j])); nonzero++;
843: }else{
844: putoa(ob2,nonzero,KpoString(TransD[n-1-(j-n)])); nonzero++;
845: }
846: putoa(ob2,nonzero,KpoUniversalNumber(newUniversalNumber(om[2*n*i+j]))); nonzero++;
847: }
848: }
849: /* printObject(ob2,0,stderr); fprintf(stderr,".\n"); */
850: putoa(obMat,i,ob2);
851: }
852: obMat = auxPruneZeroRow(obMat);
853: /* printObject(obMat,0,stderr); */
854:
855: obV = newObjectArray(2*n-3);
856: for (i=0; i<n-2; i++) putoa(obV,i,KpoString(TransX[n-1-i-1]));
857: for (i=0; i<n-1; i++) putoa(obV,i+n-2,KpoString(TransD[n-1-i-1]));
858: /* printObject(obV,0,stderr); */
859:
860: if (ringp->degreeShiftSize) {
861: /*test:
862: [(x) ring_of_differential_operators [[(x)]] weight_vector 0
863: [(weightedHomogenization) 1 (degreeShift) [[1 2 1]]] ] define_ring ;
864: (1). getRing /rr set rr (oxRingStructure) dc message
865: */
866: obShift = newObjectArray(ringp->degreeShiftN);
867: for (i=0; i<ringp->degreeShiftN; i++) {
868: obt = newObjectArray(ringp->degreeShiftSize);
869: for (j=0; j< ringp->degreeShiftSize; j++) {
870: putoa(obt,j,KpoUniversalNumber(newUniversalNumber(ringp->degreeShift[i*(ringp->degreeShiftSize)+j])));
871: }
872: putoa(obShift,i,obt);
873: }
874: /* printObject(obShift,0,stderr); */
875: }
876:
877: p = 0;
878: if (ringp->degreeShiftSize) {
879: rob = newObjectArray(3);
880: obt = newObjectArray(2);
881: putoa(obt,0,KpoString("degreeShift"));
882: putoa(obt,1,obShift);
883: putoa(rob,p, obt); p++;
884: }else {
885: rob = newObjectArray(2);
886: }
887:
888: obt = newObjectArray(2);
889: putoa(obt,0,KpoString("v"));
890: putoa(obt,1,obV);
891: putoa(rob,p, obt); p++;
892:
893: obt = newObjectArray(2);
894: putoa(obt,0,KpoString("order"));
895: putoa(obt,1,obMat);
896: putoa(rob,p, obt); p++;
897:
898: return(rob);
899: }
900: struct object oRingToOXringStructure(struct ring *ringp)
901: {
1.14 ! takayama 902: struct object rob = OINIT;
! 903: struct object tob = OINIT;
1.12 takayama 904: rob = newObjectArray(2);
905: tob = oRingToOXringStructure_short(ringp);
906: putoa(rob,0,tob);
907: tob = oRingToOXringStructure_long(ringp);
908: putoa(rob,1,tob);
1.10 takayama 909: return(rob);
910: }
911:
1.1 maekawa 912: static void warningOrder(s)
1.4 takayama 913: char *s;
1.1 maekawa 914: {
915: fprintf(stderr,"Warning in order.c: %s\n",s);
916: }
917:
918: static void errorOrder(s)
1.4 takayama 919: char *s;
1.1 maekawa 920: {
921: fprintf(stderr,"order.c: %s\n",s);
922: exit(14);
923: }
924:
925:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>