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