Index: bdm/estim/arx.h
===================================================================
--- bdm/estim/arx.h (revision 135)
+++ bdm/estim/arx.h (revision 145)
@@ -51,5 +51,7 @@
 	//! Full constructor
 	ARX (RV &rv, mat &V0, double &nu0, double frg0=1.0) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();tll=0.0;};
+	//! Set sufficient statistics
 	void set_parameters(mat &V0, double &nu0){est._V()=V0;est._nu()=nu0;last_lognc=est.lognc();tll=last_lognc;}
+	//! Returns sufficient statistics
 	void get_parameters(mat &V0, double &nu0){V0=est._V().to_mat(); nu0=est._nu();}
 	//! Here \f$dt = [y_t psi_t] \f$.
Index: bdm/itpp_ext.h
===================================================================
--- bdm/itpp_ext.h (revision 91)
+++ bdm/itpp_ext.h (revision 145)
@@ -13,5 +13,5 @@
 
 #ifndef ITEX_H
-#define ITEX_H 
+#define ITEX_H
 
 using std::cout;
@@ -19,45 +19,55 @@
 
 namespace itpp {
-Array<int> to_Arr ( const ivec &indices );
-ivec linspace ( int from, int to ); 
+	Array<int> to_Arr ( const ivec &indices );
+	ivec linspace ( int from, int to );
 
-/*!
-  \brief Gamma distribution
-  \ingroup randgen
-*/
-class Gamma_RNG {
-public:
+	bvec operator& ( const bvec &a, const bvec &b );
+	bvec operator| ( const bvec &a, const bvec &b );
+
+// template<class Num_T>
+// void set_subvector(vec &ov, ivec &iv, const Vec<Num_T> &v);
+
+	void set_subvector ( vec &ov, const ivec &iv, const vec &v );
+
+	/*!
+	  \brief Gamma distribution
+	  \ingroup randgen
+	*/
+	class Gamma_RNG {
+		public:
 //! constructor. Set lambda.
-Gamma_RNG ( double a=1.0, double b=1.0 );
-	//! Set lambda
-	void setup ( double a0, double b0 ) { alpha=a0; beta=b0;}
-	//! get lambda
-	double get_setup() const;
-	//! Get one sample.
-	double operator() () { return sample(); }
-	//! Get a sample vector.
-	vec operator() ( int n );
-	//! Get a sample matrix.
-	mat operator() ( int h, int w );
-protected:
-private:
-	//!
-	double sample();
-	//!
-	double alpha;
-	//!
-	double beta;
-	//!
-	Random_Generator RNG;
-	Normal_RNG NRNG;
-	//! compatibility with R
-	inline double exp_rand() {return -std::log ( RNG.random_01() );}
-	inline double unif_rand() {return RNG.random_01();}
-	inline double norm_rand() {return NRNG.sample();}
+			Gamma_RNG ( double a=1.0, double b=1.0 );
+			//! Set lambda
+			void setup ( double a0, double b0 ) { alpha=a0; beta=b0;}
+			//! get lambda
+			double get_setup() const;
+			//! Get one sample.
+			double operator() () { return sample(); }
+			//! Get a sample vector.
+			vec operator() ( int n );
+			//! Get a sample matrix.
+			mat operator() ( int h, int w );
+		protected:
+		private:
+			//!
+			double sample();
+			//!
+			double alpha;
+			//!
+			double beta;
+			//!
+			Random_Generator RNG;
+			Normal_RNG NRNG;
+			//! compatibility with R
+			inline double exp_rand() {return -std::log ( RNG.random_01() );}
+			inline double unif_rand() {return RNG.random_01();}
+			inline double norm_rand() {return NRNG.sample();}
 
-};
+	};
 
-bool qr ( const mat &A, mat &R );
+	bool qr ( const mat &A, mat &R );
 
 }
+
+
 #endif //ITEX_H
Index: bdm/stat/emix.cpp
===================================================================
--- bdm/stat/emix.cpp (revision 141)
+++ bdm/stat/emix.cpp (revision 145)
@@ -18,5 +18,5 @@
 
 	int i=0;
-	while ( ( w ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+	while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
 
 	return Coms ( i )->sample();
Index: bdm/stat/libBM.cpp
===================================================================
--- bdm/stat/libBM.cpp (revision 102)
+++ bdm/stat/libBM.cpp (revision 145)
@@ -25,5 +25,5 @@
 	times = in_times;
 	tsize = 0;
-	for ( i=0;i<len;i++ ) {tsize+=sizes ( i );}
+	for ( i = 0;i < len;i++ ) {tsize += sizes ( i );}
 };
 
@@ -32,14 +32,27 @@
 }
 
-RV::RV () : tsize ( 0 ),len ( 0 ) {};
+RV::RV() : tsize ( 0 ), len ( 0 ) {};
 
-void RV::add ( const RV &rv2 ) {
+bool RV::add ( const RV &rv2 ) {
 	// TODO
-	tsize+=rv2.tsize;
-	len +=rv2.len;
-	ids=concat ( ids,rv2.ids );
-	sizes=concat ( sizes,rv2.sizes );
-	times=concat ( times,rv2.times );
-	names=concat ( names,rv2.names );
+	ivec ind = rv2.findself ( *this ); //should be -1 all the time
+	ivec index = itpp::find(ind==-1); 
+	
+	if ( index.length() < rv2.len ) { //conflict
+		ids = concat ( ids, rv2.ids(index) );
+		sizes = concat ( sizes, rv2.sizes(index) );
+		times = concat ( times, rv2.times(index) );
+		names = concat ( names, rv2.names(to_Arr(index)) );
+	}
+	else {
+		ids = concat ( ids, rv2.ids );
+		sizes = concat ( sizes, rv2.sizes );
+		times = concat ( times, rv2.times );
+		names = concat ( names, rv2.names );
+	}
+	tsize = sum(sizes);
+	len = ids.length();
+	return (index.length()<rv2.len);
+
 //	return *this;
 };
@@ -58,21 +71,21 @@
 }
 
-RV RV::subselect ( ivec ind ) {
+RV RV::subselect ( ivec ind ) const {
 	return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
 }
 
-void RV::t ( int delta ) { times +=delta;}
+void RV::t ( int delta ) { times += delta;}
 
-RV RV::operator() ( ivec ind ) {
+RV RV::operator() ( ivec ind ) const {
 	return RV ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
 }
 
-bool RV::equal ( RV rv2 ) const {
-	return ( ids==rv2.ids ) && ( times == rv2.times ) && ( sizes==rv2.sizes );
+bool RV::equal ( const RV &rv2 ) const {
+	return ( ids == rv2.ids ) && ( times == rv2.times ) && ( sizes == rv2.sizes );
 }
 
 mat epdf::sampleN ( int N ) const {
-	mat X =zeros ( rv.count(),N );
-	for ( int i=0;i<N;i++ ) X.set_col ( i,this->sample() ); 
+	mat X = zeros ( rv.count(), N );
+	for ( int i = 0;i < N;i++ ) X.set_col ( i, this->sample() );
 	return X;
 };
@@ -88,12 +101,47 @@
 }
 
-ivec RV::indexlist() {
-	ivec indlist ( tsize );
+str RV::tostr() const {
+	ivec idlist ( tsize );
+	ivec tmlist ( tsize );
 	int i;
 	int pos = 0;
-	for ( i=0;i<len;i++ ) {
-		indlist.set_subvector ( pos,pos+sizes ( i )-1, ids ( i ) );
+	for ( i = 0;i < len;i++ ) {
+		idlist.set_subvector ( pos, pos + sizes ( i ) - 1, ids ( i ) );
+		tmlist.set_subvector ( pos, pos + sizes ( i ) - 1, times ( i ) );
+		pos += sizes ( i );
 	}
-	return indlist;
+	return str ( idlist, tmlist );
+}
+
+ivec RV::dataind ( RV rv2 ) const {
+	str str2 = rv2.tostr();
+	ivec res ( 0 );
+	ivec part;
+	int i;
+	for ( i = 0;i < len;i++ ) {
+		part = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
+		res = concat ( res, part );
+	}
+	return res;
+}
+
+RV RV::subt ( const RV rv2 ) const {
+	ivec res = this->findself ( rv2 ); // nonzeros
+	ivec valid = itpp::find ( res == -1 ); //-1 => value not found => it remains
+	return ( *this ) ( valid ); //keep those that were not found in rv2
+}
+
+ivec RV::findself ( const RV &rv2 ) const {
+	int i, j;
+	ivec tmp = -ones_i ( len );
+	for ( i = 0;i < len;i++ ) {
+		for ( j = 0;j < rv2.length();j++ ) {
+			if ( ( ids ( i ) == rv2.ids ( j ) ) & ( times ( i ) == rv2.times ( j ) ) ) {
+				tmp ( i ) = j;
+				break;
+			}
+		}
+	}
+	return tmp;
 }
 
Index: bdm/stat/emix.h
===================================================================
--- bdm/stat/emix.h (revision 141)
+++ bdm/stat/emix.h (revision 145)
@@ -30,6 +30,5 @@
 
 */
-class emix : public epdf
-{
+class emix : public epdf {
 	protected:
 		//! weights of the components
@@ -39,21 +38,19 @@
 	public:
 		//!Default constructor
-		emix ( RV &rv ) : epdf ( rv ) {};
+		emix(RV &rv) : epdf(rv) {};
 		//! Set weights \c w and components \c R
-		void set_parameters ( const vec &w, const Array<epdf*> &Coms );
+		void set_parameters(const vec &w, const Array<epdf*> &Coms);
 
 		vec sample() const;
-		vec mean() const
-		{
-			int i; vec mu=zeros ( rv.count() );
-			for ( i=0;i<w.length();i++ ) {mu+=w ( i ) *Coms ( i )->mean(); }
+		vec mean() const {
+			int i; vec mu = zeros(rv.count());
+			for (i = 0;i < w.length();i++) {mu += w(i) * Coms(i)->mean(); }
 			return mu;
 		}
-		double evalpdflog ( const vec &val ) const
-		{
+		double evalpdflog(const vec &val) const {
 			int i;
-			double sum=0.0;
-			for ( i=0;i<w.length();i++ ) {sum+=w ( i ) *Coms ( i )->evalpdflog ( val );}
-			return log ( sum );
+			double sum = 0.0;
+			for (i = 0;i < w.length();i++) {sum += w(i) * Coms(i)->evalpdflog(val);}
+			return log(sum);
 		};
 
@@ -65,21 +62,103 @@
 /*! \brief Chain rule decomposition of epdf
 
+Probability density in the form of Chain-rule decomposition:
+\[
+f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
+\]
+Note that
+*/
+class eprod: public epdf {
+	protected:
+		//
+		int n;
+		// pointers to epdfs
+		Array<epdf*> epdfs;
+		Array<mpdf*> mpdfs;
+		//
+		Array<ivec> rvinds;
+		Array<ivec> rvcinds;
+		//! Indicate independence of its factors
+		bool independent;
+		//! Indicate internal creation of mpdfs which must be destroyed
+		bool intermpdfs;
+	public:
+		//!Constructor from list of eFacs or list of mFacs
+		eprod(Array<mpdf*> mFacs): epdf(RV()), n(mFacs.length()), epdfs(n), mpdfs(mFacs), rvinds(n), rvcinds(n) {
+			int i;
+			intermpdfs = false;
+			for (i = 0;i < n;i++) {
+				rv.add(mpdfs(i)->_rv()); //add rv to common rvs.
+				epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
+			};
+			independent = true;
+			//test rvc of mpdfs and fill rvinds
+			for (i = 0;i < n;i++) {
+				// find ith rv in common rv
+				rvinds(i) = mpdfs(i)->_rv().dataind(rv);
+				// find ith rvc in common rv
+				rvcinds(i) = mpdfs(i)->_rvc().dataind(rv);
+				if (rvcinds(i).length()>0) {independent = false;}
+			}
+
+		};
+		eprod(Array<epdf*> eFacs): epdf(RV()), n(eFacs.length()), epdfs(eFacs), mpdfs(n), rvinds(n), rvcinds(n) {
+			int i;
+			intermpdfs = true;
+			for (i = 0;i < n;i++) {
+				if (rv.add(eFacs(i)->_rv())) {it_error("Incompatible eFacs.rv!");} //add rv to common rvs.
+				mpdfs(i) = new mepdf(*(epdfs(i)));
+				epdfs(i) = &(mpdfs(i)->_epdf()); // add pointer to epdf
+			};
+			independent = true;
+			//test rvc of mpdfs and fill rvinds
+			for (i = 0;i < n;i++) {
+				// find ith rv in common rv
+				rvinds(i) = mpdfs(i)->_rv().dataind(rv);
+			}
+		};
+
+		double evalpdflog(const vec &val) const {
+			int i;
+			double res = 0.0;
+			for (i = n - 1;i > 0;i++) {
+				if (rvcinds(i).length() > 0)
+					{mpdfs(i)->condition(val(rvcinds(i)));}
+				// add logarithms
+				res += epdfs(i)->evalpdflog(val(rvinds(i)));
+			}
+		}
+		vec sample () const{
+			vec smp=zeros(rv.count());
+			for (int i = (n - 1);i >= 0;i--) {
+				if (rvcinds(i).length() > 0)
+					{mpdfs(i)->condition(smp(rvcinds(i)));}
+				set_subvector(smp,rvinds(i), epdfs(i)->sample());
+			}			
+			return smp;
+		}
+		vec mean() const{
+			vec tmp(rv.count());
+			if (independent) {
+				for (int i=0;i<n;i++) {
+					vec pom = epdfs(i)->mean();
+					set_subvector(tmp,rvinds(i), pom);
+				}
+			}
+			else {
+				int N=50*rv.count();
+				it_warning("eprod.mean() computed by sampling");
+				tmp = zeros(rv.count());
+				for (int i=0;i<N;i++){ tmp += sample();}
+				tmp /=N;
+			}
+			return tmp;
+		}
+		~eprod(){if (intermpdfs) {for (int i=0;i<n;i++){delete mpdfs(i);}}};
+};
+
+/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type
 
 */
-class eprod: public epdf
-{
-	protected:
-		Array<epdf*> epdfs;
-		Array<mpdf*> mpdfs;
-	public:
-
-
-};
-
-/*! \brief Mixture of mpdfs with constant weights
-
-*/
-class mmix : public mpdf
-{
+class mmix : public mpdf {
 	protected:
 		//! Component (epdfs)
@@ -89,20 +168,17 @@
 	public:
 		//!Default constructor
-		mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep=&Epdf;};
+		mmix(RV &rv, RV &rvc) : mpdf(rv, rvc), Epdf(rv) {ep = &Epdf;};
 		//! Set weights \c w and components \c R
-		void set_parameters ( const vec &w, const Array<mpdf*> &Coms )
-		{
-			Array<epdf*> Eps ( Coms.length() );
+		void set_parameters(const vec &w, const Array<mpdf*> &Coms) {
+			Array<epdf*> Eps(Coms.length());
 
-			for ( int i=0;i<Coms.length();i++ )
-			{
-				Eps ( i ) =& ( Coms ( i )->_epdf() );
+			for (int i = 0;i < Coms.length();i++) {
+				Eps(i) = & (Coms(i)->_epdf());
 			}
-			Epdf.set_parameters ( w,Eps );
+			Epdf.set_parameters(w, Eps);
 		};
 
-		void condition ( const vec &cond )
-		{
-			for ( int i=0;i<Coms.length();i++ ) {Coms ( i )->condition ( cond );}
+		void condition(const vec &cond) {
+			for (int i = 0;i < Coms.length();i++) {Coms(i)->condition(cond);}
 		};
 };
Index: bdm/stat/loggers.h
===================================================================
--- bdm/stat/loggers.h (revision 102)
+++ bdm/stat/loggers.h (revision 145)
@@ -84,4 +84,5 @@
 	void step(bool final=false) {if ( ind<maxlen ) ind++; else it_error ( "memlog::ind is too high;" );}
 	void logit ( int id, vec v ) {vectors ( id ).set_row ( ind,v );}
+	//! Save values into an itfile named after \c fname.
 	void itsave(const char* fname);
 };
Index: bdm/stat/libBM.h
===================================================================
--- bdm/stat/libBM.h (revision 118)
+++ bdm/stat/libBM.h (revision 145)
@@ -18,4 +18,14 @@
 
 using namespace itpp;
+
+//! Structure of RV (used internally)
+class str{
+public:
+	ivec ids;
+	ivec times;
+	str(ivec ids0, ivec times0):ids(ids0),times(times0){
+		it_assert_debug(times0.length()==ids0.length(),"Incompatible input");
+	};
+};
 
 /*!
@@ -61,22 +71,22 @@
 	//TODO why not inline and later??
 
-	//! Find indexes of another rv in self
-	ivec find ( RV rv2 );
+	//! Find indexes of self in another rv, \return ivec of the same size as self.
+	ivec findself (const RV &rv2 ) const;
 	//! Compare if \c rv2 is identical to this \c RV
-	bool equal ( RV rv2 ) const;
-	//! Add (concat) another variable to the current one
-	void add ( const RV &rv2 );
-	//! Add (concat) another variable to the current one
-	friend RV concat ( const RV &rv1, const RV &rv2 );
+	bool equal (const RV &rv2 ) const;
+	//! Add (concat) another variable to the current one, \return 0 if all rv2 were added, 1 if rv2 is in conflict
+	bool add ( const RV &rv2 );
 	//! Subtract  another variable from the current one
-	RV subt ( RV rv2 );
+	RV subt ( const RV rv2 ) const;
 	//! Select only variables at indeces ind
-	RV subselect ( ivec ind );
+	RV subselect ( ivec ind ) const;
 	//! Select only variables at indeces ind
-	RV operator() ( ivec ind );
-	//! Generate new \c RV with \c time shifted by delta.
+	RV operator() ( ivec ind ) const;
+	//! Shift \c time shifted by delta.
 	void t ( int delta );
-	//! generate a list of indeces, i.e. which
-	ivec indexlist();
+	//! generate \c str from rv, by expanding sizes
+	str tostr() const;
+	//! generate indeces into \param crv data vector that form data vector of self.
+	ivec dataind(RV crv) const;
 
 	//!access function
@@ -92,4 +102,7 @@
 	std::string name ( int at ) {return names ( at );};
 };
+
+//! Concat two random variables
+RV concat ( const RV &rv1, const RV &rv2 );
 
 
@@ -146,6 +159,6 @@
 	//! Destructor for future use;
 	virtual ~epdf() {};
-	//! access function
-	RV _rv() const {return rv;}
+	//! access function, possibly dangerous!
+	RV& _rv() {return rv;}
 };
 
@@ -170,5 +183,5 @@
 	vec temp= ep->sample();
 	ll=ep->evalpdflog ( temp );return temp;};
-	//! Returns N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
+	//! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
 	virtual mat samplecond ( vec &cond, vec &ll, int N ) {
 		this->condition ( cond );
@@ -190,4 +203,6 @@
 	//! access function
 	RV _rvc() {return rvc;}
+	//! access function
+	RV _rv() {return rv;}
 	//!access function
 	epdf& _epdf() {return *ep;}
@@ -201,5 +216,5 @@
 public:
 	//!Default constructor
-	mepdf ( const RV &rv, const RV &rvc, epdf* em ) :mpdf ( rv,rvc ) {ep=em;};
+	mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
 };
 
Index: bdm/itpp_ext.cpp
===================================================================
--- bdm/itpp_ext.cpp (revision 141)
+++ bdm/itpp_ext.cpp (revision 145)
@@ -2,5 +2,5 @@
 // C++ Implementation: itpp_ext
 //
-// Description: 
+// Description:
 //
 //
@@ -16,33 +16,64 @@
 
 extern "C" {  /* QR factorization of a general matrix A  */
-  void dgeqrf_(int *m, int *n, double *a, int *lda, double *tau, double *work,
-	       int *lwork, int *info);
+	void dgeqrf_ ( int *m, int *n, double *a, int *lda, double *tau, double *work,
+	               int *lwork, int *info );
 };
 
 namespace itpp {
-	Array<int> to_Arr(const ivec &indices)
-    {
-    Array<int> a(indices.size());
-    for (int i = 0; i < a.size(); i++) {
-      a(i) = indices(i);
-    }
-    return a;
-  }
-
-ivec linspace ( int from, int to ) {
-	int n=to-from+1; 
-	int i;
-	it_assert_debug ( n>0,"wrong linspace" ); 
-	ivec iv ( n ); for (i=0;i<n;i++) iv(i)=from+i;
-	return iv;
-};
+	Array<int> to_Arr ( const ivec &indices ) {
+		Array<int> a ( indices.size() );
+		for ( int i = 0; i < a.size(); i++ ) {
+			a ( i ) = indices ( i );
+		}
+		return a;
+	}
+
+	ivec linspace ( int from, int to ) {
+		int n=to-from+1;
+		int i;
+		it_assert_debug ( n>0,"wrong linspace" );
+		ivec iv ( n ); for ( i=0;i<n;i++ ) iv ( i ) =from+i;
+		return iv;
+	};
+
+	void set_subvector ( vec &ov, const ivec &iv, const vec &v )
+	{
+		it_assert_debug ( ( iv.length() <=v.length() ),
+							"Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
+									"of range of v" );
+		for ( int i = 0; i < iv.length(); i++ ) {
+			it_assert_debug ( iv ( i ) <ov.length(),
+							  "Vec<>::set_subvector(ivec, vec<Num_T>): Indexing out "
+									  "of range of v" );
+			ov ( iv ( i ) ) = v ( i );
+		}
+	}
 
 //Gamma
 
-  Gamma_RNG::Gamma_RNG(double a, double b)
-  {
-    setup(a,b);
-  }
-
+	Gamma_RNG::Gamma_RNG ( double a, double b ) {
+		setup ( a,b );
+	}
+
+
+	bvec operator& ( const bvec &a, const bvec &b ) {
+		it_assert_debug ( b.size() ==a.size(), "operator&(): Vectors of different lengths" );
+
+		bvec temp ( a.size() );
+		for ( int i = 0;i < a.size();i++ ) {
+			temp ( i ) = a ( i ) & b ( i );
+		}
+		return temp;
+	}
+
+	bvec operator| ( const bvec &a, const bvec &b ) {
+		it_assert_debug ( b.size() !=a.size(), "operator&(): Vectors of different lengths" );
+
+		bvec temp ( a.size() );
+		for ( int i = 0;i < a.size();i++ ) {
+			temp ( i ) = a ( i ) | b ( i );
+		}
+		return temp;
+	}
 #define log std::log
 #define exp std::exp
@@ -50,199 +81,200 @@
 #define R_FINITE std::isfinite
 
-double  Gamma_RNG::sample()
-  {
-  //A copy of rgamma code from the R package!!
-  // 
-  
-/* Constants : */
-    const static double sqrt32 = 5.656854;
-    const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
-
-    /* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
-     * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k)
-     * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
-     */
-    const static double q1 = 0.04166669;
-    const static double q2 = 0.02083148;
-    const static double q3 = 0.00801191;
-    const static double q4 = 0.00144121;
-    const static double q5 = -7.388e-5;
-    const static double q6 = 2.4511e-4;
-    const static double q7 = 2.424e-4;
-
-    const static double a1 = 0.3333333;
-    const static double a2 = -0.250003;
-    const static double a3 = 0.2000062;
-    const static double a4 = -0.1662921;
-    const static double a5 = 0.1423657;
-    const static double a6 = -0.1367177;
-    const static double a7 = 0.1233795;
-
-    /* State variables [FIXME for threading!] :*/
-    static double aa = 0.;
-    static double aaa = 0.;
-    static double s, s2, d;    /* no. 1 (step 1) */
-    static double q0, b, si, c;/* no. 2 (step 4) */
-
-    double e, p, q, r, t, u, v, w, x, ret_val;
-    double a=alpha;
-    double scale=1.0/beta;
-
-    if (!R_FINITE(a) || !R_FINITE(scale) || a < 0.0 || scale <= 0.0)
-	{it_error("Gamma_RNG wrong parameters");}
-
-    if (a < 1.) { /* GS algorithm for parameters a < 1 */
-	if(a == 0)
-	    return 0.;
-	e = 1.0 + exp_m1 * a;
-	for(;;) {  //VS repeat
-	    p = e * unif_rand();
-	    if (p >= 1.0) {
-		x = -log((e - p) / a);
-		if (exp_rand() >= (1.0 - a) * log(x))
-		    break;
-	    } else {
-		x = exp(log(p) / a);
-		if (exp_rand() >= x)
-		    break;
-	    }
-	}
-	return scale * x;
-    }
-
-    /* --- a >= 1 : GD algorithm --- */
-
-    /* Step 1: Recalculations of s2, s, d if a has changed */
-    if (a != aa) {
-	aa = a;
-	s2 = a - 0.5;
-	s = sqrt(s2);
-	d = sqrt32 - s * 12.0;
-    }
-    /* Step 2: t = standard normal deviate,
-               x = (s,1/2) -normal deviate. */
-
-    /* immediate acceptance (i) */
-    t = norm_rand();
-    x = s + 0.5 * t;
-    ret_val = x * x;
-    if (t >= 0.0)
-	return scale * ret_val;
-
-    /* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
-    u = unif_rand();
-    if ((d * u) <= (t * t * t))
-	return scale * ret_val;
-
-    /* Step 4: recalculations of q0, b, si, c if necessary */
-
-    if (a != aaa) {
-	aaa = a;
-	r = 1.0 / a;
-	q0 = ((((((q7 * r + q6) * r + q5) * r + q4) * r + q3) * r
-	       + q2) * r + q1) * r;
-
-	/* Approximation depending on size of parameter a */
-	/* The constants in the expressions for b, si and c */
-	/* were established by numerical experiments */
-
-	if (a <= 3.686) {
-	    b = 0.463 + s + 0.178 * s2;
-	    si = 1.235;
-	    c = 0.195 / s - 0.079 + 0.16 * s;
-	} else if (a <= 13.022) {
-	    b = 1.654 + 0.0076 * s2;
-	    si = 1.68 / s + 0.275;
-	    c = 0.062 / s + 0.024;
-	} else {
-	    b = 1.77;
-	    si = 0.75;
-	    c = 0.1515 / s;
-	}
-    }
-    /* Step 5: no quotient test if x not positive */
-
-    if (x > 0.0) {
-	/* Step 6: calculation of v and quotient q */
-	v = t / (s + s);
-	if (fabs(v) <= 0.25)
-	    q = q0 + 0.5 * t * t * ((((((a7 * v + a6) * v + a5) * v + a4) * v
-				      + a3) * v + a2) * v + a1) * v;
-	else
-	    q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
-
-
-	/* Step 7: quotient acceptance (q) */
-	if (log(1.0 - u) <= q)
-	    return scale * ret_val;
-    }
-
-    for(;;) { //VS repeat
-	/* Step 8: e = standard exponential deviate
-	 *	u =  0,1 -uniform deviate
-	 *	t = (b,si)-double exponential (laplace) sample */
-	e = exp_rand();
-	u = unif_rand();
-	u = u + u - 1.0;
-	if (u < 0.0)
-	    t = b - si * e;
-	else
-	    t = b + si * e;
-	/* Step	 9:  rejection if t < tau(1) = -0.71874483771719 */
-	if (t >= -0.71874483771719) {
-	    /* Step 10:	 calculation of v and quotient q */
-	    v = t / (s + s);
-	    if (fabs(v) <= 0.25)
-		q = q0 + 0.5 * t * t *
-		    ((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v
-		      + a2) * v + a1) * v;
-	    else
-		q = q0 - s * t + 0.25 * t * t + (s2 + s2) * log(1.0 + v);
-	    /* Step 11:	 hat acceptance (h) */
-	    /* (if q not positive go to step 8) */
-	    if (q > 0.0) {
-		// TODO: w = expm1(q);
-		w = exp(q)-1;
-		/*  ^^^^^ original code had approximation with rel.err < 2e-7 */
-		/* if t is rejected sample again at step 8 */
-		if ((c * fabs(u)) <= (w * exp(e - 0.5 * t * t)))
-		    break;
-	    }
-	}
-    } /* repeat .. until  `t' is accepted */
-    x = s + 0.5 * t;
-    return scale * x * x;
+	double  Gamma_RNG::sample() {
+		//A copy of rgamma code from the R package!!
+		//
+
+		/* Constants : */
+		const static double sqrt32 = 5.656854;
+		const static double exp_m1 = 0.36787944117144232159;/* exp(-1) = 1/e */
+
+		/* Coefficients q[k] - for q0 = sum(q[k]*a^(-k))
+		 * Coefficients a[k] - for q = q0+(t*t/2)*sum(a[k]*v^k)
+		 * Coefficients e[k] - for exp(q)-1 = sum(e[k]*q^k)
+		 */
+		const static double q1 = 0.04166669;
+		const static double q2 = 0.02083148;
+		const static double q3 = 0.00801191;
+		const static double q4 = 0.00144121;
+		const static double q5 = -7.388e-5;
+		const static double q6 = 2.4511e-4;
+		const static double q7 = 2.424e-4;
+
+		const static double a1 = 0.3333333;
+		const static double a2 = -0.250003;
+		const static double a3 = 0.2000062;
+		const static double a4 = -0.1662921;
+		const static double a5 = 0.1423657;
+		const static double a6 = -0.1367177;
+		const static double a7 = 0.1233795;
+
+		/* State variables [FIXME for threading!] :*/
+		static double aa = 0.;
+		static double aaa = 0.;
+		static double s, s2, d;    /* no. 1 (step 1) */
+		static double q0, b, si, c;/* no. 2 (step 4) */
+
+		double e, p, q, r, t, u, v, w, x, ret_val;
+		double a=alpha;
+		double scale=1.0/beta;
+
+		if ( !R_FINITE ( a ) || !R_FINITE ( scale ) || a < 0.0 || scale <= 0.0 )
+			{it_error ( "Gamma_RNG wrong parameters" );}
+
+		if ( a < 1. ) { /* GS algorithm for parameters a < 1 */
+			if ( a == 0 )
+				return 0.;
+			e = 1.0 + exp_m1 * a;
+			for ( ;; ) {  //VS repeat
+				p = e * unif_rand();
+				if ( p >= 1.0 ) {
+					x = -log ( ( e - p ) / a );
+					if ( exp_rand() >= ( 1.0 - a ) * log ( x ) )
+						break;
+				}
+				else {
+					x = exp ( log ( p ) / a );
+					if ( exp_rand() >= x )
+						break;
+				}
+			}
+			return scale * x;
+		}
+
+		/* --- a >= 1 : GD algorithm --- */
+
+		/* Step 1: Recalculations of s2, s, d if a has changed */
+		if ( a != aa ) {
+			aa = a;
+			s2 = a - 0.5;
+			s = sqrt ( s2 );
+			d = sqrt32 - s * 12.0;
+		}
+		/* Step 2: t = standard normal deviate,
+		           x = (s,1/2) -normal deviate. */
+
+		/* immediate acceptance (i) */
+		t = norm_rand();
+		x = s + 0.5 * t;
+		ret_val = x * x;
+		if ( t >= 0.0 )
+			return scale * ret_val;
+
+		/* Step 3: u = 0,1 - uniform sample. squeeze acceptance (s) */
+		u = unif_rand();
+		if ( ( d * u ) <= ( t * t * t ) )
+			return scale * ret_val;
+
+		/* Step 4: recalculations of q0, b, si, c if necessary */
+
+		if ( a != aaa ) {
+			aaa = a;
+			r = 1.0 / a;
+			q0 = ( ( ( ( ( ( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r
+			         + q2 ) * r + q1 ) * r;
+
+			/* Approximation depending on size of parameter a */
+			/* The constants in the expressions for b, si and c */
+			/* were established by numerical experiments */
+
+			if ( a <= 3.686 ) {
+				b = 0.463 + s + 0.178 * s2;
+				si = 1.235;
+				c = 0.195 / s - 0.079 + 0.16 * s;
+			}
+			else if ( a <= 13.022 ) {
+				b = 1.654 + 0.0076 * s2;
+				si = 1.68 / s + 0.275;
+				c = 0.062 / s + 0.024;
+			}
+			else {
+				b = 1.77;
+				si = 0.75;
+				c = 0.1515 / s;
+			}
+		}
+		/* Step 5: no quotient test if x not positive */
+
+		if ( x > 0.0 ) {
+			/* Step 6: calculation of v and quotient q */
+			v = t / ( s + s );
+			if ( fabs ( v ) <= 0.25 )
+				q = q0 + 0.5 * t * t * ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v
+				                             + a3 ) * v + a2 ) * v + a1 ) * v;
+			else
+				q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
+
+
+			/* Step 7: quotient acceptance (q) */
+			if ( log ( 1.0 - u ) <= q )
+				return scale * ret_val;
+		}
+
+		for ( ;; ) { //VS repeat
+			/* Step 8: e = standard exponential deviate
+			 *	u =  0,1 -uniform deviate
+			 *	t = (b,si)-double exponential (laplace) sample */
+			e = exp_rand();
+			u = unif_rand();
+			u = u + u - 1.0;
+			if ( u < 0.0 )
+				t = b - si * e;
+			else
+				t = b + si * e;
+			/* Step	 9:  rejection if t < tau(1) = -0.71874483771719 */
+			if ( t >= -0.71874483771719 ) {
+				/* Step 10:	 calculation of v and quotient q */
+				v = t / ( s + s );
+				if ( fabs ( v ) <= 0.25 )
+					q = q0 + 0.5 * t * t *
+					    ( ( ( ( ( ( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v
+					        + a2 ) * v + a1 ) * v;
+				else
+					q = q0 - s * t + 0.25 * t * t + ( s2 + s2 ) * log ( 1.0 + v );
+				/* Step 11:	 hat acceptance (h) */
+				/* (if q not positive go to step 8) */
+				if ( q > 0.0 ) {
+					// TODO: w = expm1(q);
+					w = exp ( q )-1;
+					/*  ^^^^^ original code had approximation with rel.err < 2e-7 */
+					/* if t is rejected sample again at step 8 */
+					if ( ( c * fabs ( u ) ) <= ( w * exp ( e - 0.5 * t * t ) ) )
+						break;
+				}
+			}
+		} /* repeat .. until  `t' is accepted */
+		x = s + 0.5 * t;
+		return scale * x * x;
+	}
+
+
+	bool qr ( const mat &A, mat &R ) {
+		int info;
+		int m = A.rows();
+		int n = A.cols();
+		int lwork = n;
+		int k = std::min ( m, n );
+		vec tau ( k );
+		vec work ( lwork );
+
+		R = A;
+
+		// perform workspace query for optimum lwork value
+		int lwork_tmp = -1;
+		dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
+		          &info );
+		if ( info == 0 ) {
+			lwork = static_cast<int> ( work ( 0 ) );
+			work.set_size ( lwork, false );
+		}
+		dgeqrf_ ( &m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info );
+
+		// construct R
+		for ( int i=0; i<m; i++ )
+			for ( int j=0; j<std::min ( i,n ); j++ )
+				R ( i,j ) = 0;
+
+		return ( info==0 );
+	}
+
 }
-
-
-  bool qr(const mat &A, mat &R)
-  {
-    int info;
-    int m = A.rows();
-    int n = A.cols();
-    int lwork = n;
-    int k = std::min(m, n);
-    vec tau(k);
-    vec work(lwork);
-
-    R = A;
-
-    // perform workspace query for optimum lwork value
-    int lwork_tmp = -1;
-    dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
-            &info);
-    if (info == 0) {
-      lwork = static_cast<int>(work(0));
-      work.set_size(lwork, false);
-    }
-    dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
-
-    // construct R
-    for (int i=0; i<m; i++)
-      for (int j=0; j<std::min(i,n); j++)
-        R(i,j) = 0;
-
-    return (info==0);
-  }
-
-}
