--- /dev/null
+++ b/RGS/src/rgs_dirichlet.c
@@ -0,0 +1,54 @@
+/*
+ * Copyright (C) 2008 Antonio, Fabio Di Narzo <antonio.fabio _at_ gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <Rmath.h>
+#include "rgs_distrib.h"
+#include "rgs_lapack.h"
+
+/*code adapted from 'gsl_ran_dirichlet_lnpdf' in the GSL library */
+SEXP rgs_ddirichlet_work(SEXP x, SEXP in_alpha, SEXP out_ans) {
+	int K = length(x);
+	double *alpha = REAL(in_alpha);
+	double *theta = REAL(x);
+
+	size_t i;
+	double log_p = 0.0;
+	double sum_alpha = 0.0;
+
+	for (i = 0; i < K; i++)
+		log_p += (alpha[i] - 1.0) * log (theta[i]);
+
+	for (i = 0; i < K; i++)
+		sum_alpha += alpha[i];
+
+	log_p += lgammafn(sum_alpha);
+
+	for (i = 0; i < K; i++)
+		log_p -= lgammafn (alpha[i]);
+
+	REAL(out_ans)[0] = log_p;
+	return R_NilValue;
+}
+
+double rgs_ddirichlet (SEXP node) {
+	SEXP ans;
+	PROTECT(ans = allocVector(REALSXP, 1));
+	rgs_ddirichlet_work(RGS_NVALUE(node), RGS_NVALUE(RGS_CPARAMETERS(node, 1)), ans);
+	UNPROTECT(1);
+	return REAL(ans)[0];
+}

--- /dev/null
+++ b/RGS/tests/test_dirichlet.R
@@ -0,0 +1,27 @@
+source("Test.R")
+
+cwrapper <- function(x, alpha) {
+	ans <- 0.0
+	.Call("rgs_ddirichlet_work", as.real(x), as.real(alpha), ans)
+	ans
+}
+
+##Adapted from MCMCpack
+ddirichlet <- function (x, alpha) {
+	dirichlet1 <- function(x, alpha) {
+		logD <- sum(lgamma(alpha)) - lgamma(sum(alpha))
+		s <- sum((alpha - 1) * log(x))
+		exp(sum(s) - logD)
+	}
+	x <- t(x)
+	alpha <- matrix(alpha, ncol = length(alpha), nrow = nrow(x), byrow = TRUE)
+	pd <- vector(length = nrow(x))
+	for (i in 1:nrow(x)) pd[i] <- dirichlet1(x[i, ], alpha[i, ])
+	pd[apply(x, 1, function(z) any(z < 0 | z > 1))] <- 0
+	pd[apply(x, 1, function(z) all.equal(sum(z), 1) != TRUE)] <- 0
+	return(pd)
+}
+
+val.true <- log(ddirichlet(c(0.2, 0.8), c(0.8, 0.2)))
+val.check <- cwrapper(c(0.2, 0.8), c(0.8, 0.2))
+stopifnot(signif(val.true,6) == signif(val.check, 6))