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