Index: /bdm/math/libDC.h
===================================================================
--- /bdm/math/libDC.h (revision 51)
+++ /bdm/math/libDC.h (revision 75)
@@ -33,5 +33,5 @@
 	public:
 		/*!
-		 * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$.
+		 * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$.
 		 * @param  v Vector forming the outer product to be added
 		 * @param  w weight of updating; can be negative
@@ -46,10 +46,10 @@
 		virtual mat to_mat() =0;
 
-		/*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$
+		/*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$
 		@param C multiplying matrix,
 		*/
 		virtual void mult_sym ( const mat &C ) =0;
 		
-		/*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix $C$, i.e. $V = C'*V*C$
+		/*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$
 		@param C multiplying matrix,
 		*/
@@ -64,5 +64,5 @@
 
 		/*!
-		\brief Multiplies square root of $V$ by vector $x$.
+		\brief Multiplies square root of \f$V\f$ by vector \f$x\f$.
 
 		Used e.g. in generating normal samples.
@@ -71,8 +71,14 @@
 
 		/*!
-		\brief Evaluates quadratic form $x= v'*V*v$;
+		\brief Evaluates quadratic form \f$x= v'*V*v\f$;
 
 		*/
 		virtual double qform (const vec &v ) const =0;
+
+		/*!
+		\brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$;
+
+		*/
+		virtual double invqform (const vec &v ) const =0;
 
 //	//! easy version of the
@@ -112,7 +118,7 @@
 		void mult_sym ( const mat &C);
 		void mult_sym_t ( const mat &C);
-		//! store result of \c mult_sym in external matrix $U$
+		//! store result of \c mult_sym in external matrix \f$U\f$
 		void mult_sym ( const mat &C, fsqmat &U) const;
-		//! store result of \c mult_sym_t in external matrix $U$
+		//! store result of \c mult_sym_t in external matrix \f$U\f$
 		void mult_sym_t ( const mat &C, fsqmat &U) const;
 		void clear();
@@ -140,4 +146,5 @@
 		double logdet() const {return log ( det ( M ) );};
 		double qform (const  vec &v ) const {return ( v* ( M*v ) );};
+		double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );};
 		vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};
 
@@ -163,5 +170,5 @@
 /*! \brief Matrix stored in LD form, (typically known as UD) 
 
-Matrix is decomposed as follows: \f[M = L'DL\f] where only $L$ and $D$ matrices are stored.
+Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored.
 All inplace operations modifies only these and the need to compose and decompose the matrix is avoided.
 */
@@ -194,4 +201,5 @@
 		double logdet() const;
 		double qform (const vec &v ) const;
+		double invqform (const vec &v ) const;
 //	sqmat& operator -= ( const sqmat & ld2 );
 		void clear();
@@ -205,5 +213,5 @@
 		virtual void inv ( ldmat &Inv ) const;
 
-		/*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.
+		/*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class.
 		@param C matrix to multiply with
 		@param U a space where the inverse is stored.
@@ -211,5 +219,5 @@
 		void mult_sym ( const mat &C, ldmat &U) const;
 
-		/*! \brief Symmetric multiplication of $U$ by a transpose of a general matrix $C$, result of which is stored in the current class.
+		/*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class.
 		@param C matrix to multiply with
 		@param U a space where the inverse is stored.
@@ -218,7 +226,7 @@
 
 
-		/*! \brief Transforms general $A'D0 A$ into pure $L'DL$
-
-		The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$
+		/*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$
+
+		The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$
 		@param A general matrix
 		@param D0 general vector
@@ -245,7 +253,7 @@
 
 	protected:
-		//! Positive vector $D$
+		//! Positive vector \f$D\f$
 		vec D;
-		//! Lower-triangular matrix $L$
+		//! Lower-triangular matrix \f$L\f$
 		mat L;
 
Index: /bdm/math/libDC.cpp
===================================================================
--- /bdm/math/libDC.cpp (revision 37)
+++ /bdm/math/libDC.cpp (revision 75)
@@ -162,9 +162,25 @@
 	double x = 0.0, sum;
 	int i,j;
+	vec s(v.length());
+	vec S=L*v;
 
 	for ( i=0; i<D.length(); i++ ) { //rows of L
 		sum = 0.0;
-		for ( j=i; j<=i; j++ ){sum+=L( i,j )*v( i );}
+		for ( j=0; j<=i; j++ ){sum+=L( i,j )*v( j );}
 		x +=D( i )*sum*sum;
+		s(i)=sum;
+	};
+	return x;
+}
+
+double ldmat::invqform( const vec &v ) const {
+	double x = 0.0;
+	int i;
+	vec pom(v.length());
+	
+	backward_substitution(L.T(),v,pom);
+	
+	for ( i=0; i<D.length(); i++ ) { //rows of L
+		x +=pom(i)*pom(i)/D(i);
 	};
 	return x;
@@ -298,5 +314,5 @@
 			j = i + k; //change in .m 1+1=2, here 0+0+1=1
 			s = L( j, i );
-			for ( m = i + 1; m < ( j - 1 ); m++ ) {
+			for ( m = i + 1; m < ( j ); m++ ) {
 				s += L( m, i ) * Il( j, m );
 			}
