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

Diff for /OpenXM/src/kan96xx/Kan/coeff.c between version 1.1 and 1.3

version 1.1, 1999/10/08 02:12:01 version 1.3, 2001/05/04 01:06:23
Line 1 
Line 1 
   /* $OpenXM: OpenXM/src/kan96xx/Kan/coeff.c,v 1.2 2000/01/16 07:55:38 takayama Exp $ */
 #include <stdio.h>  #include <stdio.h>
 #include "datatype.h"  #include "datatype.h"
 #include "stackm.h"  #include "stackm.h"
Line 17  static void printTag(coeffType t)
Line 18  static void printTag(coeffType t)
 }  }
   
 char *intToString(i)  char *intToString(i)
 int i;       int i;
 {  {
   char work[80];    char work[80];
   char *s;    char *s;
Line 29  int i;
Line 30  int i;
 }  }
   
 char *coeffToString(cp)  char *coeffToString(cp)
 struct coeff *cp;       struct coeff *cp;
 {  {
   char *s;    char *s;
   switch(cp->tag) {    switch(cp->tag) {
Line 53  struct coeff *cp;
Line 54  struct coeff *cp;
   
 /* constructors */  /* constructors */
 struct coeff *intToCoeff(i,ringp)  struct coeff *intToCoeff(i,ringp)
 int i;       int i;
 struct ring *ringp;       struct ring *ringp;
 {  {
   struct coeff *c;    struct coeff *c;
   int p;    int p;
Line 78  struct ring *ringp;
Line 79  struct ring *ringp;
 }  }
   
 struct coeff *mpintToCoeff(b,ringp)  struct coeff *mpintToCoeff(b,ringp)
 MP_INT *b;       MP_INT *b;
 struct ring *ringp;       struct ring *ringp;
 {  {
   struct coeff *c;    struct coeff *c;
   struct coeff *c2;    struct coeff *c2;
Line 102  struct ring *ringp;
Line 103  struct ring *ringp;
     c->val.f = coeffToPoly(c2,ringp->next);      c->val.f = coeffToPoly(c2,ringp->next);
     return(c);      return(c);
     /*      /*
     errorCoeff("mpintToCoeff(): mpint --> poly_coeff has not yet be supported. Returns 0.");        errorCoeff("mpintToCoeff(): mpint --> poly_coeff has not yet be supported. Returns 0.");
     c->tag = POLY_COEFF;        c->tag = POLY_COEFF;
     c->val.f = ZERO;        c->val.f = ZERO;
     */      */
   }    }
 }  }
   
 struct coeff *polyToCoeff(f,ringp)  struct coeff *polyToCoeff(f,ringp)
 POLY f;       POLY f;
 struct ring *ringp;       struct ring *ringp;
 {  {
   struct coeff *c;    struct coeff *c;
   if (f ISZERO) {    if (f ISZERO) {
Line 132  struct ring *ringp;
Line 133  struct ring *ringp;
   
 /* returns -c */  /* returns -c */
 struct coeff *coeffNeg(c,ringp)  struct coeff *coeffNeg(c,ringp)
 struct coeff *c;       struct coeff *c;
 struct ring *ringp;       struct ring *ringp;
 {  {
   struct coeff *r;    struct coeff *r;
   r = intToCoeff(-1,ringp);    r = intToCoeff(-1,ringp);
Line 142  struct ring *ringp;
Line 143  struct ring *ringp;
 }  }
   
 int coeffToInt(cp)  int coeffToInt(cp)
 struct coeff *cp;       struct coeff *cp;
 {  {
   POLY f;    POLY f;
   switch(cp->tag) {    switch(cp->tag) {
Line 163  struct coeff *cp;
Line 164  struct coeff *cp;
 }  }
   
 void Cadd(r,a,b)  void Cadd(r,a,b)
 struct coeff *r;       struct coeff *r;
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   int p;    int p;
   POLY f;    POLY f;
Line 174  struct coeff *b;
Line 175  struct coeff *b;
     r->p = p = a->p;      r->p = p = a->p;
     r->val.i = ((a->val.i) + (b->val.i)) % p;      r->val.i = ((a->val.i) + (b->val.i)) % p;
   }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER    }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
             && r->tag == MP_INTEGER) {              && r->tag == MP_INTEGER) {
     mpz_add(r->val.bigp,a->val.bigp,b->val.bigp);      mpz_add(r->val.bigp,a->val.bigp,b->val.bigp);
   }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF    }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF
             && r->tag == POLY_COEFF) {              && r->tag == POLY_COEFF) {
     f = ppAdd(a->val.f,b->val.f);      f = ppAdd(a->val.f,b->val.f);
     r->val.f = f;      r->val.f = f;
   }else {    }else {
Line 186  struct coeff *b;
Line 187  struct coeff *b;
 }  }
   
 void Csub(r,a,b)  void Csub(r,a,b)
 struct coeff *r;       struct coeff *r;
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   int p;    int p;
   
Line 196  struct coeff *b;
Line 197  struct coeff *b;
     r->p = p = a->p;      r->p = p = a->p;
     r->val.i = ((a->val.i) - (b->val.i)) % p;      r->val.i = ((a->val.i) - (b->val.i)) % p;
   }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER    }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
             && r->tag == MP_INTEGER) {              && r->tag == MP_INTEGER) {
     mpz_sub(r->val.bigp,a->val.bigp,b->val.bigp);      mpz_sub(r->val.bigp,a->val.bigp,b->val.bigp);
   }else {    }else {
     warningCoeff("Csub(): The data type is not supported.");      warningCoeff("Csub(): The data type is not supported.");
Line 204  struct coeff *b;
Line 205  struct coeff *b;
 }  }
   
 void Cmult(r,a,b)  void Cmult(r,a,b)
 struct coeff *r;       struct coeff *r;
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   int p;    int p;
   POLY f;    POLY f;
Line 215  struct coeff *b;
Line 216  struct coeff *b;
     r->p = p = a->p;      r->p = p = a->p;
     r->val.i = ((a->val.i) * (b->val.i)) % p;      r->val.i = ((a->val.i) * (b->val.i)) % p;
   }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER    }else if (a->tag == MP_INTEGER && b->tag == MP_INTEGER
             && r->tag == MP_INTEGER) {              && r->tag == MP_INTEGER) {
     mpz_mul(r->val.bigp,a->val.bigp,b->val.bigp);      mpz_mul(r->val.bigp,a->val.bigp,b->val.bigp);
   }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF    }else if (a->tag == POLY_COEFF && b->tag == POLY_COEFF
             && r->tag == POLY_COEFF) {              && r->tag == POLY_COEFF) {
     f = ppMult(a->val.f,b->val.f);      f = ppMult(a->val.f,b->val.f);
     r->val.f = f;      r->val.f = f;
   }else {    }else {
Line 231  struct coeff *b;
Line 232  struct coeff *b;
 }  }
   
 int isZero(a)  int isZero(a)
 struct coeff *a;       struct coeff *a;
 {  {
   int r;    int r;
   if (a->tag == INTEGER) {    if (a->tag == INTEGER) {
Line 251  struct coeff *a;
Line 252  struct coeff *a;
 }  }
   
 struct coeff *coeffCopy(c)  struct coeff *coeffCopy(c)
 struct coeff *c;       struct coeff *c;
 /* Deep copy */       /* Deep copy */
 {  {
   struct coeff *r;    struct coeff *r;
   r = newCoeff();    r = newCoeff();
Line 280  struct coeff *c;
Line 281  struct coeff *c;
 }  }
   
 void CiiComb(r,p,q)  void CiiComb(r,p,q)
 struct coeff *r;       struct coeff *r;
 int p,q;       int p,q;
 /* r->val.bigp is read only. */       /* r->val.bigp is read only. */
 {  {
   switch(r->tag) {    switch(r->tag) {
   case INTEGER:    case INTEGER:
Line 300  int p,q;
Line 301  int p,q;
   
   
 MP_INT *BiiComb(p,q)  MP_INT *BiiComb(p,q)
 int p,q;       int p,q;
 /* It returns {p \choose q} when p>=0, 0<=q<=p.       /* It returns {p \choose q} when p>=0, 0<=q<=p.
    The value is read only. DON'T CHANGE IT.          The value is read only. DON'T CHANGE IT.
 p=0         0                1          p=0         0                1
 p=1       1   2            1    1          p=1       1   2            1    1
 p=2     3   4   5        1   2    1          p=2     3   4   5        1   2    1
        q=0 q=1 q=2          q=0 q=1 q=2
        [p(p+1)/2+q]          [p(p+1)/2+q]
 */       */
 {  {
   extern MP_INT *Mp_one;    extern MP_INT *Mp_one;
   int i,j;    int i,j;
Line 326  p=2     3   4   5        1   2    1
Line 327  p=2     3   4   5        1   2    1
     if (newTable == (MP_INT **) NULL) errorCoeff("BiiComb(): No more memory.");      if (newTable == (MP_INT **) NULL) errorCoeff("BiiComb(): No more memory.");
     for (i=0; i<tableSize; i++) {      for (i=0; i<tableSize; i++) {
       for (j=0; j<=i; j++) {        for (j=0; j<=i; j++) {
         newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];          newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
       }        }
     }      }
     for (i=tableSize; i< 2*p; i++) {      for (i=tableSize; i< 2*p; i++) {
       for (j=0; j<=i; j++) {        for (j=0; j<=i; j++) {
         if (j==0 || j==i) newTable[(i*(i+1))/2 + j] = Mp_one;          if (j==0 || j==i) newTable[(i*(i+1))/2 + j] = Mp_one;
         else{          else{
           newTable[(i*(i+1))/2 + j] = newMP_INT();            newTable[(i*(i+1))/2 + j] = newMP_INT();
         }          }
       }        }
     }      }
     tableSize = 2*p;      tableSize = 2*p;
Line 344  p=2     3   4   5        1   2    1
Line 345  p=2     3   4   5        1   2    1
   for (i=maxp; i<=p; i++) {    for (i=maxp; i<=p; i++) {
     for (j=1; j<i; j++) {      for (j=1; j<i; j++) {
       mpz_add(table[(i*(i+1))/2 + j],        mpz_add(table[(i*(i+1))/2 + j],
              table[((i-1)*i)/2 + j-1],                table[((i-1)*i)/2 + j-1],
              table[((i-1)*i)/2 +j]);  /* [p-1,q-1]+[p-1,q] */                table[((i-1)*i)/2 +j]);  /* [p-1,q-1]+[p-1,q] */
     }      }
   }    }
   maxp = p+1;    maxp = p+1;
   return(table[(p*(p+1))/2 +q]);    return(table[(p*(p+1))/2 +q]);
            /*   ^^^^^^ No problem for the optimizer? */    /*   ^^^^^^ No problem for the optimizer? */
 }  }
   
 int iiComb(p,q,P)  int iiComb(p,q,P)
 int p,q,P;       int p,q,P;
 {  {
 /**    /**
  foo[n_]:=Block[{p,q,ans},       foo[n_]:=Block[{p,q,ans},
     ans={};       ans={};
     For[p=0,p<=n,p++,ans=Append[ans,Table[Binomial[p,q],{q,0,n}]]];       For[p=0,p<=n,p++,ans=Append[ans,Table[Binomial[p,q],{q,0,n}]]];
     Return[ans];       Return[ans];
  ]       ]
 **/    **/
 /* We assume that int is 32 bit */    /* We assume that int is 32 bit */
 static int table[26][26]=    static int table[26][26]=
 {{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},    {{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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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},     {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},
 {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,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}};
   
   int a,b;    int a,b;
   extern MP_INT Mp_work_iiComb;    extern MP_INT Mp_work_iiComb;
Line 410  static int table[26][26]=
Line 411  static int table[26][26]=
   
   
 MP_INT *BiiPoch(p,q)  MP_INT *BiiPoch(p,q)
 int p,q;       int p,q;
 {  {
 /* It returns  p(p-1) ... (p-q+1) = [p,q] when p>=0 and 0<=q or    /* It returns  p(p-1) ... (p-q+1) = [p,q] when p>=0 and 0<=q or
                                                p<0  and 0<=q.       p<0  and 0<=q.
    [p,q] = p*[p-1,q-1].  [p,0] = 1.       [p,q] = p*[p-1,q-1].  [p,0] = 1.
    The value is read only. DON'T CHANGE IT.       The value is read only. DON'T CHANGE IT.
   
 When p>=0,       When p>=0,
 p=0         0                1       p=0         0                1
 p=1       1   2            1    1       p=1       1   2            1    1
 p=2     3   4   5        1   2    2*1       p=2     3   4   5        1   2    2*1
 p=3   6   7   8   9    1   3   3*2    3*2*1,  if p<q,then 0.       p=3   6   7   8   9    1   3   3*2    3*2*1,  if p<q,then 0.
        [p(p+1)/2+q]       [p(p+1)/2+q]
   
 When p<0, [p,q] is indexed by (-p+q)(-p+q+1)/2 + q.       When p<0, [p,q] is indexed by (-p+q)(-p+q+1)/2 + q.
 -p+q = pq.       -p+q = pq.
   
 q=3  0  (-1)(-2)(-3)  (-2)(-3)(-4)   (-3)(-4)(-5)       q=3  0  (-1)(-2)(-3)  (-2)(-3)(-4)   (-3)(-4)(-5)
 q=2  0  (-1)(-2)      (-2)(-3)       (-3)(-4)       q=2  0  (-1)(-2)      (-2)(-3)       (-3)(-4)
 q=1  0   -1            -2             -3       q=1  0   -1            -2             -3
 q=0  1    1             1              1       q=0  1    1             1              1
    -p=0  -p=2          -p=3           -p=4       -p=0  -p=2          -p=3           -p=4
 */    */
   extern MP_INT *Mp_one;    extern MP_INT *Mp_one;
   extern MP_INT *Mp_zero;    extern MP_INT *Mp_zero;
   MP_INT mp_work;    MP_INT mp_work;
Line 460  q=0  1    1             1              1
Line 461  q=0  1    1             1              1
       newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*p*(2*p+1))/2));        newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*p*(2*p+1))/2));
       if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");        if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");
       for (i=0; i<tableSize; i++) {        for (i=0; i<tableSize; i++) {
         for (j=0; j<=i; j++) {          for (j=0; j<=i; j++) {
           newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];            newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
         }          }
       }        }
       for (i=tableSize; i< 2*p; i++) {        for (i=tableSize; i< 2*p; i++) {
         for (j=0; j<=i; j++) {          for (j=0; j<=i; j++) {
           if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;            if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;
           else{            else{
             newTable[(i*(i+1))/2 + j] = newMP_INT();              newTable[(i*(i+1))/2 + j] = newMP_INT();
           }            }
         }          }
       }        }
       tableSize = 2*p;        tableSize = 2*p;
       table = newTable;        table = newTable;
Line 478  q=0  1    1             1              1
Line 479  q=0  1    1             1              1
     /* Compute the binomial coefficients up to {p \choose p} */      /* Compute the binomial coefficients up to {p \choose p} */
     for (i=maxp; i<=p; i++) {      for (i=maxp; i<=p; i++) {
       for (j=1; j<=i; j++) {        for (j=1; j<=i; j++) {
         mpz_mul_ui(table[(i*(i+1))/2 + j],          mpz_mul_ui(table[(i*(i+1))/2 + j],
                    table[((i-1)*i)/2 + j-1],                     table[((i-1)*i)/2 + j-1],
                    (unsigned long int) i); /* [i,j] = i*[i-1,j-1] */                     (unsigned long int) i); /* [i,j] = i*[i-1,j-1] */
       }        }
     }      }
     maxp = p+1;      maxp = p+1;
Line 500  q=0  1    1             1              1
Line 501  q=0  1    1             1              1
       newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*pq*(2*pq+1))/2));        newTable = (MP_INT **)sGC_malloc(sizeof(MP_INT *)*( (2*pq*(2*pq+1))/2));
       if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");        if (newTable == (MP_INT **) NULL) errorCoeff("BiiPoch(): No more memory.");
       for (i=0; i<tablepqSize; i++) {        for (i=0; i<tablepqSize; i++) {
         for (j=0; j<=i; j++) {          for (j=0; j<=i; j++) {
           newTable[(i*(i+1))/2 + j] = tablepq[(i*(i+1))/2 + j];            newTable[(i*(i+1))/2 + j] = tablepq[(i*(i+1))/2 + j];
         }          }
       }        }
       for (i=tablepqSize; i< 2*pq; i++) {        for (i=tablepqSize; i< 2*pq; i++) {
         for (j=0; j<=i; j++) {          for (j=0; j<=i; j++) {
           if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;            if (j==0) newTable[(i*(i+1))/2 + j] = Mp_one;
           else if (j==i) newTable[(i*(i+1))/2+j] = Mp_zero;            else if (j==i) newTable[(i*(i+1))/2+j] = Mp_zero;
           else { /*^^^ no problem for optimizer? */            else { /*^^^ no problem for optimizer? */
             newTable[(i*(i+1))/2 + j] = newMP_INT();              newTable[(i*(i+1))/2 + j] = newMP_INT();
           }            }
         }          }
       }        }
       tablepqSize = 2*pq;        tablepqSize = 2*pq;
       tablepq = newTable;        tablepq = newTable;
Line 520  q=0  1    1             1              1
Line 521  q=0  1    1             1              1
     mpz_init(&mp_work);      mpz_init(&mp_work);
     for (i=maxpq; i<=pq; i++) {      for (i=maxpq; i<=pq; i++) {
       for (j=1; j<i; j++) {        for (j=1; j<i; j++) {
         mpz_set_si(&mp_work,-(i-j));          mpz_set_si(&mp_work,-(i-j));
         mpz_mul(tablepq[(i*(i+1))/2 + j],          mpz_mul(tablepq[(i*(i+1))/2 + j],
                 tablepq[(i*(i+1))/2 + j-1],                  tablepq[(i*(i+1))/2 + j-1],
                 &mp_work); /* [i,j] = i*[i-1,j-1] */                  &mp_work); /* [i,j] = i*[i-1,j-1] */
       }        }
     }      }
     maxpq = pq+1;      maxpq = pq+1;
Line 537  int iiPoch(p,k,P) int p,k,P; /* p(p-1)...(p-k+1) */ {
Line 538  int iiPoch(p,k,P) int p,k,P; /* p(p-1)...(p-k+1) */ {
   int r,i;    int r,i;
   
   /*    /*
   extern MP_INT Mp_work_iiPoch;      extern MP_INT Mp_work_iiPoch;
   mpz_mod_ui(&Mp_work_iiPoch,BiiPoch(p,k),(unsigned long) P);      mpz_mod_ui(&Mp_work_iiPoch,BiiPoch(p,k),(unsigned long) P);
   return((int) mpz_get_si(&Mp_work_iiPoch));      return((int) mpz_get_si(&Mp_work_iiPoch));
   100 test takes 16ms.      100 test takes 16ms.
   */    */
   
   if (p>=0 && p < k) return(0);    if (p>=0 && p < k) return(0);
Line 552  int iiPoch(p,k,P) int p,k,P; /* p(p-1)...(p-k+1) */ {
Line 553  int iiPoch(p,k,P) int p,k,P; /* p(p-1)...(p-k+1) */ {
 }  }
   
 void CiiPoch(r,p,q)  void CiiPoch(r,p,q)
 struct coeff *r;       struct coeff *r;
 int p,q;       int p,q;
 /* r->val.bigp is read only. */       /* r->val.bigp is read only. */
 {  {
   switch(r->tag) {    switch(r->tag) {
   case INTEGER:    case INTEGER:
Line 570  int p,q;
Line 571  int p,q;
 }  }
   
 MP_INT *BiiPower(p,q)  MP_INT *BiiPower(p,q)
 int p,q;       int p,q;
 /* It returns  p^q.  q must be >=0.       /* It returns  p^q.  q must be >=0.
    p^q = [p,q] = p*[p,q-1].          p^q = [p,q] = p*[p,q-1].
    The value is read only. DON'T CHANGE IT.          The value is read only. DON'T CHANGE IT.
 */       */
   /*       /*
 [p,q] is indexed by table2D[p+q,q]=table1D[(p+q)(p+q+1)/2 + q].         [p,q] is indexed by table2D[p+q,q]=table1D[(p+q)(p+q+1)/2 + q].
 p+q = pq.         p+q = pq.
   
 q=3  0    1             8             27         q=3  0    1             8             27
 q=2  0    1             4              9         q=2  0    1             4              9
 q=1  0    1             2              3         q=1  0    1             2              3
 q=0  1    1             1              1         q=0  1    1             1              1
    -p=0  -p=1          -p=2           -p=3         -p=0  -p=1          -p=2           -p=3
   */       */
 {  {
   extern MP_INT *Mp_one;    extern MP_INT *Mp_one;
   extern MP_INT *Mp_zero;    extern MP_INT *Mp_zero;
Line 621  q=0  1    1             1              1
Line 622  q=0  1    1             1              1
       errorCoeff("BiiPower(): No more memory.");        errorCoeff("BiiPower(): No more memory.");
     for (i=0; i<tableSize; i++) {      for (i=0; i<tableSize; i++) {
       for (j=0; j<=i; j++) {        for (j=0; j<=i; j++) {
         newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];          newTable[(i*(i+1))/2 + j] = table[(i*(i+1))/2 + j];
         newTable2[(i*(i+1))/2 + j] = tableneg[(i*(i+1))/2 + j];          newTable2[(i*(i+1))/2 + j] = tableneg[(i*(i+1))/2 + j];
       }        }
     }      }
     for (i=tableSize; i< 2*pq; i++) {      for (i=tableSize; i< 2*pq; i++) {
       for (j=0; j<=i; j++) {        for (j=0; j<=i; j++) {
         if (j==0) {          if (j==0) {
           newTable[(i*(i+1))/2 + j] = Mp_one;            newTable[(i*(i+1))/2 + j] = Mp_one;
           newTable2[(i*(i+1))/2 + j] = Mp_one;            newTable2[(i*(i+1))/2 + j] = Mp_one;
         } else if (j==i) {          } else if (j==i) {
           newTable[(i*(i+1))/2+j] = Mp_zero;            newTable[(i*(i+1))/2+j] = Mp_zero;
           newTable2[(i*(i+1))/2+j] = Mp_zero;            newTable2[(i*(i+1))/2+j] = Mp_zero;
         }else { /*^^^ no problem for optimizer? */          }else { /*^^^ no problem for optimizer? */
           newTable[(i*(i+1))/2 + j] = newMP_INT();            newTable[(i*(i+1))/2 + j] = newMP_INT();
           newTable2[(i*(i+1))/2 + j] = newMP_INT();            newTable2[(i*(i+1))/2 + j] = newMP_INT();
         }          }
       }        }
     }      }
     table = newTable;      table = newTable;
Line 649  q=0  1    1             1              1
Line 650  q=0  1    1             1              1
     for (j=1; j<i; j++) {      for (j=1; j<i; j++) {
       mpz_set_si(&mp_work,-(i-j));        mpz_set_si(&mp_work,-(i-j));
       mpz_mul(tableneg[(i*(i+1))/2 + j],        mpz_mul(tableneg[(i*(i+1))/2 + j],
               tableneg[((i-1)*i)/2 + j-1],                tableneg[((i-1)*i)/2 + j-1],
               &mp_work); /* (-(i-j))^j=[i,j] = (-i+j)*[i-1,j-1] */                &mp_work); /* (-(i-j))^j=[i,j] = (-i+j)*[i-1,j-1] */
       mpz_set_si(&mp_work,i-j);        mpz_set_si(&mp_work,i-j);
       mpz_mul(table[(i*(i+1))/2 + j],        mpz_mul(table[(i*(i+1))/2 + j],
               table[((i-1)*i)/2 + j-1],                table[((i-1)*i)/2 + j-1],
               &mp_work); /* (i-j)^j=[i,j] = (i-j)*[i-1,j-1] */                &mp_work); /* (i-j)^j=[i,j] = (i-j)*[i-1,j-1] */
     }      }
   }    }
   maxp = pq+1;    maxp = pq+1;
Line 666  q=0  1    1             1              1
Line 667  q=0  1    1             1              1
 }  }
   
 int iiPower(k,s,P)  int iiPower(k,s,P)
 int k;       int k;
 int s;  /* k^s ,  s>=0*/       int s;  /* k^s ,  s>=0*/
 int P;       int P;
 {  {
   int kk,r;    int kk,r;
   r = 1;    r = 1;
Line 679  int P;
Line 680  int P;
 }  }
   
 void CiiPower(r,p,q)  void CiiPower(r,p,q)
 struct coeff *r;       struct coeff *r;
 int p,q;       int p,q;
 /* r->val.bigp is read only. */       /* r->val.bigp is read only. */
 {  {
   switch(r->tag) {    switch(r->tag) {
   case INTEGER:    case INTEGER:
Line 697  int p,q;
Line 698  int p,q;
 }  }
   
 struct coeff *stringToUniversalNumber(s,flagp)  struct coeff *stringToUniversalNumber(s,flagp)
 char *s;       char *s;
 int *flagp;       int *flagp;
 {  {
   MP_INT *value;    MP_INT *value;
   struct coeff * c;    struct coeff * c;
Line 711  int *flagp;
Line 712  int *flagp;
 }  }
   
 struct coeff *newUniversalNumber(i)  struct coeff *newUniversalNumber(i)
 int i;       int i;
 { extern struct ring *SmallRingp;  { extern struct ring *SmallRingp;
   struct coeff *c;   struct coeff *c;
   c = intToCoeff(i,SmallRingp);   c = intToCoeff(i,SmallRingp);
   return(c);   return(c);
 }  }
   
 struct coeff *newUniversalNumber2(i)  struct coeff *newUniversalNumber2(i)
 MP_INT *i;       MP_INT *i;
 { extern struct ring *SmallRingp;  { extern struct ring *SmallRingp;
   struct coeff *c;   struct coeff *c;
   c = mpintToCoeff(i,SmallRingp);   c = mpintToCoeff(i,SmallRingp);
   return(c);   return(c);
 }  }
   
 int coeffEqual(a,b)  int coeffEqual(a,b)
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   if (a->tag == INTEGER && b->tag == INTEGER) {    if (a->tag == INTEGER && b->tag == INTEGER) {
     return((a->val.i) == (b->val.i));      return((a->val.i) == (b->val.i));
Line 743  struct coeff *b;
Line 744  struct coeff *b;
 }  }
   
 int coeffGreater(a,b)  int coeffGreater(a,b)
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   if (a->tag == INTEGER && b->tag == INTEGER) {    if (a->tag == INTEGER && b->tag == INTEGER) {
     return((a->val.i) - (b->val.i));      return((a->val.i) - (b->val.i));
Line 759  struct coeff *b;
Line 760  struct coeff *b;
 }  }
   
 void universalNumberDiv(r,a,b)  void universalNumberDiv(r,a,b)
 struct coeff *r;       struct coeff *r;
 struct coeff *a;       struct coeff *a;
 struct coeff *b;       struct coeff *b;
 {  {
   if (a->tag == MP_INTEGER && b->tag == MP_INTEGER && r->tag == MP_INTEGER) {    if (a->tag == MP_INTEGER && b->tag == MP_INTEGER && r->tag == MP_INTEGER) {
     mpz_div(r->val.bigp,a->val.bigp,b->val.bigp);      mpz_div(r->val.bigp,a->val.bigp,b->val.bigp);
Line 772  struct coeff *b;
Line 773  struct coeff *b;
   
 /* Note that it does not check if c is zero or not. cf isZero(). */  /* Note that it does not check if c is zero or not. cf isZero(). */
 POLY coeffToPoly(c,ringp)  POLY coeffToPoly(c,ringp)
 struct coeff *c;       struct coeff *c;
 struct ring *ringp;       struct ring *ringp;
 { int p;  { int p;
   struct coeff *tc;   struct coeff *tc;
   
   if (c->tag == INTEGER) {   if (c->tag == INTEGER) {
     tc = intToCoeff(c->val.i,ringp);     tc = intToCoeff(c->val.i,ringp);
     return(newCell(tc,newMonomial(ringp)));     return(newCell(tc,newMonomial(ringp)));
   }else if (c->tag == MP_INTEGER) {   }else if (c->tag == MP_INTEGER) {
     tc = mpintToCoeff(c->val.bigp,ringp);     tc = mpintToCoeff(c->val.bigp,ringp);
     return(newCell(tc,newMonomial(ringp)));     return(newCell(tc,newMonomial(ringp)));
   } else if (c->tag == POLY_COEFF) {   } else if (c->tag == POLY_COEFF) {
     return(newCell(c,newMonomial(ringp)));     return(newCell(c,newMonomial(ringp)));
   }else {   }else {
     warningCoeff("coeffToPoly(): The data type is not supported. Return 0.");     warningCoeff("coeffToPoly(): The data type is not supported. Return 0.");
     return(ZERO);     return(ZERO);
   }   }
 }  }
 void errorCoeff(str)  void errorCoeff(str)
 char *str;       char *str;
 {  {
   fprintf(stderr,"Error(coeff.c): %s\n",str);    fprintf(stderr,"Error(coeff.c): %s\n",str);
   fprintf(stderr,"Type in ctrl-\\");getchar();    fprintf(stderr,"Type in ctrl-\\");getchar();
Line 799  char *str;
Line 800  char *str;
 }  }
   
 void warningCoeff(str)  void warningCoeff(str)
 char *str;       char *str;
 {  {
   fprintf(stderr,"Warning(coeff.c): %s\n",str);    fprintf(stderr,"Warning(coeff.c): %s\n",str);
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>