diff --git a/plugins/fn-r/functions.c b/plugins/fn-r/functions.c
index 066bd4b..21483c7 100644
--- a/plugins/fn-r/functions.c
+++ b/plugins/fn-r/functions.c
@@ -1152,6 +1152,35 @@ gnumeric_r_pcauchy (GnmFuncEvalInfo *ei, GnmValue const * const *args)
 /* ------------------------------------------------------------------------- */
 
 
+static GnmFuncHelp const help_r_ptukey[] = {
+	{ GNM_FUNC_HELP_NAME, F_("R.PTUKEY:cumulative distribution function of the Studentized range distribution") },
+	{ GNM_FUNC_HELP_ARG, F_("x:observation") },
+	{ GNM_FUNC_HELP_ARG, F_("nmeans:the number of means") },
+	{ GNM_FUNC_HELP_ARG, F_("df:the number of degrees of freedom of the distribution") },
+	{ GNM_FUNC_HELP_ARG, F_("nranges:the number of ranges") },
+	{ GNM_FUNC_HELP_ARG, F_("lower_tail:if true (the default), the lower tail of the distribution is considered") },
+	{ GNM_FUNC_HELP_ARG, F_("log_p:if true, log of the probability is used") },
+	{ GNM_FUNC_HELP_DESCRIPTION, F_("This function returns the cumulative distribution function of the Studentized range distribution.") },
+	{ GNM_FUNC_HELP_SEEALSO, "" },
+	{ GNM_FUNC_HELP_END }
+};
+
+static GnmValue *
+gnumeric_r_ptukey (GnmFuncEvalInfo *ei, GnmValue const * const *args)
+{
+	gnm_float x = value_get_as_float (args[0]);
+	gnm_float nmeans = value_get_as_float (args[1]);
+	gnm_float df = value_get_as_float (args[2]);
+	gnm_float nranges = args[3] ? value_get_as_float (args[3]) : 1;
+	gboolean lower_tail = args[4] ? value_get_as_checked_bool (args[4]) : TRUE;
+	gboolean log_p = args[5] ? value_get_as_checked_bool (args[5]) : FALSE;
+
+	return value_new_float (ptukey (x, nmeans, df, nranges, lower_tail, log_p));
+}
+
+/* ------------------------------------------------------------------------- */
+
+
 static GnmFuncHelp const help_r_qcauchy[] = {
 	{ GNM_FUNC_HELP_NAME, F_("R.QCAUCHY:probability quantile function of the Cauchy distribution") },
 	{ GNM_FUNC_HELP_ARG, F_("p:probability") },
@@ -1664,6 +1693,13 @@ GnmFuncDescriptor const stat_functions[] = {
 		GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
 	},
 	{
+		"r.ptukey",
+		"fff|fbb",
+		help_r_ptukey,
+		gnumeric_r_ptukey, NULL, NULL, NULL,
+		GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE,
+	},
+	{
 		"r.qcauchy",
 		"fff|bb",
 		help_r_qcauchy,
diff --git a/plugins/fn-r/generate b/plugins/fn-r/generate
index c136a8d..a09cca2 100644
--- a/plugins/fn-r/generate
+++ b/plugins/fn-r/generate
@@ -8,6 +8,7 @@ sub distribution;
 
 my %funcs = ();
 my %argoverride = ();
+my %defaults;
 
 {
     my $of = "of the distribution";
@@ -132,6 +133,15 @@ my %argoverride = ();
 	 'skew-t',
 	 ({ 'n' => "the number of degrees of freedom $of",
 	    @common })];
+
+    $funcs{'ptukey'} =
+	[\&distribution,
+	 'Studentized range',
+	 ({ 'nranges' => "the number of ranges",
+	    'nmeans' => "the number of means",
+	    'df' => "the number of degrees of freedom $of",
+	    @common })];
+    $defaults{'ptukey:nranges'} = 1;
 }
 
 
@@ -297,7 +307,7 @@ sub distribution {
     my $n = 0;
     foreach (@args) {
 	my ($type,$name) = @{ $_ };
-	my $def = undef;
+	my $def = $defaults{"$func:$name"};
 	$def = 'TRUE' if $name eq 'lower_tail';
 	$def = 'FALSE' if $name eq 'give_log' || $name eq 'log_p';
 
diff --git a/plugins/fn-r/plugin.xml.in b/plugins/fn-r/plugin.xml.in
index 901f254..f3eb79f 100644
--- a/plugins/fn-r/plugin.xml.in
+++ b/plugins/fn-r/plugin.xml.in
@@ -44,6 +44,7 @@
 				<function name="r.psnorm"/>
 				<function name="r.pst"/>
 				<function name="r.pt"/>
+				<function name="r.ptukey"/>
 				<function name="r.pweibull"/>
 				<function name="r.qbeta"/>
 				<function name="r.qbinom"/>
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 308fabb..fb22bda 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -5109,8 +5109,432 @@ gnm_float qgeom(gnm_float p, gnm_float prob, gboolean lower_tail, gboolean log_p
 }
 
 /* ------------------------------------------------------------------------ */
+/* Imported src/nmath/ptukey.c from R -- by hand for now.  */
+/*
+ *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 1998       Ross Ihaka
+ *  Copyright (C) 2000--2007 The R Core Team
+ *
+ *  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 2 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, a copy is available at
+ *  http://www.r-project.org/Licenses/
+ *
+ *  SYNOPSIS
+ *
+ *    #include <Rmath.h>
+ *    double ptukey(q, rr, cc, df, lower_tail, log_p);
+ *
+ *  DESCRIPTION
+ *
+ *    Computes the probability that the maximum of rr studentized
+ *    ranges, each based on cc means and with df degrees of freedom
+ *    for the standard error, is less than q.
+ *
+ *    The algorithm is based on that of the reference.
+ *
+ *  REFERENCE
+ *
+ *    Copenhaver, Margaret Diponzio & Holland, Burt S.
+ *    Multiple comparisons of simple effects in
+ *    the two-way analysis of variance with fixed effects.
+ *    Journal of Statistical Computation and Simulation,
+ *    Vol.30, pp.1-15, 1988.
+ */
+
+static gnm_float ptukey_wprob(gnm_float w, gnm_float rr, gnm_float cc)
+{
+/*  wprob() :
+
+	This function calculates probability integral of Hartley's
+	form of the range.
+
+	w     = value of range
+	rr    = no. of rows or groups
+	cc    = no. of columns or treatments
+	ir    = error flag = 1 if pr_w probability > 1
+	pr_w = returned probability integral from (0, w)
+
+	program will not terminate if ir is raised.
+
+	bb = upper limit of legendre integration
+	iMax = maximum acceptable value of integral
+	nleg = order of legendre quadrature
+	ihalf = int ((nleg + 1) / 2)
+	wlar = value of range above which wincr1 intervals are used to
+	       calculate second part of integral,
+	       else wincr2 intervals are used.
+	C1, C2, C3 = values which are used as cutoffs for terminating
+	or modifying a calculation.
+
+	M_1_SQRT_2PI = 1 / sqrt(2 * pi);  from abramowitz & stegun, p. 3.
+	M_SQRT2 = sqrt(2)
+	xleg = legendre 12-point nodes
+	aleg = legendre 12-point coefficients
+ */
+#define nleg	12
+#define ihalf	6
+
+    /* looks like this is suboptimal for gnm_float precision.
+       (see how C1-C3 are used) <MM>
+    */
+    /* const gnm_float iMax  = 1.; not used if = 1*/
+    const static gnm_float C1 = -30.;
+    const static gnm_float C2 = -50.;
+    const static gnm_float C3 = 60.;
+    const static gnm_float bb   = 8.;
+    const static gnm_float wlar = 3.;
+    const static gnm_float wincr1 = 2.;
+    const static gnm_float wincr2 = 3.;
+    const static gnm_float xleg[ihalf] = {
+	GNM_const(0.981560634246719250690549090149),
+	GNM_const(0.904117256370474856678465866119),
+	GNM_const(0.769902674194304687036893833213),
+	GNM_const(0.587317954286617447296702418941),
+	GNM_const(0.367831498998180193752691536644),
+	GNM_const(0.125233408511468915472441369464)
+    };
+    const static gnm_float aleg[ihalf] = {
+	GNM_const(0.047175336386511827194615961485),
+	GNM_const(0.106939325995318430960254718194),
+	GNM_const(0.160078328543346226334652529543),
+	GNM_const(0.203167426723065921749064455810),
+	GNM_const(0.233492536538354808760849898925),
+	GNM_const(0.249147045813402785000562436043)
+    };
+    gnm_float a, ac, pr_w, b, binc, c, cc1,
+	pminus, pplus, qexpo, qsqz, rinsum, wi, wincr, xx;
+    /*LDOUBLE*/ gnm_float blb, bub, einsum, elsum;
+    int j, jj;
+
+
+    qsqz = w * 0.5;
+
+    /* if w >= 16 then the integral lower bound (occurs for c=20) */
+    /* is 0.99999999999995 so return a value of 1. */
+
+    if (qsqz >= bb)
+	return 1.0;
+
+    /* find (f(w/2) - 1) ^ cc */
+    /* (first term in integral of hartley's form). */
+
+    pr_w = 2 * pnorm(qsqz, 0.,1., 1,0) - 1.; /* erf(qsqz / M_SQRT2) */
+    /* if pr_w ^ cc < 2e-22 then set pr_w = 0 */
+    if (pr_w >= exp(C2 / cc))
+	pr_w = gnm_pow(pr_w, cc);
+    else
+	pr_w = 0.0;
+
+    /* if w is large then the second component of the */
+    /* integral is small, so fewer intervals are needed. */
+
+    if (w > wlar)
+	wincr = wincr1;
+    else
+	wincr = wincr2;
+
+    /* find the integral of second term of hartley's form */
+    /* for the integral of the range for equal-length */
+    /* intervals using legendre quadrature.  limits of */
+    /* integration are from (w/2, 8).  two or three */
+    /* equal-length intervals are used. */
+
+    /* blb and bub are lower and upper limits of integration. */
+
+    blb = qsqz;
+    binc = (bb - qsqz) / wincr;
+    bub = blb + binc;
+    einsum = 0.0;
+
+    /* integrate over each interval */
+
+    cc1 = cc - 1.0;
+    for (wi = 1; wi <= wincr; wi++) {
+	elsum = 0.0;
+	a = (gnm_float)(0.5 * (bub + blb));
+
+	/* legendre quadrature with order = nleg */
+
+	b = (gnm_float)(0.5 * (bub - blb));
+
+	for (jj = 1; jj <= nleg; jj++) {
+	    if (ihalf < jj) {
+		j = (nleg - jj) + 1;
+		xx = xleg[j-1];
+	    } else {
+		j = jj;
+		xx = -xleg[j-1];
+	    }
+	    c = b * xx;
+	    ac = a + c;
+
+	    /* if exp(-qexpo/2) < 9e-14, */
+	    /* then doesn't contribute to integral */
+
+	    qexpo = ac * ac;
+	    if (qexpo > C3)
+		break;
+
+	    pplus = 2 * pnorm(ac, 0., 1., 1,0);
+	    pminus= 2 * pnorm(ac, w,  1., 1,0);
+
+	    /* if rinsum ^ (cc-1) < 9e-14, */
+	    /* then doesn't contribute to integral */
+
+	    rinsum = (pplus * 0.5) - (pminus * 0.5);
+	    if (rinsum >= exp(C1 / cc1)) {
+		rinsum = (aleg[j-1] * exp(-(0.5 * qexpo))) * pow(rinsum, cc1);
+		elsum += rinsum;
+	    }
+	}
+	elsum *= (((2.0 * b) * cc) * M_1_SQRT_2PI);
+	einsum += elsum;
+	blb = bub;
+	bub += binc;
+    }
+
+    /* if pr_w ^ rr < 9e-14, then return 0 */
+    pr_w += (gnm_float) einsum;
+    if (pr_w <= exp(C1 / rr))
+	return 0.;
+
+    pr_w = pow(pr_w, rr);
+    if (pr_w >= 1.)/* 1 was iMax was eps */
+	return 1.;
+    return pr_w;
+} /* wprob() */
+
+
+static gnm_float
+R_ptukey(gnm_float q, gnm_float rr, gnm_float cc, gnm_float df,
+	 gboolean lower_tail, gboolean log_p)
+{
+/*  function ptukey() [was qprob() ]:
+
+	q = value of studentized range
+	rr = no. of rows or groups
+	cc = no. of columns or treatments
+	df = degrees of freedom of error term
+	ir[0] = error flag = 1 if wprob probability > 1
+	ir[1] = error flag = 1 if qprob probability > 1
+
+	qprob = returned probability integral over [0, q]
+
+	The program will not terminate if ir[0] or ir[1] are raised.
+
+	All references in wprob to Abramowitz and Stegun
+	are from the following reference:
+
+	Abramowitz, Milton and Stegun, Irene A.
+	Handbook of Mathematical Functions.
+	New York:  Dover publications, Inc. (1970).
+
+	All constants taken from this text are
+	given to 25 significant digits.
+
+	nlegq = order of legendre quadrature
+	ihalfq = int ((nlegq + 1) / 2)
+	eps = max. allowable value of integral
+	eps1 & eps2 = values below which there is
+		      no contribution to integral.
+
+	d.f. <= dhaf:	integral is divided into ulen1 length intervals.  else
+	d.f. <= dquar:	integral is divided into ulen2 length intervals.  else
+	d.f. <= deigh:	integral is divided into ulen3 length intervals.  else
+	d.f. <= dlarg:	integral is divided into ulen4 length intervals.
+
+	d.f. > dlarg:	the range is used to calculate integral.
+
+	M_LN2 = log(2)
+
+	xlegq = legendre 16-point nodes
+	alegq = legendre 16-point coefficients
+
+	The coefficients and nodes for the legendre quadrature used in
+	qprob and wprob were calculated using the algorithms found in:
+
+	Stroud, A. H. and Secrest, D.
+	Gaussian Quadrature Formulas.
+	Englewood Cliffs,
+	New Jersey:  Prentice-Hall, Inc, 1966.
+
+	All values matched the tables (provided in same reference)
+	to 30 significant digits.
+
+	f(x) = .5 + erf(x / sqrt(2)) / 2      for x > 0
+
+	f(x) = erfc( -x / sqrt(2)) / 2	      for x < 0
+
+	where f(x) is standard normal c. d. f.
+
+	if degrees of freedom large, approximate integral
+	with range distribution.
+ */
+#define nlegq	16
+#define ihalfq	8
+
+/*  const gnm_float eps = 1.0; not used if = 1 */
+    const static gnm_float eps1 = -30.0;
+    const static gnm_float eps2 = GNM_const(1.0e-14);
+    const static gnm_float dhaf  = 100.0;
+    const static gnm_float dquar = 800.0;
+    const static gnm_float deigh = 5000.0;
+    const static gnm_float dlarg = 25000.0;
+    const static gnm_float ulen1 = 1.0;
+    const static gnm_float ulen2 = 0.5;
+    const static gnm_float ulen3 = 0.25;
+    const static gnm_float ulen4 = 0.125;
+    const static gnm_float xlegq[ihalfq] = {
+	GNM_const(0.989400934991649932596154173450),
+	GNM_const(0.944575023073232576077988415535),
+	GNM_const(0.865631202387831743880467897712),
+	GNM_const(0.755404408355003033895101194847),
+	GNM_const(0.617876244402643748446671764049),
+	GNM_const(0.458016777657227386342419442984),
+	GNM_const(0.281603550779258913230460501460),
+	GNM_const(0.950125098376374401853193354250e-1)
+    };
+    const static gnm_float alegq[ihalfq] = {
+	GNM_const(0.271524594117540948517805724560e-1),
+	GNM_const(0.622535239386478928628438369944e-1),
+	GNM_const(0.951585116824927848099251076022e-1),
+	GNM_const(0.124628971255533872052476282192),
+	GNM_const(0.149595988816576732081501730547),
+	GNM_const(0.169156519395002538189312079030),
+	GNM_const(0.182603415044923588866763667969),
+	GNM_const(0.189450610455068496285396723208)
+    };
+    gnm_float ans, f2, f21, f2lf, ff4, otsum, qsqz, rotsum, t1, twa1, ulen, wprb;
+    int i, j, jj;
+
+#ifdef IEEE_754
+    if (gnm_isnan(q) || gnm_isnan(rr) || gnm_isnan(cc) || gnm_isnan(df))
+	ML_ERR_return_NAN;
+#endif
+
+    if (q <= 0)
+	return R_DT_0;
+
+    /* df must be > 1 */
+    /* there must be at least two values */
+
+    if (df < 2 || rr < 1 || cc < 2) ML_ERR_return_NAN;
+
+    if(!gnm_finite(q))
+	return R_DT_1;
+
+    if (df > dlarg)
+	return R_DT_val(ptukey_wprob(q, rr, cc));
+
+    /* calculate leading constant */
+
+    f2 = df * 0.5;
+    /* gnm_lgamma(u) = log(gamma(u)) */
+    f2lf = ((f2 * gnm_log(df)) - (df * M_LN2gnum)) - gnm_lgamma(f2);
+    f21 = f2 - 1.0;
+
+    /* integral is divided into unit, half-unit, quarter-unit, or */
+    /* eighth-unit length intervals depending on the value of the */
+    /* degrees of freedom. */
+
+    ff4 = df * 0.25;
+    if	    (df <= dhaf)	ulen = ulen1;
+    else if (df <= dquar)	ulen = ulen2;
+    else if (df <= deigh)	ulen = ulen3;
+    else			ulen = ulen4;
+
+    f2lf += gnm_log(ulen);
+
+    /* integrate over each subinterval */
+
+    ans = 0.0;
+
+    for (i = 1; i <= 50; i++) {
+	otsum = 0.0;
+
+	/* legendre quadrature with order = nlegq */
+	/* nodes (stored in xlegq) are symmetric around zero. */
+
+	twa1 = (2 * i - 1) * ulen;
+
+	for (jj = 1; jj <= nlegq; jj++) {
+	    if (ihalfq < jj) {
+		j = jj - ihalfq - 1;
+		t1 = (f2lf + (f21 * gnm_log(twa1 + (xlegq[j] * ulen))))
+		    - (((xlegq[j] * ulen) + twa1) * ff4);
+	    } else {
+		j = jj - 1;
+		t1 = (f2lf + (f21 * gnm_log(twa1 - (xlegq[j] * ulen))))
+		    + (((xlegq[j] * ulen) - twa1) * ff4);
+
+	    }
+
+	    /* if exp(t1) < 9e-14, then doesn't contribute to integral */
+	    if (t1 >= eps1) {
+		if (ihalfq < jj) {
+		    qsqz = q * gnm_sqrt(((xlegq[j] * ulen) + twa1) * 0.5);
+		} else {
+		    qsqz = q * gnm_sqrt(((-(xlegq[j] * ulen)) + twa1) * 0.5);
+		}
+
+		/* call wprob to find integral of range portion */
+
+		wprb = ptukey_wprob(qsqz, rr, cc);
+		rotsum = (wprb * alegq[j]) * gnm_exp(t1);
+		otsum += rotsum;
+	    }
+	    /* end legendre integral for interval i */
+	    /* L200: */
+	}
+
+	/* if integral for interval i < 1e-14, then stop.
+	 * However, in order to avoid small area under left tail,
+	 * at least  1 / ulen  intervals are calculated.
+	 */
+	if (i * ulen >= 1.0 && otsum <= eps2)
+	    break;
+
+	/* end of interval i */
+	/* L330: */
+
+	ans += otsum;
+    }
+
+    if(otsum > eps2) { /* not converged */
+	ML_ERROR(ME_PRECISION);
+    }
+    if (ans > 1.)
+	ans = 1.;
+    return R_DT_val(ans);
+}
+
+/* Cleaning up done by tools/import-R:  */
+#undef nlegq
+#undef ihalfq
+#undef nleg
+#undef ihalf
+
+/* ------------------------------------------------------------------------ */
 /* --- END MAGIC R SOURCE MARKER --- */
 
+/* Silly order-of-arguments wrapper.  */
+gnm_float
+ptukey(gnm_float x, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail, gboolean log_p)
+{
+	return R_ptukey (x, nranges, nmeans, df, lower_tail, log_p);
+}
+
+
 static gnm_float
 logspace_signed_add (gnm_float logx, gnm_float logabsy, gboolean ypos)
 {
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 038850e..6f6096e 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -132,6 +132,9 @@ gnm_float pcauchy (gnm_float x, gnm_float location, gnm_float scale, gboolean lo
 gnm_float random_exppow_pdf     (gnm_float x, gnm_float a, gnm_float b);
 gnm_float random_laplace_pdf    (gnm_float x, gnm_float a);
 
+/* Studentized range distribution */
+gnm_float ptukey(gnm_float x, gnm_float nmeans, gnm_float df, gnm_float nranges, gboolean lower_tail, gboolean log_p);
+
 /* ------------------------------------------------------------------------- */
 
 typedef gnm_float (*GnmPFunc) (gnm_float x, const gnm_float shape[],
