Annotation of OpenXM/src/kan96xx/Kan/poly2.c, Revision 1.8
1.8 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/poly2.c,v 1.7 2005/07/03 11:08:54 ohara Exp $ */
1.1 maekawa 2: #include <stdio.h>
1.7 ohara 3: #include <stdlib.h>
1.1 maekawa 4: #include "datatype.h"
5: #include "stackm.h"
6: #include "extern.h"
7: #include "extern2.h"
8:
9: static POLY mapZmonom(POLY f,struct ring *ringp);
10:
11: POLY ppAdd(f,g)
1.3 takayama 12: POLY f; POLY g; /* The result is read only. */
1.1 maekawa 13: {
14: POLY node;
15: struct listPoly nod;
16: POLY h;
17: struct coeff *c;
18: int gt;
19:
20: node = &nod;
21: node->next = POLYNULL;
22: h = node;
23: if (f == POLYNULL) return(g);
24: if (g == POLYNULL) return(f);
25: checkRing(f,g);
26:
27: while (f != POLYNULL && g != POLYNULL) {
28: /*printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
29: checkRing2(f,g); /* for debug */
30: gt = (*mmLarger)(f,g);
31: switch (gt) {
32: case 1: /* f > g */
33: h -> next = newCell(f->coeffp,f->m);
34: h = h->next;
35: f = f->next;
36: if (f == POLYNULL) {
1.3 takayama 37: h->next = g;
38: return(node->next);
1.1 maekawa 39: }
40: break;
41: case 0: /* f < g */
42: h->next = newCell(g->coeffp,g->m);
43: h = h->next;
44: g = g->next;
45: if (g == POLYNULL) {
1.3 takayama 46: h->next = f;
47: return(node->next);
1.1 maekawa 48: }
49: break;
50: case 2:/* f == g */
51: c = coeffCopy(f->coeffp);
52: Cadd(c,c,g->coeffp);
53: if (!isZero(c)) {
1.3 takayama 54: h->next = newCell(c,f->m);
55: h = h->next;
56: f = f->next;
57: g = g->next;
58: if (f == POLYNULL) {
59: h->next = g;
60: return(node->next);
61: }
62: if (g == POLYNULL) {
63: h->next = f;
64: return(node->next);
65: }
1.1 maekawa 66: }else{
1.3 takayama 67: f = f->next;
68: g = g->next;
69: if (f == POLYNULL) {
70: h->next = g;
71: return(node->next);
72: }
73: if (g == POLYNULL) {
74: h->next = f;
75: return(node->next);
76: }
1.1 maekawa 77: }
78: break;
79: default:
80: errorPoly("ppAdd(). Internal error. Invalid value of mmLarger.");
81: break;
82: }
83: }
84: return(node->next);
85: }
86:
87: POLY ppSub(f,g)
1.3 takayama 88: POLY f; POLY g; /* The result is read only. */
1.1 maekawa 89: {
90: POLY h;
91: struct coeff *c;
92:
93: if (g == POLYNULL) return(f);
94: if (f == POLYNULL) return(cpMult(intToCoeff(-1,g->m->ringp),g));
95: checkRing(f,g);
96:
97: h = cpMult(intToCoeff(-1,g->m->ringp),g);
98: return(ppAdd(f,h));
99: }
100:
101:
102: POLY cpMult(c,f)
1.3 takayama 103: struct coeff *c;
104: POLY f;
1.1 maekawa 105: {
106: POLY node;
107: struct listPoly nod;
108: POLY h;
109: struct coeff *newc;
110: int p;
111: node = &nod;
112: if (f == POLYNULL || isZero(c)) return(POLYNULL);
113: p = f->m->ringp->p;
114: node ->next = POLYNULL;
115: h = node;
116: while (f != POLYNULL) {
117: newc = coeffCopy(c);
118: Cmult(newc,newc,f->coeffp);
119: if ((p==0) || !isZero(newc)) {
120: h->next = newCell(newc,f->m);
121: h = h->next;
122: }
123: f = f->next;
124: }
125: return(node->next);
126: }
127:
128: MONOMIAL monomialAdd_poly(m,m2)
1.3 takayama 129: MONOMIAL m,m2;
1.1 maekawa 130: {
131: extern int Msize;
132: MONOMIAL f;
133: int i;
134: int n;
135: n = m->ringp->n;
136: f = (MONOMIAL) sGC_malloc(sizeof(struct smallMonomial)+n*Msize);
137:
138: if (f == (MONOMIAL) NULL) errorPoly("No more memory.");
139: f->ringp = m->ringp;
140: for (i=0; i<n; i++) {
141: (f->e)[i].x = (m->e)[i].x + (m2->e)[i].x;
142: (f->e)[i].D = (m->e)[i].D + (m2->e)[i].D;
143: }
144: return(f);
145: }
146:
147: /* Note that mpMult_poly is called from mmLarger_tower! */
148: POLY mpMult_poly(f,g)
1.3 takayama 149: POLY f;
150: POLY g;
1.1 maekawa 151: {
152: POLY node;
153: struct listPoly nod;
154: struct coeff *c;
155: int n,i;
156: POLY h;
157: MONOMIAL m;
158: int p;
159: node = &nod;
160: if (f == POLYNULL || g == POLYNULL) return(POLYNULL);
161: node->next = POLYNULL;
162: h = node;
163: checkRing(f,g);
164: n = f->m->ringp->n; p = f->m->ringp->p;
165: while(g != POLYNULL) {
166: checkRing2(f,g);
167: c = coeffCopy(f->coeffp);
168: Cmult(c,c,g->coeffp);
169: if ((p==0) || !isZero(c)) {
170: m = (*monomialAdd)(f->m,g->m);
171: h->next = newCell(c,m);
172: h = h->next;
173: }
174: g = g->next;
175: }
176: return(node->next);
177: }
178:
179: POLY ppMult_old(f,g)
1.3 takayama 180: POLY f,g;
1.1 maekawa 181: {
182: POLY r;
183: POLY tmp;
184: r = POLYNULL;
185: while( f != POLYNULL) {
186: tmp = (*mpMult)(f,g);
187: r = ppAddv(r,tmp); /* r and tmp will be broken */
188: f = f->next;
189: }
190: return(r);
191: }
192:
193: POLY ppAddv(f,g)
1.3 takayama 194: POLY f; POLY g; /* It breaks f and g. Use it just after calling mpMult() */
1.1 maekawa 195: {
196: POLY node;
197: struct listPoly nod;
198: POLY h;
199: struct coeff *c;
200: int gt;
201:
202: node = &nod;
203: /*printf("ppAddv:1%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
204: node->next = POLYNULL;
205: h = node;
206: if (f == POLYNULL) return(g);
207: if (g == POLYNULL) return(f);
208: checkRing(f,g);
209:
210: while (f != POLYNULL && g != POLYNULL) {
211: checkRing2(f,g); /* for debug */
212: /*printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1));*/
213: gt = (*mmLarger)(f,g);
214: switch (gt) {
215: case 1: /* f > g */
216: h->next = f;
217: h = h->next; f = f->next;;
218: if (f == POLYNULL) {
1.3 takayama 219: h->next = g;
220: return(node->next);
1.1 maekawa 221: }
222: break;
223: case 0: /* f < g */
224: h->next = g;
225: h = h->next; g = g->next;
226: if (g == POLYNULL) {
1.3 takayama 227: h->next = f;
228: return(node->next);
1.1 maekawa 229: }
230: break;
231: case 2:/* f == g */
232: c = f->coeffp;
233: Cadd(c,c,g->coeffp);
234: if (!isZero(c)) {
1.3 takayama 235: h->next = f;
236: h = h->next; f = f->next;;
237: g = g->next;
238: if (f == POLYNULL) {
239: h->next = g;
240: return(node->next);
241: }
242: if (g == POLYNULL) {
243: h->next = f;
244: return(node->next);
245: }
1.1 maekawa 246: }else{
1.3 takayama 247: f = f->next;
248: g = g->next;
249: if (f == POLYNULL) {
250: h->next = g;
251: return(node->next);
252: }
253: if (g == POLYNULL) {
254: h->next = f;
255: return(node->next);
256: }
1.1 maekawa 257: }
258: break;
259: default:
260: errorPoly("ppAddv(). Internal error. Invalid value of mmLarger.");
261: break;
262: }
263: }
264: return(node->next);
265: }
266:
267: POLY pPower(f,k)
1.3 takayama 268: POLY f;
269: int k;
1.1 maekawa 270: {
271: POLY r;
272: int i,n;
273: if (f == POLYNULL) return(POLYNULL); /* Is it ok? 0^0 = 0.*/
274: if (k == 0) return(cxx(1,0,0,f->m->ringp));
275: if (f->next == POLYNULL && k<0) {
276: /* when f is monomial. */
277: r = newCell(coeffCopy(f->coeffp),monomialCopy(f->m));
278: n = r->m->ringp->n;
279: for (i=0; i<n; i++) {
280: r->m->e[i].x *= k;
281: r->m->e[i].D *= k;
282: }
283: if (!isOne(r->coeffp)) {
284: warningPoly("pPower(poly,negative integer) not implemented yet. Returns 1.");
285: r = cxx(1,0,0,f->m->ringp);
286: }
287: return(r);
288: }
289: r = cxx(1,0,0,f->m->ringp);
290: if (k < 0) {
291: warningPoly("pPower(poly,negative integer) not implemented yet. Returns 1.");
292: }
293: for (i=0; i<k; i++) {
294: r = ppMult(f,r);
295: }
296: return(r);
297: }
298:
299: POLY pPower_poly(f,k)
1.3 takayama 300: POLY f;
301: int k;
1.1 maekawa 302: {
303: POLY r;
304: int i,n;
305: if (f == POLYNULL) return(POLYNULL); /* Is it ok? 0^0 = 0.*/
306: if (k == 0) return(cxx(1,0,0,f->m->ringp));
307: if (f->next == POLYNULL && k<0) {
308: /* when f is monomial. */
309: r = newCell(coeffCopy(f->coeffp),monomialCopy(f->m));
310: n = r->m->ringp->n;
311: for (i=0; i<n; i++) {
312: r->m->e[i].x *= k;
313: r->m->e[i].D *= k;
314: }
315: if (!isOne(r->coeffp)) {
316: warningPoly("pPower_poly(poly,negative integer) not implemented yet. Returns 1.");
317: r = cxx(1,0,0,f->m->ringp);
318: }
319: return(r);
320: }
321: r = cxx(1,0,0,f->m->ringp);
322: if (k < 0) {
323: warningPoly("pPower_poly(poly,negative integer) not implemented yet. Returns 1.");
324: }
325: for (i=0; i<k; i++) {
326: r = ppMult_poly(f,r);
327: }
328: return(r);
329: }
330:
331: POLY modulop_trash(f,ringp)
1.3 takayama 332: POLY f;
333: struct ring *ringp;
1.1 maekawa 334: {
335: int p;
336: POLY h;
337: MP_INT *c;
338: int cc;
339: POLY node;
340: struct ring *nextRing;
341: POLY fc;
342:
343: if (f == POLYNULL) return(f);
344: p = ringp->p;
345: if (f->m->ringp->p != 0) {
346: warningPoly("modulop(f,ringp) must be called with f in the characteristic 0 ring. Returns 0.");
347: return(POLYNULL);
348: }
349: if (f->m->ringp->n != ringp->n) {
350: warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
351: return(POLYNULL);
352: }
353:
354: /* The case of ringp->next != NULL */
355: if (ringp->next != (struct ring *)NULL) {
356: nextRing = ringp->next;
357: node = newCell(newCoeff(),newMonomial(ringp));
358: node->next = POLYNULL;
359: h = node;
360:
361: while (f != POLYNULL) {
362: fc = bxx(f->coeffp->val.bigp,0,0,nextRing);
363: h->next = newCell(newCoeff(),monomialCopy(f->m));
364: h = h->next;
365: h->m->ringp = ringp;
366: h->coeffp->tag = POLY_COEFF;
367: h->coeffp->p = p;
368: h->coeffp->val.f = fc;
369: f = f->next;
370: }
371: return(node->next);
372: }
373:
374:
375: /* In case of ringp->next == NULL */
376: if (p != 0) {
377: c = newMP_INT();
378: node = newCell(newCoeff(),newMonomial(ringp));
379: node->next = POLYNULL;
380: h = node;
381:
382: while (f != POLYNULL) {
383: mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
384: cc = (int) mpz_get_si(c);
385: if (cc != 0) {
1.3 takayama 386: h->next = newCell(newCoeff(),monomialCopy(f->m));
387: h = h->next;
388: h->m->ringp = ringp;
389: h->coeffp->tag = INTEGER;
390: h->coeffp->p = p;
391: h->coeffp->val.i = cc;
1.1 maekawa 392: }
393: f = f->next;
394: }
395: return(node->next);
396: }else{
397: h = f = pcmCopy(f);
398: while (f != POLYNULL) {
399: f->m->ringp = ringp;
400: f = f->next;
401: }
402: return(h);
403: }
404:
405: }
406:
407: POLY modulop(f,ringp)
1.3 takayama 408: POLY f;
409: struct ring *ringp;
410: /* Z[x] ---> R[x] where R=Z, Z/Zp, ringp->next. */
1.1 maekawa 411: {
412: int p;
413: POLY h;
414: MP_INT *c;
415: int cc;
416: POLY node;
417: POLY fc;
418:
419: if (f == POLYNULL) return(f);
420: p = ringp->p;
421: if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
422: warningPoly("modulop(f,ringp) must be called with f in the characteristic 0 ring Z[x]. Returns 0.");
423: return(POLYNULL);
424: }
425: if (f->m->ringp->n != ringp->n) {
426: warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
427: return(POLYNULL);
428: }
429:
430: /* [1] The case of R = ringp->next */
431: if (ringp->next != (struct ring *)NULL) {
432: h = ZERO;
433: while (f != POLYNULL) {
434: h = ppAdd(h,mapZmonom(f,ringp));
435: f = f->next;
436: }
437: return(h);
438: }
439:
440: /* [2] The case of R = Z/Zp */
441: if (p != 0) {
442: c = newMP_INT();
443: node = newCell(newCoeff(),newMonomial(ringp));
444: node->next = POLYNULL;
445: h = node;
446:
447: while (f != POLYNULL) {
448: mpz_mod_ui(c,f->coeffp->val.bigp,(unsigned long int)p);
449: cc = (int) mpz_get_si(c);
450: if (cc != 0) {
1.3 takayama 451: h->next = newCell(newCoeff(),monomialCopy(f->m));
452: h = h->next;
453: h->m->ringp = ringp;
454: h->coeffp->tag = INTEGER;
455: h->coeffp->p = p;
456: h->coeffp->val.i = cc;
1.1 maekawa 457: }
458: f = f->next;
459: }
460: return(node->next);
461: }
462:
463: /* [3] The case of R = Z */
464: h = f = pcmCopy(f);
465: while (f != POLYNULL) {
466: f->m->ringp = ringp;
467: f = f->next;
468: }
469: return(h);
470:
471:
472: }
473:
474: POLY modulopZ(f,pcoeff)
1.3 takayama 475: POLY f;
476: struct coeff *pcoeff;
477: /* Z[x] ---> Z[x] , f ---> f mod pcoeff*/
1.1 maekawa 478: {
479: int p;
480: POLY h;
481: struct coeff *c;
482: int cc;
483: POLY node;
484: POLY fc;
485: MP_INT *bigp;
486: MP_INT *tmp;
487: struct ring *ringp;
488:
489: if (f == POLYNULL) return(f);
490: ringp = f->m->ringp;
491: if (pcoeff->tag != MP_INTEGER) {
492: warningPoly("modulopZ(): pcoeff must be a universalNumber.");
493: warningPoly("Returns 0.");
494: return(POLYNULL);
495: }
496: bigp = pcoeff->val.bigp;
497: if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
498: warningPoly("modulopZ(f,p) must be called with f in the characteristic 0 ring Z[x]. Returns 0.");
499: return(POLYNULL);
500: }
501: if (f->m->ringp->n != ringp->n) {
502: warningPoly("modulop(f,ringp): f->m->ringp->n must be equal to ringp->n. Returns 0.");
503: return(POLYNULL);
504: }
505:
506: /* [1] The case of R = ringp->next */
507: if (ringp->next != (struct ring *)NULL) {
508: warningPoly("modulopZ workds only for flat polynomials. Returns 0.");
509: return(POLYNULL);
510: }
511:
512: /* [2] The case of R = Z */
513: node = newCell(newCoeff(),newMonomial(ringp));
514: node->next = POLYNULL;
515: h = node;
516:
517: c = newCoeff();
518: tmp = newMP_INT();
519: while (f != POLYNULL) {
520: mpz_mod(tmp,f->coeffp->val.bigp,bigp);
521: if (mpz_sgn(tmp) != 0) {
522: c->tag = MP_INTEGER;
523: c->p = 0;
524: c->val.bigp = tmp;
525: h->next = newCell(c,monomialCopy(f->m));
526: h = h->next;
527: h->m->ringp = ringp;
528:
529: c = newCoeff();
530: tmp = newMP_INT();
531: }
532: f = f->next;
533: }
534: return(node->next);
535:
536: }
537:
538: struct pairOfPOLY quotientByNumber(f,pcoeff)
1.3 takayama 539: POLY f;
540: struct coeff *pcoeff;
541: /* Z[x] ---> Z[x],Z[x] , f = first*pcoeff + second */
1.1 maekawa 542: {
543: int p;
544: POLY h;
545: POLY h2;
546: struct coeff *c;
547: int cc;
548: POLY node;
549: struct coeff *c2;
550: POLY node2;
551: POLY fc;
552: MP_INT *bigp;
553: MP_INT *tmp;
554: MP_INT *tmp2;
555: struct ring *ringp;
556: struct pairOfPOLY r;
557:
558: if (f == POLYNULL) {
559: r.first = f; r.second = f;
560: return(r);
561: }
562: ringp = f->m->ringp;
563: if (pcoeff->tag != MP_INTEGER) {
564: warningPoly("quotientByNumber(): pcoeff must be a universalNumber.");
565: warningPoly("Returns (0,0).");
566: r.first = f; r.second = f;
567: return(r);
568: }
569: bigp = pcoeff->val.bigp;
570: if (f->m->ringp->p != 0 || f->m->ringp->next != (struct ring *)NULL) {
571: warningPoly("quotientByNumber(f,p) must be called with f in the characteristic 0 ring Z[x]. Returns (0,0).");
572: r.first = f; r.second = f;
573: return(r);
574: }
575: if (f->m->ringp->n != ringp->n) {
576: warningPoly("quotientByNumber(f,p): f->m->ringp->n must be equal to ringp->n. Returns 0.");
577: r.first = f; r.second = f;
578: return(r);
579: }
580:
581: /* [1] The case of R = ringp->next */
582: if (ringp->next != (struct ring *)NULL) {
583: warningPoly("quotientByNumber() workds only for flat polynomials. Returns 0.");
584: r.first = f; r.second = f;
585: return(r);
586: }
587:
588: /* [2] The case of R = Z */
589: node = newCell(newCoeff(),newMonomial(ringp));
590: node->next = POLYNULL;
591: h = node;
592: node2 = newCell(newCoeff(),newMonomial(ringp));
593: node2->next = POLYNULL;
594: h2 = node2;
595:
596: c = newCoeff();
597: tmp = newMP_INT();
598: c2 = newCoeff();
599: tmp2 = newMP_INT();
600: while (f != POLYNULL) {
601: mpz_mod(tmp,f->coeffp->val.bigp,bigp);
602: if (mpz_sgn(tmp) != 0) {
603: c->tag = MP_INTEGER;
604: c->p = 0;
605: c->val.bigp = tmp;
606: h->next = newCell(c,monomialCopy(f->m));
607: h = h->next;
608: h->m->ringp = ringp;
609:
610: c = newCoeff();
611: tmp = newMP_INT();
612: }
613: mpz_tdiv_q(tmp2,f->coeffp->val.bigp,bigp);
614: if (mpz_sgn(tmp2) != 0) {
615: c2->tag = MP_INTEGER;
616: c2->p = 0;
617: c2->val.bigp = tmp2;
618: h2->next = newCell(c2,monomialCopy(f->m));
619: h2 = h2->next;
620: h2->m->ringp = ringp;
621:
622: c2 = newCoeff();
623: tmp2 = newMP_INT();
624: }
625: f = f->next;
626: }
627: r.first = node2->next;
628: r.second = node->next;
629: return(r);
630:
631: }
632:
633:
634: POLY modulo0(f,ringp)
1.3 takayama 635: POLY f;
636: struct ring *ringp;
1.1 maekawa 637: {
638: int p;
639: POLY h;
640: struct coeff *c;
641: POLY node;
642: if (f == POLYNULL) return(f);
643: p = ringp->p;
644: if (p != 0) {
645: warningPoly("modulo0(f,ringp) must be called with the characteristic 0 ring*ringp. Returns 0.");
646: return(POLYNULL);
647: }
648: switch (f->coeffp->tag) {
649: case MP_INTEGER:
650: if (f->m->ringp->p == 0) {
651: node = pcmCopy(f);
652: f = node;
653: while (f != POLYNULL) {
1.3 takayama 654: f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
655: f = f->next;
1.1 maekawa 656: }
657: return(node);
658: }
659: break;
660: case POLY_COEFF:
661: node = pcmCopy(f);
662: f = node;
663: while (f != POLYNULL) {
664: f->m->ringp = ringp; /* Touch the monomial "ringp" field. */
665: f = f->next;
666: }
667: return(node);
668: break;
669: case INTEGER:
670: node = newCell(newCoeff(),newMonomial(ringp));
671: node->next = POLYNULL;
672: h = node;
673:
674: while (f != POLYNULL) {
675: c = newCoeff();
676: c->tag = MP_INTEGER;
677: c->p = 0;
678: c->val.bigp = newMP_INT();
679: mpz_set_si(c->val.bigp,f->coeffp->val.i);
680: h->next = newCell(c,monomialCopy(f->m));
681: h = h->next;
682: h->m->ringp = ringp;
683: f = f->next;
684: }
685: return(node->next);
686: break;
687: default:
688: warningPoly("modulo0(): coefficients have to be MP_INTEGER or INTEGER. Returns 0");
689: return(POLYNULL);
690: break;
691: }
692: }
693:
694:
695: struct object test(ob) /* test3 */
1.3 takayama 696: struct object ob;
1.1 maekawa 697: {
1.6 takayama 698: struct object rob = OINIT;
1.1 maekawa 699: int k;
700: static POLY f0;
701: static POLY f1;
702:
703: POLY addNode;
704: POLY f;
705: int i;
706: static int s=0;
707: MP_INT *mp;
708: extern struct ring *SmallRingp;
709: extern struct ring *CurrentRingp;
710: addNode = pMalloc(SmallRingp);
711: k = ob.lc.ival;
712: switch(s) {
713: case 0:
714: f0 = addNode;
715: for (i=k; i>=0; i--) {
716: f0->next = bxx(BiiPower(-k,i),0,i,CurrentRingp);
717: if (f0->next != POLYNULL) {
1.3 takayama 718: f0 = f0->next;
1.1 maekawa 719: }
720: }
721: f0 = addNode->next;
722: s++;
723: rob.lc.poly = f0;
724: break;
725: case 1:
726: f1 = addNode;
727: for (i=k; i>=0; i--) {
728: f1->next = bxx(BiiPower(k,i),0,i,CurrentRingp);
729: if (f1->next != POLYNULL) {
1.3 takayama 730: f1 = f1->next;
1.1 maekawa 731: }
732: }
733: f1 = addNode->next;
734: s = 0;
735: rob.lc.poly = f1;
736: break;
737: default:
738: rob.lc.poly = POLYNULL;
739: s = 0;
740: break;
741: }
742:
743:
744: rob.tag = Spoly;
745: return(rob);
746: }
747:
748:
749: int pLength(f)
1.3 takayama 750: POLY f;
1.1 maekawa 751: {
752: int c=0;
753: if (f ISZERO) return(0);
754: while (f != POLYNULL) {
755: c++;
756: f = f->next;
757: }
758: return(c);
759: }
760:
761:
762: POLY ppAddv2(f,g,top,nexttop)
1.3 takayama 763: POLY f; POLY g; /* It breaks f and g. Use it just after calling mpMult() */
764: POLY top;
765: POLY *nexttop;
766: /* top is the starting address in the list f.
767: if top == POLYNULL, start from f.
768:
769: *nexttop == 0
770: == g
771: == h or 0
1.1 maekawa 772:
1.3 takayama 773: It must be called as r = ppAddv2(r,g,...);
1.1 maekawa 774: */
775: {
776: POLY node;
777: struct listPoly nod;
778: POLY h;
779: struct coeff *c;
780: int gt;
781: POLY g0;
782: POLY f0; /* for debug */
783:
784: node = &nod;
785: /* printf("ppAddv:1%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
786: node->next = POLYNULL;
787: h = node;
788: *nexttop = POLYNULL;
789: if (f == POLYNULL) return(g);
790: if (g == POLYNULL) return(f);
791: checkRing(f,g);
792:
793: f0 = f;
794: if (top != POLYNULL) {
795: while (f != top) {
796: if (f == POLYNULL) {
1.3 takayama 797: fprintf(stderr,"\nppAddv2(): Internal error.\n");fflush(stderr);
798: fprintf(stderr,"f = %s\n",POLYToString(f0,'*',0));
799: fprintf(stderr,"g = %s\n",POLYToString(g0,'*',0));
800: fprintf(stderr,"top=%s\n",POLYToString(top,'*',0));
801: errorPoly("ppAddv2(). Internal error=1.");
1.1 maekawa 802: }
803: h->next = f;
804: h = h->next;
805: f = f->next;
806: }
807: }
808: g0 = g;
809: *nexttop = g0;
810:
811: while (f != POLYNULL && g != POLYNULL) {
812: checkRing2(f,g); /* for debug */
813: /* printf("%s + %s\n",POLYToString(f,'*',1),POLYToString(g,'*',1)); */
814: gt = (*mmLarger)(f,g);
815: switch (gt) {
816: case 1: /* f > g */
817: h->next = f;
818: h = h->next; f = f->next;;
819: if (f == POLYNULL) {
1.3 takayama 820: h->next = g;
821: return(node->next);
1.1 maekawa 822: }
823: break;
824: case 0: /* f < g */
825: h->next = g;
826: h = h->next; g = g->next;
827: if (g == POLYNULL) {
1.3 takayama 828: h->next = f;
829: return(node->next);
1.1 maekawa 830: }
831: break;
832: case 2:/* f == g */
833: c = g->coeffp;
834: Cadd(c,f->coeffp,c);
835: if (!isZero(c)) {
1.3 takayama 836: h->next = g;
837: h = h->next;
838: f = f->next;;
839: g = g->next;
840: if (f == POLYNULL) {
841: h->next = g;
842: return(node->next);
843: }
844: if (g == POLYNULL) {
845: h->next = f;
846: return(node->next);
847: }
1.1 maekawa 848: }else{
1.3 takayama 849: if (g == g0) {
850: if (h != node) {
851: *nexttop = h;
852: }else{
853: *nexttop = POLYNULL;
854: }
855: }
1.1 maekawa 856:
1.3 takayama 857: f = f->next;
858: g = g->next;
1.1 maekawa 859:
1.3 takayama 860: if (f == POLYNULL) {
861: h->next = g;
862: return(node->next);
863: }
864: if (g == POLYNULL) {
865: h->next = f;
866: return(node->next);
867: }
1.1 maekawa 868: }
869: break;
870: default:
871: errorPoly("ppAddv(). Internal error. Invalid value of mmLarger.");
872: break;
873: }
874: }
875: return(node->next);
876: }
877:
878: POLY ppMult(f,g)
1.3 takayama 879: POLY f,g;
1.1 maekawa 880: {
881: POLY r;
882: POLY tmp;
883: POLY top;
884: POLY nexttop;
885: r = POLYNULL; top = POLYNULL;
886: while( f != POLYNULL) {
887: /* tmp = (*mpMult)(f,g); (*mpMult) is no more used. */
888: tmp = (*(f->m->ringp->multiplication))(f,g);
889: /*printf("mpMult(%s,%s) ->%s\n",POLYToString(f,'*',1),POLYToString(g,'*',1),POLYToString(tmp,'*',1)); */
890: r = ppAddv2(r,tmp,top,&nexttop); /* r and tmp will be broken */
891: top = nexttop;
892: f = f->next;
893: }
894: return(r);
895: }
896:
897: POLY ppMult_poly(f,g)
1.3 takayama 898: POLY f,g;
1.1 maekawa 899: {
900: POLY r;
901: POLY tmp;
902: POLY top;
903: POLY nexttop;
904: r = POLYNULL; top = POLYNULL;
905: while( f != POLYNULL) {
906: tmp = mpMult_poly(f,g);
907: r = ppAddv2(r,tmp,top,&nexttop); /* r and tmp will be broken */
908: top = nexttop;
909: f = f->next;
910: }
911: return(r);
912: }
913:
914: POLY mapZmonom(f,ringp)
1.3 takayama 915: POLY f; /* assumes monomial. f \in Z[x] */
916: struct ring *ringp; /* R[x] */
1.1 maekawa 917: {
918: struct ring *nextRing;
919: struct ring nextRing0;
920: POLY ff;
921: POLY node;
922: POLY gg;
923: int l,c,d;
924:
925: nextRing = ringp->next;
926: nextRing0 = *nextRing; nextRing0.p = 0; nextRing0.next = 0;
927: /* nextRing0 = Z[y] where y is the variables of R. */
928:
929: ff = bxx(f->coeffp->val.bigp,0,0,&nextRing0);
930: ff = modulop(ff,nextRing);
931:
932: node = newCell(newCoeff(),monomialCopy(f->m));
933: node->next = POLYNULL;
934: node->m->ringp = ringp;
935:
936: node->coeffp->p = ringp->p;
937:
938: l = ringp->l; c = ringp->c;
939: /* If q-analog q x ---> (q) x. */
940: if (l-c > 0) {
941: d = node->m->e[0].x; /* degree of q in R[x].*/
942: node->m->e[0].x = 0;
943: gg = cxx(1,0,d,nextRing); /* q^d = x[0]^d */
944: ff = ppMult(gg,ff);
945: }
946:
947: node->coeffp->tag = POLY_COEFF;
948: node->coeffp->val.f = ff;
949: return(node);
950: }
1.4 takayama 951:
952: POLY reduceContentOfPoly(POLY f,struct coeff **contp) {
953: struct coeff *cont;
954: struct coeff *cOne = NULL;
955: extern struct ring *SmallRingp;
956: if (cOne == NULL) cOne = intToCoeff(1,SmallRingp);
1.5 takayama 957: *contp = cOne;
1.4 takayama 958:
959: if (f == POLYNULL) return f;
960: if (f->m->ringp->p != 0) return f;
961: if (f->coeffp->tag != MP_INTEGER) return f;
962: cont = gcdOfCoeff(f);
963: *contp = cont;
964: if (coeffGreater(cont,cOne)) {
965: f = quotientByNumber(f,cont).first;
966: }
967: return f;
968: }
969:
970: int coeffSizeMin(POLY f) {
971: int size;
972: int t;
973: if (f == POLYNULL) return 0;
974: if (f->m->ringp->p != 0) return 0;
975: if (f->coeffp->tag != MP_INTEGER) return 0;
976: size = mpz_size(f->coeffp->val.bigp);
977: while (f != POLYNULL) {
978: t = mpz_size(f->coeffp->val.bigp);
979: if (t < size) size = t;
980: if (size == 1) return size;
981: f = f->next;
982: }
983: }
984:
985: struct coeff *gcdOfCoeff(POLY f) {
986: extern struct ring *SmallRingp;
987: struct coeff *t;
988: MP_INT *tmp;
989: MP_INT *tmp2;
990: static MP_INT *cOne = NULL;
991: if (cOne == NULL) {
992: cOne = newMP_INT();
993: mpz_set_si(cOne,(long) 1);
994: }
1.5 takayama 995: if (f == POLYNULL) return intToCoeff(0,SmallRingp);
1.4 takayama 996: if (f->m->ringp->p != 0) return intToCoeff(0,SmallRingp);
997: if (f->coeffp->tag != MP_INTEGER) return intToCoeff(0,SmallRingp);
998: tmp = f->coeffp->val.bigp;
999: tmp2 = newMP_INT();
1000: while (f != POLYNULL) {
1001: mpz_gcd(tmp2,tmp,f->coeffp->val.bigp); /* tmp = tmp2 OK? */
1002: tmp = tmp2;
1003: if (mpz_cmp(tmp,cOne)==0) return intToCoeff(1,SmallRingp);
1004: f = f->next;
1005: }
1006: return mpintToCoeff(tmp,SmallRingp);
1.1 maekawa 1007:
1.5 takayama 1008: }
1009:
1010: int shouldReduceContent(POLY f,int ss) {
1.8 ! takayama 1011: extern int DoCancel;
1.5 takayama 1012: static int prevSize = 1;
1013: int size;
1014: if (f == POLYNULL) return 0;
1015: if (f->m->ringp->p != 0) return 0;
1016: if (f->coeffp->tag != MP_INTEGER) return 0;
1017: if (DoCancel & 2) return 1;
1018: /* Apply the Noro strategy to reduce content. */
1019: size = mpz_size(f->coeffp->val.bigp);
1020: if (ss > 0) {
1021: prevSize = size;
1022: return 0;
1023: }
1024: if (size > 2*prevSize) {
1025: return 1;
1026: }else{
1027: return 0;
1028: }
1.4 takayama 1029: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to poly2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan