[BACK]Return to coeff.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Kan

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>