binomial.c

   1 #include "binomial.h"
   2
   3 struct binom {
   4     Vector      **binom;
   5     unsigned    size;
   6     unsigned    n;
   7 };
   8 static struct binom binom;
   9
  10 Value *binomial(unsigned n, unsigned k)
  11 {
  12     int i, j;
  13
  14     if (n < binom.n)
  15         return &binom.binom[n]->p[k];
  16
  17     if (n >= binom.size) {
  18         int size = 3*(n + 5)/2;
  19
  20         binom.binom = (Vector **)realloc(binom.binom, size*sizeof(Vector *));
  21         binom.size = size;
  22     }
  23     for (i = binom.n; i <= n; ++i) {
  24         binom.binom[i] = Vector_Alloc(i+1);
  25         if (!i)
  26             value_set_si(binom.binom[0]->p[0], 1);
  27         else {
  28             value_set_si(binom.binom[i]->p[0], 1);
  29             value_set_si(binom.binom[i]->p[i], 1);
  30             for (j = 1; j < i; ++j)
  31                 value_addto(binom.binom[i]->p[j],
  32                             binom.binom[i-1]->p[j-1], binom.binom[i-1]->p[j]);
  33         }
  34     }
  35     binom.n = n+1;
  36     return &binom.binom[n]->p[k];
  37 }
  38
  39 struct fact {
  40     Value       *fact;
  41     unsigned    size;
  42     unsigned    n;
  43 };
  44 static struct fact fact;
  45
  46 Value *factorial(unsigned n)
  47 {
  48     int i;
  49
  50     if (n < fact.n)
  51         return &fact.fact[n];
  52
  53     if (n >= fact.size) {
  54         int size = 3*(n + 5)/2;
  55
  56         fact.fact = (Value *)realloc(fact.fact, size*sizeof(Value));
  57         fact.size = size;
  58     }
  59     for (i = fact.n; i <= n; ++i) {
  60         value_init(fact.fact[i]);
  61         if (!i)
  62             value_set_si(fact.fact[0], 1);
  63         else
  64             mpz_mul_ui(fact.fact[i], fact.fact[i-1], i);
  65     }
  66     fact.n = n+1;
  67     return &fact.fact[n];
  68 }