Annotation of OpenXM/src/kan96xx/Kan/coeff.c, Revision 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