Annotation of OpenXM/src/kan96xx/Kan/coeff.c, Revision 1.4
1.4 ! ohara 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/coeff.c,v 1.3 2001/05/04 01:06:23 takayama Exp $ */
1.1 maekawa 2: #include <stdio.h>
1.4 ! ohara 3: #include <stdlib.h>
! 4: #include <string.h>
1.1 maekawa 5: #include "datatype.h"
6: #include "stackm.h"
7: #include "extern.h"
8: #include "extern2.h"
9:
10:
11: static void printTag(coeffType t)
12: {
13: switch(t) {
14: case UNKNOWN: fprintf(stderr," UNKNOWN "); break;
15: case INTEGER: fprintf(stderr," INTEGER "); break;
16: case MP_INTEGER: fprintf(stderr," MP_INTEGER "); break;
17: case POLY_COEFF: fprintf(stderr," POLY_COEFF "); break;
18: default: fprintf(stderr," Unknown tag of coeffType."); break;
19: }
20: }
21:
22: char *intToString(i)
1.3 takayama 23: int i;
1.1 maekawa 24: {
25: char work[80];
26: char *s;
27: sprintf(work,"%d",i);
28: s = (char *)sGC_malloc((strlen(work)+2)*sizeof(char));
29: if (s == (char *)NULL) errorCoeff("No more memory.");
30: strcpy(s,work);
31: return(s);
32: }
33:
34: char *coeffToString(cp)
1.3 takayama 35: struct coeff *cp;
1.1 maekawa 36: {
37: char *s;
38: switch(cp->tag) {
39: case INTEGER:
40: return(intToString((cp->val).i));
41: break;
42: case MP_INTEGER:
43: s = mpz_get_str((char *)NULL,10,(cp->val).bigp);
44: return(s);
45: break;
46: case POLY_COEFF:
47: s = POLYToString((cp->val).f,'*',1);
48: return(s);
49: break;
50: default:
51: warningCoeff("coeffToString: Unknown tag.");
52: return(" Unknown-coeff ");
53: break;
54: }
55: }
56:
57: /* constructors */
58: struct coeff *intToCoeff(i,ringp)
1.3 takayama 59: int i;
60: struct ring *ringp;
1.1 maekawa 61: {
62: struct coeff *c;
63: int p;
64: c =newCoeff();
65: if (ringp->next != (struct ring *)NULL) {
66: c->tag = POLY_COEFF;
67: c->val.f = cxx(i,0,0,ringp->next);
68: c->p = ringp->p;
69: }else{
70: p = ringp->p;
71: if (p) {
72: c->tag = INTEGER; c->p = p;
73: c->val.i = i % p;
74: }else{
75: c->tag = MP_INTEGER; c->p = 0;
76: c->val.bigp = newMP_INT();
77: mpz_set_si(c->val.bigp,(long) i);
78: }
79: }
80: return(c);
81: }
82:
83: struct coeff *mpintToCoeff(b,ringp)
1.3 takayama 84: MP_INT *b;
85: struct ring *ringp;
1.1 maekawa 86: {
87: struct coeff *c;
88: struct coeff *c2;
89: int p;
90: c =newCoeff();
91: p = ringp->p;
92: if (ringp->next == (struct ring *)NULL) {
93: if (p) {
94: c->tag = INTEGER; c->p = p;
95: c->val.i = ((int) mpz_get_si(b)) % p;
96: }else{
97: c->tag = MP_INTEGER; c->p = 0;
98: c->val.bigp = b;
99: }
100: return(c);
101: }else{
102: c2 = mpintToCoeff(b,ringp->next);
103: c->tag = POLY_COEFF;
104: c->p = ringp->p;
105: c->val.f = coeffToPoly(c2,ringp->next);
106: return(c);
107: /*
1.3 takayama 108: errorCoeff("mpintToCoeff(): mpint --> poly_coeff has not yet be supported. Returns 0.");
109: c->tag = POLY_COEFF;
110: c->val.f = ZERO;
1.1 maekawa 111: */
112: }
113: }
114:
115: struct coeff *polyToCoeff(f,ringp)
1.3 takayama 116: POLY f;
117: struct ring *ringp;
1.1 maekawa 118: {
119: struct coeff *c;
120: if (f ISZERO) {
121: c =newCoeff();
122: c->tag = POLY_COEFF; c->p = ringp->p;
123: c->val.f = f;
124: return(c);
125: }
126: if (f->m->ringp != ringp->next) {
127: errorCoeff("polyToCoeff(): f->m->ringp != ringp->next. Returns 0.");
128: f = 0;
129: }
130: c =newCoeff();
131: c->tag = POLY_COEFF; c->p = ringp->p;
132: c->val.f = f;
133: return(c);
134: }
135:
136: /* returns -c */
137: struct coeff *coeffNeg(c,ringp)
1.3 takayama 138: struct coeff *c;
139: struct ring *ringp;
1.1 maekawa 140: {
141: struct coeff *r;
142: r = intToCoeff(-1,ringp);
143: Cmult(r,r,c);
144: return(r);
145: }
146:
147: int coeffToInt(cp)
1.3 takayama 148: struct coeff *cp;
1.1 maekawa 149: {
150: POLY f;
151: switch(cp->tag) {
152: case INTEGER:
153: return(cp->val.i);
154: break;
155: case MP_INTEGER:
156: return(mpz_get_si(cp->val.bigp));
157: break;
158: case POLY_COEFF:
159: f = cp->val.f;
160: if (f == ZERO) return(0);
161: else return(coeffToInt(f->coeffp));
162: break;
163: default:
164: errorCoeff("Unknown tag in coeffToInt().\n");
165: }
166: }
167:
168: void Cadd(r,a,b)
1.3 takayama 169: struct coeff *r;
170: struct coeff *a;
171: struct coeff *b;
1.1 maekawa 172: {
173: int p;
174: POLY f;
175:
176: if (a->tag == INTEGER && b->tag == INTEGER && r->tag == INTEGER) {
177: r->p = p = a->p;
178: r->val.i = ((a->val.i) + (b->val.i)) % p;
179: }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
1.3 takayama 180: && r->tag == MP_INTEGER) {
1.1 maekawa 181: mpz_add(r->val.bigp,a->val.bigp,b->val.bigp);
182: }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF
1.3 takayama 183: && r->tag == POLY_COEFF) {
1.1 maekawa 184: f = ppAdd(a->val.f,b->val.f);
185: r->val.f = f;
186: }else {
187: warningCoeff("Cadd(): The data type is not supported.");
188: }
189: }
190:
191: void Csub(r,a,b)
1.3 takayama 192: struct coeff *r;
193: struct coeff *a;
194: struct coeff *b;
1.1 maekawa 195: {
196: int p;
197:
198: if (a->tag == INTEGER && b->tag == INTEGER && r->tag == INTEGER) {
199: r->p = p = a->p;
200: r->val.i = ((a->val.i) - (b->val.i)) % p;
201: }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
1.3 takayama 202: && r->tag == MP_INTEGER) {
1.1 maekawa 203: mpz_sub(r->val.bigp,a->val.bigp,b->val.bigp);
204: }else {
205: warningCoeff("Csub(): The data type is not supported.");
206: }
207: }
208:
209: void Cmult(r,a,b)
1.3 takayama 210: struct coeff *r;
211: struct coeff *a;
212: struct coeff *b;
1.1 maekawa 213: {
214: int p;
215: POLY f;
216:
217: if (a->tag == INTEGER && b->tag == INTEGER && r->tag == INTEGER) {
218: r->p = p = a->p;
219: r->val.i = ((a->val.i) * (b->val.i)) % p;
220: }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
1.3 takayama 221: && r->tag == MP_INTEGER) {
1.1 maekawa 222: mpz_mul(r->val.bigp,a->val.bigp,b->val.bigp);
223: }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF
1.3 takayama 224: && r->tag == POLY_COEFF) {
1.1 maekawa 225: f = ppMult(a->val.f,b->val.f);
226: r->val.f = f;
227: }else {
228: warningCoeff("Cmult(): The data type is not supported.");
229: printTag(a->tag); printTag(b->tag); printTag(r->tag); fprintf(stderr,"\n");
230: warningCoeff("Returns coeffCopy(a) ------------------");
231: printf("Type in return.\n");getchar(); getchar();
232: *r = *(coeffCopy(a));
233: }
234: }
235:
236: int isZero(a)
1.3 takayama 237: struct coeff *a;
1.1 maekawa 238: {
239: int r;
240: if (a->tag == INTEGER) {
241: if (a->val.i) return(0);
242: else return(1);
243: }else if (a->tag == MP_INTEGER) {
244: r = mpz_cmp_si(a->val.bigp,(long)0);
245: if (r == 0) return(1);
246: else return(0);
247: }else if (a->tag == POLY_COEFF) {
248: if (a->val.f ISZERO) return(1);
249: else return(0);
250: }else{
251: warningCoeff("CisZero(): The data type is not supported.");
252: return(1);
253: }
254: }
255:
256: struct coeff *coeffCopy(c)
1.3 takayama 257: struct coeff *c;
258: /* Deep copy */
1.1 maekawa 259: {
260: struct coeff *r;
261: r = newCoeff();
262: switch(c->tag) {
263: case INTEGER:
264: r->tag = INTEGER;
265: r->p = c->p;
266: r->val.i = c->val.i;
267: break;
268: case MP_INTEGER:
269: r->tag = MP_INTEGER; r->p = 0;
270: r->val.bigp = newMP_INT();
271: mpz_set(r->val.bigp,c->val.bigp);
272: break;
273: case POLY_COEFF:
274: r->tag = POLY_COEFF;
275: r->p = c->p;
276: r->val.f = pcmCopy(c->val.f); /* The polynomial is deeply copied. */
277: break;
278: default:
279: warningCoeff("coeffCopy(): Unknown tag for the argument.");
280: break;
281: }
282: return(r);
283: }
284:
285: void CiiComb(r,p,q)
1.3 takayama 286: struct coeff *r;
287: int p,q;
288: /* r->val.bigp is read only. */
1.1 maekawa 289: {
290: switch(r->tag) {
291: case INTEGER:
292: r->val.i = iiComb(p,q,r->p);
293: break;
294: case MP_INTEGER:
295: r->val.bigp = BiiComb(p,q);
296: break;
297: default:
298: warningCoeff("CiiComb(): Unknown tag.");
299: break;
300: }
301: }
302:
303:
304:
305: MP_INT *BiiComb(p,q)
1.3 takayama 306: int p,q;
307: /* It returns {p \choose q} when p>=0, 0<=q<=p.
308: The value is read only. DON'T CHANGE IT.
309: p=0 0 1
310: p=1 1 2 1 1
311: p=2 3 4 5 1 2 1
312: q=0 q=1 q=2
313: [p(p+1)/2+q]
314: */
1.1 maekawa 315: {
316: extern MP_INT *Mp_one;
317: int i,j;
318: MP_INT **newTable;
319: static MP_INT **table;
320: static int tableSize = 0;
321: static int maxp = 0; /* You have {p \choose q} in the table when p<maxp. */
322: if (p<=0 || q<=0) return(Mp_one);
323: if (p<=q) return(Mp_one);
324: if (p<maxp) return(table[(p*(p+1))/2+q]);
325: /* Enlarge the table if it's small. */
326: if ( !(p < tableSize) ) {
327: /* The new size is 2*p. */
328: newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*p*(2*p+1))/2));
329: if (newTable == (MP_INT **) NULL) errorCoeff("BiiComb(): No more memory.");
330: for (i=0; i<tableSize; i++) {
331: for (j=0; j<=i; j++) {
1.3 takayama 332: newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
1.1 maekawa 333: }
334: }
335: for (i=tableSize; i< 2*p; i++) {
336: for (j=0; j<=i; j++) {
1.3 takayama 337: if (j==0 || j==i) newTable[(i*(i+1))/2 + j] = Mp_one;
338: else{
339: newTable[(i*(i+1))/2 + j] = newMP_INT();
340: }
1.1 maekawa 341: }
342: }
343: tableSize = 2*p;
344: table = newTable;
345: }
346: /* Compute the binomial coefficients up to {p \choose p} */
347: for (i=maxp; i<=p; i++) {
348: for (j=1; j<i; j++) {
349: mpz_add(table[(i*(i+1))/2 + j],
1.3 takayama 350: table[((i-1)*i)/2 + j-1],
351: table[((i-1)*i)/2 +j]); /* [p-1,q-1]+[p-1,q] */
1.1 maekawa 352: }
353: }
354: maxp = p+1;
355: return(table[(p*(p+1))/2 +q]);
1.3 takayama 356: /* ^^^^^^ No problem for the optimizer? */
1.1 maekawa 357: }
358:
359: int iiComb(p,q,P)
1.3 takayama 360: int p,q,P;
1.1 maekawa 361: {
1.3 takayama 362: /**
363: foo[n_]:=Block[{p,q,ans},
364: ans={};
365: For[p=0,p<=n,p++,ans=Append[ans,Table[Binomial[p,q],{q,0,n}]]];
366: Return[ans];
367: ]
368: **/
369: /* We assume that int is 32 bit */
370: static int table[26][26]=
371: {{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
372: {1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
373: {1,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
374: {1,3,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
375: {1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
376: {1,5,10,10,5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
377: {1,6,15,20,15,6,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
378: {1,7,21,35,35,21,7,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
379: {1,8,28,56,70,56,28,8,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
380: {1,9,36,84,126,126,84,36,9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
381: {1,10,45,120,210,252,210,120,45,10,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
382: {1,11,55,165,330,462,462,330,165,55,11,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
383: {1,12,66,220,495,792,924,792,495,220,66,12,1,0,0,0,0,0,0,0,0,0,0,0,0,0},
384: {1,13,78,286,715,1287,1716,1716,1287,715,286,78,13,1,0,0,0,0,0,0,0,0,0,0,0,0},
385: {1,14,91,364,1001,2002,3003,3432,3003,2002,1001,364,91,14,1,0,0,0,0,0,0,0,0,0,0,0},
386: {1,15,105,455,1365,3003,5005,6435,6435,5005,3003,1365,455,105,15,1,0,0,0,0,0,0,0,0,0,0},
387: {1,16,120,560,1820,4368,8008,11440,12870,11440,8008,4368,1820,560,120,16,1,0,0,0,0,0,0,0,0,0},
388: {1,17,136,680,2380,6188,12376,19448,24310,24310,19448,12376,6188,2380,680,136,17,1,0,0,0,0,0,0,0,0},
389: {1,18,153,816,3060,8568,18564,31824,43758,48620,43758,31824,18564,8568,3060,816,153,18,1,0,0,0,0,0,0,0},
390: {1,19,171,969,3876,11628,27132,50388,75582,92378,92378,75582,50388,27132,11628,3876,969,171,19,1,0,0,0,0,0,0},
391: {1,20,190,1140,4845,15504,38760,77520,125970,167960,184756,167960,125970,77520,38760,15504,4845,1140,190,20,1,0,0,0,0,0},
392: {1,21,210,1330,5985,20349,54264,116280,203490,293930,352716,352716,293930,203490,116280,54264,20349,5985,1330,210,21,1,0,0,0,0},
393: {1,22,231,1540,7315,26334,74613,170544,319770,497420,646646,705432,646646,497420,319770,170544,74613,26334,7315,1540,231,22,1,0,0,0},
394: {1,23,253,1771,8855,33649,100947,245157,490314,817190,1144066,1352078,1352078,1144066,817190,490314,245157,100947,33649,8855,1771,253,23,1,0,0},
395: {1,24,276,2024,10626,42504,134596,346104,735471,1307504,1961256,2496144,2704156,2496144,1961256,1307504,735471,346104,134596,42504,10626,2024,276,24,1,0},
396: {1,25,300,2300,12650,53130,177100,480700,1081575,2042975,3268760,4457400,5200300,5200300,4457400,3268760,2042975,1081575,480700,177100,53130,12650,2300,300,25,1}};
1.1 maekawa 397:
398: int a,b;
399: extern MP_INT Mp_work_iiComb;
400:
401: if ((p<=0) || (q<=0)) return( 1 );
402: if (p <= q) return( 1 );
403: if (p<26) return( table[p][q] % P);
404: else {
405: mpz_mod_ui(&Mp_work_iiComb,BiiComb(p,q),(unsigned long) P);
406: return((int) mpz_get_si(&Mp_work_iiComb));
407: }
408: /*a=iiComb(p-1,q-1,P); b=iiComb(p-1,q,P);
409: a = a+b; a %= P;
410: return(a);
411: */
412: }
413:
414:
415: MP_INT *BiiPoch(p,q)
1.3 takayama 416: int p,q;
1.1 maekawa 417: {
1.3 takayama 418: /* It returns p(p-1) ... (p-q+1) = [p,q] when p>=0 and 0<=q or
419: p<0 and 0<=q.
420: [p,q] = p*[p-1,q-1]. [p,0] = 1.
421: The value is read only. DON'T CHANGE IT.
422:
423: When p>=0,
424: p=0 0 1
425: p=1 1 2 1 1
426: p=2 3 4 5 1 2 2*1
427: p=3 6 7 8 9 1 3 3*2 3*2*1, if p<q,then 0.
428: [p(p+1)/2+q]
429:
430: When p<0, [p,q] is indexed by (-p+q)(-p+q+1)/2 + q.
431: -p+q = pq.
432:
433: q=3 0 (-1)(-2)(-3) (-2)(-3)(-4) (-3)(-4)(-5)
434: q=2 0 (-1)(-2) (-2)(-3) (-3)(-4)
435: q=1 0 -1 -2 -3
436: q=0 1 1 1 1
437: -p=0 -p=2 -p=3 -p=4
438: */
1.1 maekawa 439: extern MP_INT *Mp_one;
440: extern MP_INT *Mp_zero;
441: MP_INT mp_work;
442: int i,j;
443: MP_INT **newTable;
444: static MP_INT **table;
445: static int tableSize = 0;
446: static int maxp = 0;
447: static MP_INT **tablepq;
448: static int tablepqSize = 0;
449: static int maxpq = 0;
450: int pq;
451:
452: if (p>=0) {
453: if (q < 0) {
454: warningCoeff("Negative argument to BiiPoch(). Returns zero.");
455: return(Mp_zero);
456: }
457: if (q == 0) return(Mp_one);
458: if (p<q) return(Mp_zero);
459: if (p<maxp) return(table[(p*(p+1))/2+q]);
460: /* Enlarge the table if it's small. */
461: if ( !(p < tableSize) ) {
462: /* The new size is 2*p. */
463: newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*p*(2*p+1))/2));
464: if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");
465: for (i=0; i<tableSize; i++) {
1.3 takayama 466: for (j=0; j<=i; j++) {
467: newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
468: }
1.1 maekawa 469: }
470: for (i=tableSize; i< 2*p; i++) {
1.3 takayama 471: for (j=0; j<=i; j++) {
472: if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;
473: else{
474: newTable[(i*(i+1))/2 + j] = newMP_INT();
475: }
476: }
1.1 maekawa 477: }
478: tableSize = 2*p;
479: table = newTable;
480: }
481: /* Compute the binomial coefficients up to {p \choose p} */
482: for (i=maxp; i<=p; i++) {
483: for (j=1; j<=i; j++) {
1.3 takayama 484: mpz_mul_ui(table[(i*(i+1))/2 + j],
485: table[((i-1)*i)/2 + j-1],
486: (unsigned long int) i); /* [i,j] = i*[i-1,j-1] */
1.1 maekawa 487: }
488: }
489: maxp = p+1;
490: return(table[(p*(p+1))/2 +q]);
491: /* ^^^^^^ No problem for the optimizer? */
492: }else{
493: if (q < 0) {
494: warningCoeff("Negative argument to BiiPoch(). Returns zero.");
495: return(Mp_zero);
496: }
497: if (q == 0) return(Mp_one);
498: pq = -p+q;
499: if (pq<maxpq) return(tablepq[(pq*(pq+1))/2+q]);
500: /* Enlarge the table if it's small. */
501: if ( !(pq < tablepqSize) ) {
502: /* The new size is 2*pq. */
503: newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*pq*(2*pq+1))/2));
504: if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");
505: for (i=0; i<tablepqSize; i++) {
1.3 takayama 506: for (j=0; j<=i; j++) {
507: newTable[(i*(i+1))/2 + j] = tablepq[(i*(i+1))/2 + j];
508: }
1.1 maekawa 509: }
510: for (i=tablepqSize; i< 2*pq; i++) {
1.3 takayama 511: for (j=0; j<=i; j++) {
512: if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;
513: else if (j==i) newTable[(i*(i+1))/2+j] = Mp_zero;
514: else { /*^^^ no problem for optimizer? */
515: newTable[(i*(i+1))/2 + j] = newMP_INT();
516: }
517: }
1.1 maekawa 518: }
519: tablepqSize = 2*pq;
520: tablepq = newTable;
521: }
522: /* Compute the Poch up to [p,q], -p+q<=pq */
523: mpz_init(&mp_work);
524: for (i=maxpq; i<=pq; i++) {
525: for (j=1; j<i; j++) {
1.3 takayama 526: mpz_set_si(&mp_work,-(i-j));
527: mpz_mul(tablepq[(i*(i+1))/2 + j],
528: tablepq[(i*(i+1))/2 + j-1],
529: &mp_work); /* [i,j] = i*[i-1,j-1] */
1.1 maekawa 530: }
531: }
532: maxpq = pq+1;
533: return(tablepq[(pq*(pq+1))/2 +q]);
534: /* ^^^^^^ No problem for the optimizer? */
535: }
536: }
537:
538:
539: int iiPoch(p,k,P) int p,k,P; /* p(p-1)...(p-k+1) */ {
540: int r,i;
541:
542: /*
1.3 takayama 543: extern MP_INT Mp_work_iiPoch;
544: mpz_mod_ui(&Mp_work_iiPoch,BiiPoch(p,k),(unsigned long) P);
545: return((int) mpz_get_si(&Mp_work_iiPoch));
546: 100 test takes 16ms.
1.1 maekawa 547: */
548:
549: if (p>=0 && p < k) return(0);
550: r = 1;
551: for (i=0;i<k;i++) {
552: r = r*(p-i); r %= P;
553: }
554: return(r);
555: }
556:
557: void CiiPoch(r,p,q)
1.3 takayama 558: struct coeff *r;
559: int p,q;
560: /* r->val.bigp is read only. */
1.1 maekawa 561: {
562: switch(r->tag) {
563: case INTEGER:
564: r->val.i = iiPoch(p,q,r->p);
565: break;
566: case MP_INTEGER:
567: r->val.bigp = BiiPoch(p,q);
568: break;
569: default:
570: warningCoeff("CiiPoch(): Unknown tag.");
571: break;
572: }
573: }
574:
575: MP_INT *BiiPower(p,q)
1.3 takayama 576: int p,q;
577: /* It returns p^q. q must be >=0.
578: p^q = [p,q] = p*[p,q-1].
579: The value is read only. DON'T CHANGE IT.
580: */
581: /*
582: [p,q] is indexed by table2D[p+q,q]=table1D[(p+q)(p+q+1)/2 + q].
583: p+q = pq.
584:
585: q=3 0 1 8 27
586: q=2 0 1 4 9
587: q=1 0 1 2 3
588: q=0 1 1 1 1
589: -p=0 -p=1 -p=2 -p=3
590: */
1.1 maekawa 591: {
592: extern MP_INT *Mp_one;
593: extern MP_INT *Mp_zero;
594: MP_INT mp_work;
595: int i,j;
596: MP_INT **newTable;
597: MP_INT **newTable2;
598: static MP_INT **table;
599: static int tableSize = 0;
600: static int maxp = 0;
601: static MP_INT **tableneg;
602: int pq;
603:
604:
605: if (q < 0) {
606: warningCoeff("Negative argument to BiiPower(). Returns zero.");
607: return(Mp_zero);
608: }
609: if (q == 0) return(Mp_one);
610: if (p>=0) {
611: pq = p+q;
612: }else {
613: pq = -p+q;
614: }
615: if (pq<maxp && p>=0) return(table[(pq*(pq+1))/2+q]);
616: if (pq<maxp && p<0) return(tableneg[(pq*(pq+1))/2+q]);
617:
618: /* Enlarge the table if it's small. */
619: if ( !(pq < tableSize) ) {
620: /* The new size is 2*pq. */
621: newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*pq*(2*pq+1))/2));
622: newTable2 = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*pq*(2*pq+1))/2));
623: if (newTable == (MP_INT **) NULL || newTable2 == (MP_INT **)NULL)
624: errorCoeff("BiiPower(): No more memory.");
625: for (i=0; i<tableSize; i++) {
626: for (j=0; j<=i; j++) {
1.3 takayama 627: newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
628: newTable2[(i*(i+1))/2 + j] = tableneg[(i*(i+1))/2 + j];
1.1 maekawa 629: }
630: }
631: for (i=tableSize; i< 2*pq; i++) {
632: for (j=0; j<=i; j++) {
1.3 takayama 633: if (j==0) {
634: newTable[(i*(i+1))/2 + j] = Mp_one;
635: newTable2[(i*(i+1))/2 + j] = Mp_one;
636: } else if (j==i) {
637: newTable[(i*(i+1))/2+j] = Mp_zero;
638: newTable2[(i*(i+1))/2+j] = Mp_zero;
639: }else { /*^^^ no problem for optimizer? */
640: newTable[(i*(i+1))/2 + j] = newMP_INT();
641: newTable2[(i*(i+1))/2 + j] = newMP_INT();
642: }
1.1 maekawa 643: }
644: }
645: table = newTable;
646: tableneg = newTable2;
647: tableSize = 2*pq;
648: }
649: /* Compute the Power up to [p,q], p+q<=pq */
650: mpz_init(&mp_work);
651: for (i=maxp; i<=pq; i++) {
652: for (j=1; j<i; j++) {
653: mpz_set_si(&mp_work,-(i-j));
654: mpz_mul(tableneg[(i*(i+1))/2 + j],
1.3 takayama 655: tableneg[((i-1)*i)/2 + j-1],
656: &mp_work); /* (-(i-j))^j=[i,j] = (-i+j)*[i-1,j-1] */
1.1 maekawa 657: mpz_set_si(&mp_work,i-j);
658: mpz_mul(table[(i*(i+1))/2 + j],
1.3 takayama 659: table[((i-1)*i)/2 + j-1],
660: &mp_work); /* (i-j)^j=[i,j] = (i-j)*[i-1,j-1] */
1.1 maekawa 661: }
662: }
663: maxp = pq+1;
664: if (p>=0) {
665: return(table[(pq*(pq+1))/2 +q]);
666: }else{
667: return(tableneg[(pq*(pq+1))/2 +q]);
668: }
669: }
670:
671: int iiPower(k,s,P)
1.3 takayama 672: int k;
673: int s; /* k^s , s>=0*/
674: int P;
1.1 maekawa 675: {
676: int kk,r;
677: r = 1;
678: for (kk=0; kk<s; kk++ ) {
679: r *= k; r %= P;
680: }
681: return(r);
682: }
683:
684: void CiiPower(r,p,q)
1.3 takayama 685: struct coeff *r;
686: int p,q;
687: /* r->val.bigp is read only. */
1.1 maekawa 688: {
689: switch(r->tag) {
690: case INTEGER:
691: r->val.i = iiPower(p,q,r->p);
692: break;
693: case MP_INTEGER:
694: r->val.bigp = BiiPower(p,q);
695: break;
696: default:
697: warningCoeff("CiiPower(): Unknown tag.");
698: break;
699: }
700: }
701:
702: struct coeff *stringToUniversalNumber(s,flagp)
1.3 takayama 703: char *s;
704: int *flagp;
1.1 maekawa 705: {
706: MP_INT *value;
707: struct coeff * c;
708: value =newMP_INT();
709: *flagp = mpz_set_str(value,s,10);
710: c = newCoeff();
711: c->tag = MP_INTEGER; c->p = 0;
712: c->val.bigp = value;
713: return(c);
714: }
715:
716: struct coeff *newUniversalNumber(i)
1.3 takayama 717: int i;
1.1 maekawa 718: { extern struct ring *SmallRingp;
1.3 takayama 719: struct coeff *c;
720: c = intToCoeff(i,SmallRingp);
721: return(c);
1.1 maekawa 722: }
723:
724: struct coeff *newUniversalNumber2(i)
1.3 takayama 725: MP_INT *i;
1.1 maekawa 726: { extern struct ring *SmallRingp;
1.3 takayama 727: struct coeff *c;
728: c = mpintToCoeff(i,SmallRingp);
729: return(c);
1.1 maekawa 730: }
731:
732: int coeffEqual(a,b)
1.3 takayama 733: struct coeff *a;
734: struct coeff *b;
1.1 maekawa 735: {
736: if (a->tag == INTEGER && b->tag == INTEGER) {
737: return((a->val.i) == (b->val.i));
738: }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER) {
739: if (mpz_cmp(a->val.bigp,b->val.bigp)==0) return(1);
740: else return(0);
741: }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF) {
742: return(ppSub(a->val.f,b->val.f) ISZERO);
743: }else {
744: warningCoeff("coeffEqual(): The data type is not supported.");
745: }
746: }
747:
748: int coeffGreater(a,b)
1.3 takayama 749: struct coeff *a;
750: struct coeff *b;
1.1 maekawa 751: {
752: if (a->tag == INTEGER && b->tag == INTEGER) {
753: return((a->val.i) - (b->val.i));
754: }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER) {
755: return(mpz_cmp(a->val.bigp,b->val.bigp));
756: }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF) {
757: warningCoeff("coeffGreater(POLY,POLY) always returns 0.");
758: return(0);
759: }else {
760: warningCoeff("coeffGreater(): The data type is not supported.");
761: }
762: }
763:
764: void universalNumberDiv(r,a,b)
1.3 takayama 765: struct coeff *r;
766: struct coeff *a;
767: struct coeff *b;
1.1 maekawa 768: {
769: if (a->tag == MP_INTEGER && b->tag == MP_INTEGER && r->tag == MP_INTEGER) {
770: mpz_div(r->val.bigp,a->val.bigp,b->val.bigp);
771: }else {
772: warningCoeff("universalNumberDiv(): The data type is not supported.");
773: }
774: }
775:
776: /* Note that it does not check if c is zero or not. cf isZero(). */
777: POLY coeffToPoly(c,ringp)
1.3 takayama 778: struct coeff *c;
779: struct ring *ringp;
1.1 maekawa 780: { int p;
1.3 takayama 781: struct coeff *tc;
1.1 maekawa 782:
1.3 takayama 783: if (c->tag == INTEGER) {
784: tc = intToCoeff(c->val.i,ringp);
785: return(newCell(tc,newMonomial(ringp)));
786: }else if (c->tag == MP_INTEGER) {
787: tc = mpintToCoeff(c->val.bigp,ringp);
788: return(newCell(tc,newMonomial(ringp)));
789: } else if (c->tag == POLY_COEFF) {
790: return(newCell(c,newMonomial(ringp)));
791: }else {
792: warningCoeff("coeffToPoly(): The data type is not supported. Return 0.");
793: return(ZERO);
794: }
1.1 maekawa 795: }
796: void errorCoeff(str)
1.3 takayama 797: char *str;
1.1 maekawa 798: {
799: fprintf(stderr,"Error(coeff.c): %s\n",str);
800: fprintf(stderr,"Type in ctrl-\\");getchar();
801: exit(15);
802: }
803:
804: void warningCoeff(str)
1.3 takayama 805: char *str;
1.1 maekawa 806: {
807: fprintf(stderr,"Warning(coeff.c): %s\n",str);
808: }
809:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to coeff.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan