numbpart.cc

Go to the documentation of this file.

// compute number of partitions of n
// written by Thorsten Reinecke, 1999-12-06
// last change: 2005-06-05

// implementation is now threadsafe   

#include <gmp.h>
#include <stdexcept>
#include "utils.H"
#include <cmath>

namespace numbpart
{

inline int omega(const int m)
{
  return (3*m*m-m)/2;
} 

class Cnumbpart : private ForbidAssignment
{
 // this class can be utilized to provide access to a bunch of
 // (possible cached) partition numbers
private:
 const int size;
 typedef mpz_t* Pmpz_t;
 Pmpz_t* const p;
 void numbpart_recurse(const int n);

public:
 Cnumbpart(const int n) : size(n+1), p(new Pmpz_t [size])
  {
    for (int i=0; i<size; ++i) p[i]=NULL;
  }
 ~Cnumbpart()
  {
    for (int i=0; i<size; ++i) 
     if (p[i]) { mpz_clear(*p[i]); delete p[i]; p[i]=NULL; }
    delete [] p;
  }
 void get_numbpart(mpz_t res, const int n)
 {
   if (n<0) throw std::invalid_argument("negative argument in numbpart");
   if (n>=size) throw std::out_of_range("argument overflows provided range in numbpart");
   // if numbpart(n) has not yet been computed, then compute it...
   numbpart_recurse(n);
   mpz_set(res,*p[n]);
 } 
};

void Cnumbpart::numbpart_recurse(const int n)
{
  if (p[n]) return; // value is cached
  if (n==0) { p[n]=new mpz_t[1]; mpz_init_set_si(*p[n],1); return; }

  p[n]=new mpz_t[1];
  mpz_init(*p[n]);

#if 0
  // this one uses the stack massively!
  int m,i;
  m=1;
  while ((i=n-omega(m))>=0)
   {
     numbpart_recurse(i);
     if (m%2) mpz_add(*p[n],*p[n],*p[i]); else mpz_sub(*p[n],*p[n],*p[i]);
     m++;
   }
  m=1;
  while ((i=n-omega(-m))>=0)
   {
     numbpart_recurse(i);
     if (m%2) mpz_add(*p[n],*p[n],*p[i]); else mpz_sub(*p[n],*p[n],*p[i]);
     m++;
   }
#else
  // this consumes less stack space...
  int m = static_cast<int>(std::floor( -(1.0/6.0) + std::sqrt((2.0*n/3.0) + (1.0/36.0)) ));
  while (omega(-m)<=n) ++m;

  int i=n-omega(m);
  if (i>=0)
   {
     numbpart_recurse(i);
     if (m&1) mpz_add(*p[n],*p[n],*p[i]); else mpz_sub(*p[n],*p[n],*p[i]);
   }

  while (--m)
   {
     i=n-omega(m);
     numbpart_recurse(i);
     if (m&1) mpz_add(*p[n],*p[n],*p[i]); else mpz_sub(*p[n],*p[n],*p[i]);
     i-=m; // that is i=n-omega(-m)
     numbpart_recurse(i);
     if (m&1) mpz_add(*p[n],*p[n],*p[i]); else mpz_sub(*p[n],*p[n],*p[i]);
   }
#endif
}


void numbpart(mpz_t res, const int n)
{
  Cnumbpart capsule(n+1);
  capsule.get_numbpart(res,n);
}

}


#if 0
#include <iostream>
int main()
{
  using namespace std;
  mpz_t res;
  mpz_init(res);

  const unsigned int min = 0;
  const unsigned int max = 300000;
  const unsigned int step = 501;

#if 0
  cout << "Testing numbpart-function:" << endl;
  for (unsigned int i=min; i<=max; i+=step)
   {
     numbpart::numbpart(res,i);
     cout << "p(" << i << ") = ";
     mpz_out_str(stdout,10,res);
     cout << endl;
   }
#endif

  cout << "Testing numbpart cached access:" << endl;
  numbpart::Cnumbpart capsule(max);
  for (unsigned int i=min; i<=max; i+=step)
   {
     capsule.get_numbpart(res,i);
     cout << "p(" << i << ") = ";
     mpz_out_str(stdout,10,res);
     cout << endl;
   }

  mpz_clear(res);
}
#endif

---