From patchwork Fri Apr  8 20:36:20 2016
Content-Type: text/plain; charset="utf-8"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
X-Patchwork-Submitter: Ed Smith-Rowland <3dw4rd@verizon.net>
X-Patchwork-Id: 608181
Return-Path: 
 <gcc-patches-return-424585-incoming=patchwork.ozlabs.org@gcc.gnu.org>
X-Original-To: incoming@patchwork.ozlabs.org
Delivered-To: patchwork-incoming@bilbo.ozlabs.org
Received: from sourceware.org (server1.sourceware.org [209.132.180.131])
	(using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256
	bits)) (No client certificate requested)
	by ozlabs.org (Postfix) with ESMTPS id 3qhWVQ3JLxz9t4n
	for <incoming@patchwork.ozlabs.org>;
	Sat,  9 Apr 2016 06:37:33 +1000 (AEST)
Authentication-Results: ozlabs.org; dkim=pass (1024-bit key;
	unprotected) header.d=gcc.gnu.org header.i=@gcc.gnu.org
	header.b=v3FcBXhg; dkim-atps=neutral
DomainKey-Signature: a=rsa-sha1; c=nofws; d=gcc.gnu.org; h=list-id
	:list-unsubscribe:list-archive:list-post:list-help:sender:to
	:from:subject:message-id:date:mime-version:content-type; q=dns;
	s=default; b=BLdqZGps1wq2cSfFjc7M4rY0BLHzVmZ55rDlSc5eaCqcmOdMT0
	zjJF88X6duPWGnqNcooTRFLCSZ4TgqU5qB/auhoOuwGG665AKtMHp9XaTu6IILLO
	11TAa+/s2WP85UVa0Ni9KvQK6Ylcm1+8zqBUvLX9lvVpczTzQIw0Ym+f0=
DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=gcc.gnu.org; h=list-id
	:list-unsubscribe:list-archive:list-post:list-help:sender:to
	:from:subject:message-id:date:mime-version:content-type; s=
	default; bh=Y5Kl035L75QmZDfaOTxPX0U3OOE=; b=v3FcBXhg66IuOpMfgZGc
	xFFEXPcQHpwe09wxRnJ+FD7dQUpCK12p2NmCueqFw480yVt10iXt0gWN93msPnu2
	tBuiUKTZwZlMH/5I0b4QQP0gUAnFgI+p/K9LF6imYGrsrF3MHxVD2P7G8YoA9Bcs
	+tJHHugQ3NqXeVMCdSD914c=
Received: (qmail 96532 invoked by alias); 8 Apr 2016 20:37:24 -0000
Mailing-List: contact gcc-patches-help@gcc.gnu.org; run by ezmlm
Precedence: bulk
List-Id: <gcc-patches.gcc.gnu.org>
List-Unsubscribe: 
 <mailto:gcc-patches-unsubscribe-incoming=patchwork.ozlabs.org@gcc.gnu.org>
List-Archive: <http://gcc.gnu.org/ml/gcc-patches/>
List-Post: <mailto:gcc-patches@gcc.gnu.org>
List-Help: <mailto:gcc-patches-help@gcc.gnu.org>
Sender: gcc-patches-owner@gcc.gnu.org
Delivered-To: mailing list gcc-patches@gcc.gnu.org
Received: (qmail 96499 invoked by uid 89); 8 Apr 2016 20:37:22 -0000
Authentication-Results: sourceware.org; auth=none
X-Virus-Found: No
X-Spam-SWARE-Status: No, score=0.2 required=5.0 tests=AWL, BAYES_50,
	NO_RDNS_DOTCOM_HELO, RCVD_IN_DNSWL_LOW, RP_MATCHES_RCVD,
	SPF_PASS autolearn=ham version=3.3.2 spammy=Electronic, lm,
	Beta, irene
X-Spam-User: qpsmtpd, 2 recipients
X-HELO: vms173023pub.verizon.net
Received: from vms173023pub.verizon.net (HELO vms173023pub.verizon.net)
	(206.46.173.23) by sourceware.org
	(qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-SHA encrypted)
	ESMTPS; Fri, 08 Apr 2016 20:37:12 +0000
Received: from vz-proxy-m003.mx.aol.com ([64.236.83.10]) by
	vms173023.mailsrvcs.net (Oracle Communications Messaging
	Server 7.0.5.32.0 64bit (built Jul 16 2014)) with ESMTPA id
	<0O5C00HDQ18LFO90@vms173023.mailsrvcs.net>;
	Fri, 08 Apr 2016 15:36:22 -0500 (CDT)
X-CMAE-Score: 0
X-CMAE-Analysis: v=2.1 cv=EdU1O6SC c=1 sm=1
	tr=0	a=pKq81LsWol7b8sI3Hpopew==:117
	a=hTZS6uLwhHpx5LOQAS/v3KalXBI=:19	a=kziv93cY1bsA:10
	a=r77TgQKjGQsHNAKrUKIA:9
	a=pVbYQFakkQVIKPZwQHIA:9	a=QEXdDO2ut3YA:10 a=U9JUAQm4G5AA:10
	a=o1OHuDzbAAAA:8	a=LPZBJLCeC_td05IT8QsA:9 a=sXYZbZiifVsA:10
	a=FoLSdCaC-3YA:10	a=b4LDLZbEAAAA:8 a=mDV3o1hIAAAA:8
	a=AIPdC31oAAAA:8 a=U1ntXoZbAAAA:8	a=PYnjg3YJAAAA:8
	a=dEJNbJKXEDIeK__103IA:9
	a=Ch0kvYTPxYI1iL6r:21	a=iTR2pfzJBJbwoCfO:21
	a=xM6gWm99s215vdAj:21 a=NWVoK91CQyQA:10
Received: by 71.121.253.69 with SMTP id 76178f38;
	Fri, 08 Apr 2016 20:36:21 GMT
To: "libstdc++@gcc.gnu.org" <libstdc++@gcc.gnu.org>,
	gcc-patches <gcc-patches@gcc.gnu.org>,
	Jonathan Wakely <jwakely.gcc@gmail.com>
From: Ed Smith-Rowland <3dw4rd@verizon.net>
Subject: Document C++ Special Math Functions.
Message-id: <57081644.3090904@verizon.net>
Date: Fri, 08 Apr 2016 16:36:20 -0400
User-Agent: Mozilla/5.0 (X11; Linux x86_64;
	rv:38.0) Gecko/20100101 Thunderbird/38.6.0
MIME-version: 1.0
Content-type: multipart/mixed; boundary=------------040901030905070600010708

I wanted to ship the TR29124 special math functions in libstdc++ with 
some documentation.
More could be done but this covers the function definitions, argument 
ranges, template parms, and arguments.
There is a little mainpage with some overview, history, biblio, and links.
2016-04-08  Edward Smith-Rowland  <3dw4rd@verizon.net>

	Document C++19/TR29124 C++ Special Math Functions.
	* include/bits/specfun.h: Add Doxygen markup.

Index: include/bits/specfun.h
===================================================================
--- include/bits/specfun.h	(revision 234525)
+++ include/bits/specfun.h	(working copy)
@@ -75,16 +75,182 @@
    * @{
    */
 
+  /**
+   * @mainpage Mathematical Special Functions
+   *
+   * @section intro Introduction and History
+   * The first significant library upgrade on the road to C++2011,
+   * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2005/n1836.pdf">
+   * TR1</a>, included a set of 23 mathematical functions that significntly
+   * extended the standard trancendental functions inherited from C and declared
+   * in @<cmath@>.
+   *
+   * Although most components from TR1 were eventually adopted for C++11 these
+   * math function were left behind out of concern for implementability.
+   * The math functions were published as a separate international standard
+   * <a href="http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2010/n3060.pdf">
+   * IS 29124 - Extensions to the C++ Library to Support Mathematical Special
+   * Functions</a>.
+   *
+   * For C++17 these functions were incorporated into the main standard.
+   *
+   * @section contents Contents
+   * The folowing functions are implemented in namespace @c std:
+   * - @ref assoc_laguerre "assoc_laguerre - Associated Laguerre functions"
+   * - @ref assoc_legendre "assoc_legendre - Associated Legendre functions"
+   * - @ref beta "beta - Beta functions"
+   * - @ref comp_ellint_1 "comp_ellint_1 - Complete elliptic functions of the first kind"
+   * - @ref comp_ellint_2 "comp_ellint_2 - Complete elliptic functions of the second kind"
+   * - @ref comp_ellint_3 "comp_ellint_3 - Complete elliptic functions of the third kind"
+   * - @ref cyl_bessel_i "cyl_bessel_i - Regular modified cylindrical Bessel functions"
+   * - @ref cyl_bessel_j "cyl_bessel_j - Cylindrical Bessel functions of the first kind"
+   * - @ref cyl_bessel_k "cyl_bessel_k - Irregular modified cylindrical Bessel functions"
+   * - @ref cyl_neumann "cyl_neumann - Cylindrical Neumann functions or Cylindrical Bessel functions of the second kind"
+   * - @ref ellint_1 "ellint_1 - Incomplete elliptic functions of the first kind"
+   * - @ref ellint_2 "ellint_2 - Incomplete elliptic functions of the second kind"
+   * - @ref ellint_3 "ellint_3 - Incomplete elliptic functions of the third kind"
+   * - @ref expint "expint - The exponential integral"
+   * - @ref hermite "hermite - Hermite polynomials"
+   * - @ref laguerre "laguerre - Laguerre functions"
+   * - @ref legendre "legendre - Legendre polynomials"
+   * - @ref riemann_zeta "riemann_zeta - The Riemann zeta function"
+   * - @ref sph_bessel "sph_bessel - Spherical Bessel functions"
+   * - @ref sph_legendre "sph_legendre - Spherical Legendre functions"
+   * - @ref sph_neumann "sph_neumann - Spherical Neumann functions"
+   *
+   * The hypergeometric functions were stricken from the TR29124 and C++17
+   * versions of this math library because of implementation concerns.
+   * However, since they were in the TR1 version and since they are popular
+   * we kept them as an extension in namespace @c __gnu_cxx:
+   * - @ref conf_hyperg "conf_hyperg - Confluent hypergeometric functions"
+   * - @ref hyperg "hyperg - Hypergeometric functions"
+   *
+   * @section general General Features
+   *
+   * @subsection "Argument Promotion"
+   * The arguments suppled to the non-suffixed functions will be promoted
+   * according to the following rules:
+   * 1. If any argument intended to be floating opint is given an integral value
+   * That integral value is promoted to double.
+   * 2. All floating point arguments are promoted up to the largest floating
+   *    point precision among them.
+   *
+   * @subsection NaN NaN Arguments
+   * If any of the floating point arguments supplied to these functions is
+   * invalid or NaN (std::numeric_limits<Tp>::quiet_NaN),
+   * the value NaN is returned.
+   *
+   * @section impl Implementation
+   *
+   * We strive to implement the underlying math with type generic algorithms
+   * to the greatest extent possible.  In practice, the function are thin
+   * wrappers that dispatch to function templates. Type dependence is
+   * controlled with std::numeric_limits and functions thereof.
+   *
+   * We don't promote *c float to *c double or *c double to <tt>long double</tt>
+   * reflexively.  The goal is for float functions to operate more quickly,
+   * at the cost of float accuracy and possibly a smaller domain of validity.
+   * Similaryly, <tt>long double</tt> should give you more dynamic range
+   * and slightly more pecision than @c double on many systems.
+   *
+   * @section testing Testing
+   *
+   * These functions have been tested against equivalent implementations
+   * from the <a href="http://www.gnu.org/software/gsl">
+   * Gnu Scientific Library, GSL</a> and
+   * <a href="http://www.boost.org/doc/libs/1_60_0/libs/math/doc/html/index.html>Boost</a>
+   * and the ratio
+   * @f[
+   *   \frac{|f - f_{test}|}{|f_{test}|}
+   * @f]
+   * is generally found to be within 10^-15 for 64-bit double on linux-x86_64 systems
+   * over most of the ranges of validity.
+   * 
+   * @todo Provide accuracy comparisons on a per-function basis for a small
+   *       number of targets.
+   *
+   * @section bibliography General Bibliography
+   *
+   * @see Abramowitz and Stegun: Handbook of Mathematical Functions,
+   * with Formulas, Graphs, and Mathematical Tables
+   * Edited by Milton Abramowitz and Irene A. Stegun,
+   * National Bureau of Standards  Applied Mathematics Series - 55
+   * Issued June 1964, Tenth Printing, December 1972, with corrections
+   * Electronic versions of A&S abound including both pdf and navigable html.
+   * @see for example  http://people.math.sfu.ca/~cbm/aands/
+   *
+   * @see The old A&S has been redone as the
+   * NIST Digital Library of Mathematical Functions: http://dlmf.nist.gov/
+   * This version is far more navigable and includes more recent work.
+   *
+   * @see An Atlas of Functions: with Equator, the Atlas Function Calculator
+   * 2nd Edition, by Oldham, Keith B., Myland, Jan, Spanier, Jerome
+   *
+   * @see Asymptotics and Special Functions by Frank W. J. Olver,
+   * Academic Press, 1974
+   *
+   * @see Numerical Recipes in C, The Art of Scientific Computing,
+   * by William H. Press, Second Ed., Saul A. Teukolsky,
+   * William T. Vetterling, and Brian P. Flannery,
+   * Cambridge University Press, 1992
+   *
+   * @see The Special Functions and Their Approximations: Volumes 1 and 2,
+   * by Yudell L. Luke, Academic Press, 1969
+   */
+
   // Associated Laguerre polynomials
 
+  /**
+   * Return the associated Laguerre polynomial of order @c n,
+   * degree @c m: @f$ L_n^m(x) @f$ for @c float argument.
+   *
+   * @see assoc_laguerre for more details.
+   */
   inline float
   assoc_laguerref(unsigned int __n, unsigned int __m, float __x)
   { return __detail::__assoc_laguerre<float>(__n, __m, __x); }
 
+  /**
+   * Return the associated Laguerre polynomial of order @c n,
+   * degree @c m: @f$ L_n^m(x) @f$.
+   *
+   * @see assoc_laguerre for more details.
+   */
   inline long double
   assoc_laguerrel(unsigned int __n, unsigned int __m, long double __x)
   { return __detail::__assoc_laguerre<long double>(__n, __m, __x); }
 
+  /**
+   * Return the associated Laguerre polynomial of nonnegative order @c n,
+   * nonnegative degree @c m and real argument @c x: @f$ L_n^m(x) @f$.
+   *
+   * The associated Laguerre function of real degree @f$ \alpha @f$,
+   * @f$ L_n^\alpha(x) @f$, is defined by
+   * @f[
+   * 	 L_n^\alpha(x) = \frac{(\alpha + 1)_n}{n!}
+   * 			 {}_1F_1(-n; \alpha + 1; x)
+   * @f]
+   * where @f$ (\alpha)_n @f$ is the Pochhammer symbol and
+   * @f$ {}_1F_1(a; c; x) @f$ is the confluent hypergeometric function.
+   *
+   * The associated Laguerre polynomial is defined for integral
+   * degree @f$ \alpha = m @f$ by:
+   * @f[
+   * 	 L_n^m(x) = (-1)^m \frac{d^m}{dx^m} L_{n + m}(x)
+   * @f]
+   * where the Laguerre polynomial is defined by:
+   * @f[
+   * 	 L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x})
+   * @f]
+   * and @f$ x >= 0 @f$.
+   * @see laguerre for details of the Laguerre function of degree @c n
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param __n The order of the Laguerre function, <tt>__n >= 0</tt>.
+   * @param __m The degree of the Laguerre function, <tt>__m >= 0</tt>.
+   * @param __x The argument of the Laguerre function, <tt>__x >= 0</tt>.
+   * @throw std::domain_error if <tt>__x < 0</tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     assoc_laguerre(unsigned int __n, unsigned int __m, _Tp __x)
@@ -95,14 +261,42 @@
 
   // Associated Legendre functions
 
+  /**
+   * Return the associated Legendre function of degree @c l and order @c m
+   * for @c float argument.
+   *
+   * @see assoc_legendre for more details.
+   */
   inline float
   assoc_legendref(unsigned int __l, unsigned int __m, float __x)
   { return __detail::__assoc_legendre_p<float>(__l, __m, __x); }
 
+  /**
+   * Return the associated Legendre function of degree @c l and order @c m.
+   *
+   * @see assoc_legendre for more details.
+   */
   inline long double
   assoc_legendrel(unsigned int __l, unsigned int __m, long double __x)
   { return __detail::__assoc_legendre_p<long double>(__l, __m, __x); }
 
+
+  /**
+   * Return the associated Legendre function of degree @c l and order @c m.
+   *
+   * The associated Legendre function is derived from the Legendre function
+   * @f$ P_l(x) @f$ by the Rodrigues formula:
+   * @f[
+   *   P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)
+   * @f]
+   * @see legendre for details of the Legendre function of degree @c l
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __l  The degree <tt>__l >= 0</tt>.
+   * @param  __m  The order <tt>__m <= l</tt>.
+   * @param  __x  The argument, <tt>abs(__x) <= 1</tt>.
+   * @throw std::domain_error if <tt>abs(__x) > 1</tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     assoc_legendre(unsigned int __l, unsigned int __m, _Tp __x)
@@ -113,32 +307,89 @@
 
   // Beta functions
 
+  /**
+   * Return the beta function, @f$ B(a,b) @f$, for @c float parameters @c a, @c b.
+   *
+   * @see beta for more details.
+   */
   inline float
-  betaf(float __x, float __y)
-  { return __detail::__beta<float>(__x, __y); }
+  betaf(float __a, float __b)
+  { return __detail::__beta<float>(__a, __b); }
 
+  /**
+   * Return the beta function, @f$B(a,b)@f$, for long double
+   * parameters @c a, @c b.
+   *
+   * @see beta for more details.
+   */
   inline long double
-  betal(long double __x, long double __y)
-  { return __detail::__beta<long double>(__x, __y); }
+  betal(long double __a, long double __b)
+  { return __detail::__beta<long double>(__a, __b); }
 
-  template<typename _Tpx, typename _Tpy>
-    inline typename __gnu_cxx::__promote_2<_Tpx, _Tpy>::__type
-    beta(_Tpx __x, _Tpy __y)
+  /**
+   * Return the beta function, @f$B(a,b)@f$, for real parameters @c a, @c b.
+   *
+   * The beta function is defined by
+   * @f[
+   *   B(a,b) = \int_0^1 t^{a - 1} (1 - t)^{b - 1} dt
+   *          = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
+   * @f]
+   * where @f$ x > 0 @f$ and @f$ y > 0 @f$
+   *
+   * @tparam _Tpa The floating-point type of the parameter @c __a.
+   * @tparam _Tpb The floating-point type of the parameter @c __b.
+   * @param __a The first argument of the beta function, <tt> __a > 0 </tt>.
+   * @param __b The second argument of the beta function, <tt> __b > 0 </tt>.
+   * @throw std::domain_error if <tt> __a < 0 </tt> or <tt> __b < 0 </tt>.
+   */
+  template<typename _Tpa, typename _Tpb>
+    inline typename __gnu_cxx::__promote_2<_Tpa, _Tpb>::__type
+    beta(_Tpa __a, _Tpb __b)
     {
-      typedef typename __gnu_cxx::__promote_2<_Tpx, _Tpy>::__type __type;
-      return __detail::__beta<__type>(__x, __y);
+      typedef typename __gnu_cxx::__promote_2<_Tpa, _Tpb>::__type __type;
+      return __detail::__beta<__type>(__a, __b);
     }
 
   // Complete elliptic integrals of the first kind
 
+  /**
+   * Return the complete elliptic integral of the first kind @f$ E(k) @f$
+   * for @c float modulus @c k.
+   *
+   * @see comp_ellint_1 for details.
+   */
   inline float
   comp_ellint_1f(float __k)
   { return __detail::__comp_ellint_1<float>(__k); }
 
+  /**
+   * Return the complete elliptic integral of the first kind @f$ E(k) @f$
+   * for long double modulus @c k.
+   *
+   * @see comp_ellint_1 for details.
+   */
   inline long double
   comp_ellint_1l(long double __k)
   { return __detail::__comp_ellint_1<long double>(__k); }
 
+  /**
+   * Return the complete elliptic integral of the first kind
+   * @f$ K(k) @f$ for real modulus @c k.
+   *
+   * The complete elliptic integral of the first kind is defined as
+   * @f[
+   *   K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta}
+   * 					     {\sqrt{1 - k^2 sin^2\theta}}
+   * @f]
+   * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the
+   * first kind and the modulus @f$ |k| <= 1 @f$.
+   * @see ellint_1 for details of the incomplete elliptic function
+   * of the first kind.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @param  __k  The modulus, <tt> abs(__k) <= 1 </tt>
+   * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     comp_ellint_1(_Tp __k)
@@ -149,14 +400,43 @@
 
   // Complete elliptic integrals of the second kind
 
+  /**
+   * Return the complete elliptic integral of the second kind @f$ E(k) @f$
+   * for @c float modulus @c k.
+   *
+   * @see comp_ellint_2 for details.
+   */
   inline float
   comp_ellint_2f(float __k)
   { return __detail::__comp_ellint_2<float>(__k); }
 
+  /**
+   * Return the complete elliptic integral of the second kind @f$ E(k) @f$
+   * for long double modulus @c k.
+   *
+   * @see comp_ellint_2 for details.
+   */
   inline long double
   comp_ellint_2l(long double __k)
   { return __detail::__comp_ellint_2<long double>(__k); }
 
+  /**
+   * Return the complete elliptic integral of the second kind @f$ E(k) @f$
+   * for real modulus @c k.
+   *
+   * The complete elliptic integral of the second kind is defined as
+   * @f[
+   *   E(k) = E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
+   * @f]
+   * where @f$ E(k,\phi) @f$ is the incomplete elliptic integral of the
+   * second kind and the modulus @f$ |k| <= 1 @f$.
+   * @see ellint_2 for details of the incomplete elliptic function
+   * of the second kind.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @param  __k  The modulus, @c abs(__k) <= 1
+   * @throw std::domain_error if @c abs(__k) > 1.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     comp_ellint_2(_Tp __k)
@@ -167,14 +447,47 @@
 
   // Complete elliptic integrals of the third kind
 
+  /**
+   * @brief Return the complete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu) @f$ for @c float modulus @c k.
+   *
+   * @see comp_ellint_3 for details.
+   */
   inline float
   comp_ellint_3f(float __k, float __nu)
   { return __detail::__comp_ellint_3<float>(__k, __nu); }
 
+  /**
+   * @brief Return the complete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu) @f$ for <tt>long double</tt> modulus @c k.
+   *
+   * @see comp_ellint_3 for details.
+   */
   inline long double
   comp_ellint_3l(long double __k, long double __nu)
   { return __detail::__comp_ellint_3<long double>(__k, __nu); }
 
+  /**
+   * Return the complete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ for real modulus @c k.
+   *
+   * The complete elliptic integral of the third kind is defined as
+   * @f[
+   *   \Pi(k,\nu) = \Pi(k,\nu,\pi/2) = \int_0^{\pi/2}
+   * 		     \frac{d\theta}
+   * 		   {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}}
+   * @f]
+   * where @f$ \Pi(k,\nu,\phi) @f$ is the incomplete elliptic integral of the
+   * second kind and the modulus @f$ |k| <= 1 @f$.
+   * @see ellint_3 for details of the incomplete elliptic function
+   * of the third kind.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @tparam _Tpn The floating-point type of the argument @c __nu.
+   * @param  __k  The modulus, @c abs(__k) <= 1
+   * @param  __nu  The argument
+   * @throw std::domain_error if @c abs(__k) > 1.
+   */
   template<typename _Tp, typename _Tpn>
     inline typename __gnu_cxx::__promote_2<_Tp, _Tpn>::__type
     comp_ellint_3(_Tp __k, _Tpn __nu)
@@ -185,14 +498,42 @@
 
   // Regular modified cylindrical Bessel functions
 
+  /**
+   * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
+   * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_i for setails.
+   */
   inline float
   cyl_bessel_if(float __nu, float __x)
   { return __detail::__cyl_bessel_i<float>(__nu, __x); }
 
+  /**
+   * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
+   * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_i for setails.
+   */
   inline long double
   cyl_bessel_il(long double __nu, long double __x)
   { return __detail::__cyl_bessel_i<long double>(__nu, __x); }
 
+  /**
+   * Return the regular modified Bessel function @f$ I_{\nu}(x) @f$
+   * for real order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * The regular modified cylindrical Bessel function is:
+   * @f[
+   *  I_{\nu}(x) = i^{-\nu}J_\nu(ix) = \sum_{k=0}^{\infty}
+   * 		\frac{(x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}
+   * @f]
+   *
+   * @tparam _Tpnu The floating-point type of the order @c __nu.
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __nu  The order
+   * @param  __x   The argument, <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tpnu, typename _Tp>
     inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
     cyl_bessel_i(_Tpnu __nu, _Tp __x)
@@ -203,14 +544,42 @@
 
   // Cylindrical Bessel functions (of the first kind)
 
+  /**
+   * Return the Bessel function of the first kind @f$ J_{\nu}(x) @f$
+   * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_j for setails.
+   */
   inline float
   cyl_bessel_jf(float __nu, float __x)
   { return __detail::__cyl_bessel_j<float>(__nu, __x); }
 
+  /**
+   * Return the Bessel function of the first kind @f$ J_{\nu}(x) @f$
+   * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_j for setails.
+   */
   inline long double
   cyl_bessel_jl(long double __nu, long double __x)
   { return __detail::__cyl_bessel_j<long double>(__nu, __x); }
 
+  /**
+   * Return the Bessel function @f$ J_{\nu}(x) @f$ of real order @f$ \nu @f$
+   * and argument @f$ x >= 0 @f$.
+   *
+   * The cylindrical Bessel function is:
+   * @f[
+   *    J_{\nu}(x) = \sum_{k=0}^{\infty}
+   *              \frac{(-1)^k (x/2)^{\nu + 2k}}{k!\Gamma(\nu+k+1)}
+   * @f]
+   *
+   * @tparam _Tpnu The floating-point type of the order @c __nu.
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __nu  The order
+   * @param  __x   The argument, <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tpnu, typename _Tp>
     inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
     cyl_bessel_j(_Tpnu __nu, _Tp __x)
@@ -221,14 +590,48 @@
 
   // Irregular modified cylindrical Bessel functions
 
+  /**
+   * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
+   * for @c float order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_k for setails.
+   */
   inline float
   cyl_bessel_kf(float __nu, float __x)
   { return __detail::__cyl_bessel_k<float>(__nu, __x); }
 
+  /**
+   * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
+   * for <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * @see cyl_bessel_k for setails.
+   */
   inline long double
   cyl_bessel_kl(long double __nu, long double __x)
   { return __detail::__cyl_bessel_k<long double>(__nu, __x); }
 
+  /**
+   * Return the irregular modified Bessel function @f$ K_{\nu}(x) @f$
+   * of real order @f$ \nu @f$ and argument @f$ x @f$.
+   *
+   * The irregular modified Bessel function is defined by:
+   * @f[
+   * 	K_{\nu}(x) = \frac{\pi}{2}
+   * 		     \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin \nu\pi}
+   * @f]
+   * where for integral @f$ \nu = n @f$ a limit is taken:
+   * @f$ lim_{\nu \to n} @f$.
+   * For negative argument we have simply:
+   * @f[
+   * 	K_{-\nu}(x) = K_{\nu}(x)
+   * @f]
+   *
+   * @tparam _Tpnu The floating-point type of the order @c __nu.
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __nu  The order
+   * @param  __x   The argument, <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tpnu, typename _Tp>
     inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
     cyl_bessel_k(_Tpnu __nu, _Tp __x)
@@ -239,14 +642,44 @@
 
   // Cylindrical Neumann functions
 
+  /**
+   * Return the Neumann function @f$ N_{\nu}(x) @f$
+   * of @c float order @f$ \nu @f$ and argument @f$ x @f$.
+   *
+   * @see cyl_neumann for setails.
+   */
   inline float
   cyl_neumannf(float __nu, float __x)
   { return __detail::__cyl_neumann_n<float>(__nu, __x); }
 
+  /**
+   * Return the Neumann function @f$ N_{\nu}(x) @f$
+   * of <tt>long double</tt> order @f$ \nu @f$ and argument @f$ x @f$.
+   *
+   * @see cyl_neumann for setails.
+   */
   inline long double
   cyl_neumannl(long double __nu, long double __x)
   { return __detail::__cyl_neumann_n<long double>(__nu, __x); }
 
+  /**
+   * Return the Neumann function @f$ N_{\nu}(x) @f$
+   * of real order @f$ \nu @f$ and argument @f$ x >= 0 @f$.
+   *
+   * The Neumann function is defined by:
+   * @f[
+   *    N_{\nu}(x) = \frac{J_{\nu}(x) \cos \nu\pi - J_{-\nu}(x)}
+   *                      {\sin \nu\pi}
+   * @f]
+   * where @f$ x >= 0 @f$ and for integral order @f$ \nu = n @f$
+   * a limit is taken: @f$ lim_{\nu \to n} @f$.
+   *
+   * @tparam _Tpnu The floating-point type of the order @c __nu.
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __nu  The order
+   * @param  __x   The argument, <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tpnu, typename _Tp>
     inline typename __gnu_cxx::__promote_2<_Tpnu, _Tp>::__type
     cyl_neumann(_Tpnu __nu, _Tp __x)
@@ -257,14 +690,44 @@
 
   // Incomplete elliptic integrals of the first kind
 
+  /**
+   * Return the incomplete elliptic integral of the first kind @f$ E(k,\phi) @f$
+   * for @c float modulus @f$ k @f$ and angle @f$ \phi @f$.
+   *
+   * @see ellint_1 for details.
+   */
   inline float
   ellint_1f(float __k, float __phi)
   { return __detail::__ellint_1<float>(__k, __phi); }
 
+  /**
+   * Return the incomplete elliptic integral of the first kind @f$ E(k,\phi) @f$
+   * for <tt>long double</tt> modulus @f$ k @f$ and angle @f$ \phi @f$.
+   *
+   * @see ellint_1 for details.
+   */
   inline long double
   ellint_1l(long double __k, long double __phi)
   { return __detail::__ellint_1<long double>(__k, __phi); }
 
+  /**
+   * Return the incomplete elliptic integral of the first kind @f$ F(k,\phi) @f$
+   * for @c real modulus @f$ k @f$ and angle @f$ \phi @f$.
+   *
+   * The incomplete elliptic integral of the first kind is defined as
+   * @f[
+   *   F(k,\phi) = \int_0^{\phi}\frac{d\theta}
+   * 				     {\sqrt{1 - k^2 sin^2\theta}}
+   * @f]
+   * For  @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
+   * the first kind, @f$ K(k) @f$.  @see comp_ellint_1.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @tparam _Tpp The floating-point type of the angle @c __phi.
+   * @param  __k  The modulus, <tt> abs(__k) <= 1 </tt>
+   * @param  __phi  The integral limit argument in radians
+   * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
+   */
   template<typename _Tp, typename _Tpp>
     inline typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type
     ellint_1(_Tp __k, _Tpp __phi)
@@ -275,14 +738,44 @@
 
   // Incomplete elliptic integrals of the second kind
 
+  /**
+   * @brief Return the incomplete elliptic integral of the second kind
+   * @f$ E(k,\phi) @f$ for @c float argument.
+   *
+   * @see ellint_2 for details.
+   */
   inline float
   ellint_2f(float __k, float __phi)
   { return __detail::__ellint_2<float>(__k, __phi); }
 
+  /**
+   * @brief Return the incomplete elliptic integral of the second kind
+   * @f$ E(k,\phi) @f$.
+   *
+   * @see ellint_2 for details.
+   */
   inline long double
   ellint_2l(long double __k, long double __phi)
   { return __detail::__ellint_2<long double>(__k, __phi); }
 
+  /**
+   * Return the incomplete elliptic integral of the second kind
+   * @f$ E(k,\phi) @f$.
+   *
+   * The incomplete elliptic integral of the second kind is defined as
+   * @f[
+   *   E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
+   * @f]
+   * For  @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
+   * the second kind, @f$ E(k) @f$.  @see comp_ellint_2.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @tparam _Tpp The floating-point type of the angle @c __phi.
+   * @param  __k  The modulus, <tt> abs(__k) <= 1 </tt>
+   * @param  __phi  The integral limit argument in radians
+   * @return  The elliptic function of the second kind.
+   * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
+   */
   template<typename _Tp, typename _Tpp>
     inline typename __gnu_cxx::__promote_2<_Tp, _Tpp>::__type
     ellint_2(_Tp __k, _Tpp __phi)
@@ -293,14 +786,49 @@
 
   // Incomplete elliptic integrals of the third kind
 
+  /**
+   * @brief Return the incomplete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu,\phi) @f$ for @c float argument.
+   *
+   * @see ellint_3 for details.
+   */
   inline float
   ellint_3f(float __k, float __nu, float __phi)
   { return __detail::__ellint_3<float>(__k, __nu, __phi); }
 
+  /**
+   * @brief Return the incomplete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu,\phi) @f$.
+   *
+   * @see ellint_3 for details.
+   */
   inline long double
   ellint_3l(long double __k, long double __nu, long double __phi)
   { return __detail::__ellint_3<long double>(__k, __nu, __phi); }
 
+  /**
+   * @brief Return the incomplete elliptic integral of the third kind
+   * @f$ \Pi(k,\nu,\phi) @f$.
+   *
+   * The incomplete elliptic integral of the third kind is defined by:
+   * @f[
+   *   \Pi(k,\nu,\phi) = \int_0^{\phi}
+   * 			 \frac{d\theta}
+   * 			 {(1 - \nu \sin^2\theta)
+   * 			  \sqrt{1 - k^2 \sin^2\theta}}
+   * @f]
+   * For  @f$ \phi= \pi/2 @f$ this becomes the complete elliptic integral of
+   * the third kind, @f$ \Pi(k,\nu) @f$.  @see comp_ellint_3.
+   *
+   * @tparam _Tp The floating-point type of the modulus @c __k.
+   * @tparam _Tpn The floating-point type of the argument @c __nu.
+   * @tparam _Tpp The floating-point type of the angle @c __phi.
+   * @param  __k  The modulus, <tt> abs(__k) <= 1 </tt>
+   * @param  __nu  The second argument
+   * @param  __phi  The integral limit argument in radians
+   * @return  The elliptic function of the third kind.
+   * @throw std::domain_error if <tt> abs(__k) > 1 </tt>.
+   */
   template<typename _Tp, typename _Tpn, typename _Tpp>
     inline typename __gnu_cxx::__promote_3<_Tp, _Tpn, _Tpp>::__type
     ellint_3(_Tp __k, _Tpn __nu, _Tpp __phi)
@@ -311,14 +839,36 @@
 
   // Exponential integrals
 
+  /**
+   * Return the exponential integral @f$ Ei(x) @f$ for @c float argument @c x.
+   *
+   * @see expint for details.
+   */
   inline float
   expintf(float __x)
   { return __detail::__expint<float>(__x); }
 
+  /**
+   * Return the exponential integral @f$ Ei(x) @f$
+   * for <tt>long double</tt> argument @c x.
+   *
+   * @see expint for details.
+   */
   inline long double
   expintl(long double __x)
   { return __detail::__expint<long double>(__x); }
 
+  /**
+   * Return the exponential integral @f$ Ei(x) @f$ for @c real argument @c x.
+   *
+   * The exponential integral is given by
+   * \f[
+   *   Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
+   * \f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __x  The argument of the exponential integral function.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     expint(_Tp __x)
@@ -329,14 +879,44 @@
 
   // Hermite polynomials
 
+  /**
+   * Return the Hermite polynomial @f$ H_n(x) @f$ of nonnegative order n
+   * and float argument @c x.
+   *
+   * @see hermite for details.
+   */
   inline float
   hermitef(unsigned int __n, float __x)
   { return __detail::__poly_hermite<float>(__n, __x); }
 
+  /**
+   * Return the Hermite polynomial @f$ H_n(x) @f$ of nonnegative order n
+   * and <tt>long double</tt> argument @c x.
+   *
+   * @see hermite for details.
+   */
   inline long double
   hermitel(unsigned int __n, long double __x)
   { return __detail::__poly_hermite<long double>(__n, __x); }
 
+  /**
+   * Return the Hermite polynomial @f$ H_n(x) @f$ of order n
+   * and @c real argument @c x.
+   *
+   * The Hermite polynomial is defined by:
+   * @f[
+   *   H_n(x) = (-1)^n e^{x^2} \frac{d^n}{dx^n} e^{-x^2}
+   * @f]
+   *
+   * The Hermite polynomial obeys a reflection formula:
+   * @f[
+   *   H_n(-x) = (-1)^n H_n(x)
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param __n The order
+   * @param __x The argument
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     hermite(unsigned int __n, _Tp __x)
@@ -347,14 +927,40 @@
 
   // Laguerre polynomials
 
+  /**
+   * Returns the Laguerre polynomial @f$ L_n(x) @f$ of nonnegative degree @c n
+   * and @c float argument  @f$ x >= 0 @f$.
+   *
+   * @see laguerre for more details.
+   */
   inline float
   laguerref(unsigned int __n, float __x)
   { return __detail::__laguerre<float>(__n, __x); }
 
+  /**
+   * Returns the Laguerre polynomial @f$ L_n(x) @f$ of nonnegative degree @c n
+   * and <tt>long double</tt> argument @f$ x >= 0 @f$.
+   *
+   * @see laguerre for more details.
+   */
   inline long double
   laguerrel(unsigned int __n, long double __x)
   { return __detail::__laguerre<long double>(__n, __x); }
 
+  /**
+   * Returns the Laguerre polynomial @f$ L_n(x) @f$
+   * of nonnegative degree @c n and real argument @f$ x >= 0 @f$.
+   *
+   * The Laguerre polynomial is defined by:
+   * @f[
+   * 	 L_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n} (x^ne^{-x})
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param __n The nonnegative order
+   * @param __x The argument <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     laguerre(unsigned int __n, _Tp __x)
@@ -365,32 +971,92 @@
 
   // Legendre polynomials
 
+  /**
+   * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
+   * degree @f$ l @f$ and @c float argument @f$ |x| <= 0 @f$.
+   *
+   * @see legendre for more details.
+   */
   inline float
-  legendref(unsigned int __n, float __x)
-  { return __detail::__poly_legendre_p<float>(__n, __x); }
+  legendref(unsigned int __l, float __x)
+  { return __detail::__poly_legendre_p<float>(__l, __x); }
 
+  /**
+   * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
+   * degree @f$ l @f$ and <tt>long double</tt> argument @f$ |x| <= 0 @f$.
+   *
+   * @see legendre for more details.
+   */
   inline long double
-  legendrel(unsigned int __n, long double __x)
-  { return __detail::__poly_legendre_p<long double>(__n, __x); }
+  legendrel(unsigned int __l, long double __x)
+  { return __detail::__poly_legendre_p<long double>(__l, __x); }
 
+  /**
+   * Return the Legendre polynomial @f$ P_l(x) @f$ of nonnegative
+   * degree @f$ l @f$ and real argument @f$ |x| <= 0 @f$.
+   *
+   * The Legendre function of order @f$ l @f$ and argument @f$ x @f$,
+   * @f$ P_l(x) @f$, is defined by:
+   * @f[
+   *   P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l}
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param __l The degree @f$ l >= 0 @f$
+   * @param __x The argument @c abs(__x) <= 1
+   * @throw std::domain_error if @c abs(__x) > 1
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
-    legendre(unsigned int __n, _Tp __x)
+    legendre(unsigned int __l, _Tp __x)
     {
       typedef typename __gnu_cxx::__promote<_Tp>::__type __type;
-      return __detail::__poly_legendre_p<__type>(__n, __x);
+      return __detail::__poly_legendre_p<__type>(__l, __x);
     }
 
   // Riemann zeta functions
 
+  /**
+   * Return the Riemann zeta function @f$ \zeta(s) @f$
+   * for @c float argument @f$ s @f$.
+   *
+   * @see riemann_zeta for more details.
+   */
   inline float
   riemann_zetaf(float __s)
   { return __detail::__riemann_zeta<float>(__s); }
 
+  /**
+   * Return the Riemann zeta function @f$ \zeta(s) @f$
+   * for <tt>long double</tt> argument @f$ s @f$.
+   *
+   * @see riemann_zeta for more details.
+   */
   inline long double
   riemann_zetal(long double __s)
   { return __detail::__riemann_zeta<long double>(__s); }
 
+  /**
+   * Return the Riemann zeta function @f$ \zeta(s) @f$
+   * for real argument @f$ s @f$.
+   *
+   * The Riemann zeta function is defined by:
+   * @f[
+   * 	\zeta(s) = \sum_{k=1}^{\infty} k^{-s} \hbox{ for } s > 1
+   * @f]
+   * and
+   * @f[
+   * 	\zeta(s) = \frac{1}{1-2^{1-s}}\sum_{k=1}^{\infty}(-1)^{k-1}k^{-s}
+   *              \hbox{ for } 0 <= s <= 1
+   * @f]
+   * For s < 1 use the reflection formula:
+   * @f[
+   * 	\zeta(s) = 2^s \pi^{s-1} \sin(\frac{\pi s}{2}) \Gamma(1-s) \zeta(1-s)
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __s.
+   * @param __s The argument <tt> s != 1 </tt>
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     riemann_zeta(_Tp __s)
@@ -401,14 +1067,40 @@
 
   // Spherical Bessel functions
 
+  /**
+   * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
+   * and @c float argument @f$ x >= 0 @f$.
+   *
+   * @see sph_bessel for more details.
+   */
   inline float
   sph_besself(unsigned int __n, float __x)
   { return __detail::__sph_bessel<float>(__n, __x); }
 
+  /**
+   * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
+   * and <tt>long double</tt> argument @f$ x >= 0 @f$.
+   *
+   * @see sph_bessel for more details.
+   */
   inline long double
   sph_bessell(unsigned int __n, long double __x)
   { return __detail::__sph_bessel<long double>(__n, __x); }
 
+  /**
+   * Return the spherical Bessel function @f$ j_n(x) @f$ of nonnegative order n
+   * and real argument @f$ x >= 0 @f$.
+   *
+   * The spherical Bessel function is defined by:
+   * @f[
+   *  j_n(x) = \left(\frac{\pi}{2x} \right) ^{1/2} J_{n+1/2}(x)
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __n  The integral order <tt> n >= 0 </tt>
+   * @param  __x  The real argument <tt> x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     sph_bessel(unsigned int __n, _Tp __x)
@@ -419,14 +1111,43 @@
 
   // Spherical associated Legendre functions
 
+  /**
+   * Return the spherical Legendre function of nonnegative integral
+   * degree @c l and order @c m and float angle @f$ \theta @f$ in radians.
+   *
+   * @see sph_legendre for details.
+   */
   inline float
   sph_legendref(unsigned int __l, unsigned int __m, float __theta)
   { return __detail::__sph_legendre<float>(__l, __m, __theta); }
 
+  /**
+   * Return the spherical Legendre function of nonnegative integral
+   * degree @c l and order @c m and <tt>long double</tt> angle @f$ \theta @f$
+   * in radians.
+   *
+   * @see sph_legendre for details.
+   */
   inline long double
   sph_legendrel(unsigned int __l, unsigned int __m, long double __theta)
   { return __detail::__sph_legendre<long double>(__l, __m, __theta); }
 
+  /**
+   * Return the spherical Legendre function of nonnegative integral
+   * degree @c l and order @c m and real angle @f$ \theta @f$ in radians.
+   *
+   * The spherical Legendre function is defined by
+   * @f[
+   *  Y_l^m(\theta,\phi) = (-1)^m[\frac{(2l+1)}{4\pi}
+   *                              \frac{(l-m)!}{(l+m)!}]
+   *                   P_l^m(\cos\theta) \exp^{im\phi}
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the angle @c __theta.
+   * @param __l The order <tt> __l >= 0 </tt>
+   * @param __m The degree <tt> __m >= 0 </tt> and <tt> __m <= __l </tt>
+   * @param __theta The radian polar angle argument
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta)
@@ -437,14 +1158,40 @@
 
   // Spherical Neumann functions
 
+  /**
+   * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
+   * and @c float argument @f$ x >= 0 @f$.
+   *
+   * @see sph_neumann for details.
+   */
   inline float
   sph_neumannf(unsigned int __n, float __x)
   { return __detail::__sph_neumann<float>(__n, __x); }
 
+  /**
+   * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
+   * and <tt>long double</tt> @f$ x >= 0 @f$.
+   *
+   * @see sph_neumann for details.
+   */
   inline long double
   sph_neumannl(unsigned int __n, long double __x)
   { return __detail::__sph_neumann<long double>(__n, __x); }
 
+  /**
+   * Return the spherical Neumann function of integral order @f$ n >= 0 @f$
+   * and real argument @f$ x >= 0 @f$.
+   *
+   * The spherical Neumann function is defined by
+   * @f[
+   *    n_n(x) = \left(\frac{\pi}{2x} \right) ^{1/2} N_{n+1/2}(x)
+   * @f]
+   *
+   * @tparam _Tp The floating-point type of the argument @c __x.
+   * @param  __n  The integral order <tt> n >= 0 @f$ </tt>
+   * @param  __x  The real argument <tt> __x >= 0 </tt>
+   * @throw std::domain_error if <tt> __x < 0 </tt>.
+   */
   template<typename _Tp>
     inline typename __gnu_cxx::__promote<_Tp>::__type
     sph_neumann(unsigned int __n, _Tp __x)
@@ -463,14 +1210,44 @@
 
   // Confluent hypergeometric functions
 
+  /**
+   * Return the confluent hypergeometric function of @c float
+   * numeratorial parameter @c a, denominatorial parameter @c c,
+   * and argument @c x.
+   *
+   * @see conf_hyperg for details.
+   */
   inline float
   conf_hypergf(float __a, float __c, float __x)
   { return std::__detail::__conf_hyperg<float>(__a, __c, __x); }
 
+  /**
+   * Return the confluent hypergeometric function of @c long double
+   * numeratorial parameter @c a, denominatorial parameter @c c,
+   * and argument @c x.
+   *
+   * @see conf_hyperg for details.
+   */
   inline long double
   conf_hypergl(long double __a, long double __c, long double __x)
   { return std::__detail::__conf_hyperg<long double>(__a, __c, __x); }
 
+  /**
+   * Return the confluent hypergeometric function @f$ {}_1F_1(a;c;x) @f$
+   * of real numeratorial parameter @c a, denominatorial parameter @c c,
+   * and argument @c x.
+   *
+   * The confluent hypergeometric function is defined by
+   * @f[
+   *    {}_1F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n x^n}{(c)_n n!}
+   * @f]
+   * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$,
+   * @f$ (x)_0 = 1 @f$
+   *
+   * @param __a The numeratorial parameter
+   * @param __c The denominatorial parameter
+   * @param __x The argument
+   */
   template<typename _Tpa, typename _Tpc, typename _Tp>
     inline typename __gnu_cxx::__promote_3<_Tpa, _Tpc, _Tp>::__type
     conf_hyperg(_Tpa __a, _Tpc __c, _Tp __x)
@@ -481,14 +1258,45 @@
 
   // Hypergeometric functions
 
+  /**
+   * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
+   * of @ float numeratorial parameters @c a and @c b,
+   * denominatorial parameter @c c, and argument @c x.
+   *
+   * @see hyperg for details.
+   */
   inline float
   hypergf(float __a, float __b, float __c, float __x)
   { return std::__detail::__hyperg<float>(__a, __b, __c, __x); }
 
+  /**
+   * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
+   * of @ long double numeratorial parameters @c a and @c b,
+   * denominatorial parameter @c c, and argument @c x.
+   *
+   * @see hyperg for details.
+   */
   inline long double
   hypergl(long double __a, long double __b, long double __c, long double __x)
   { return std::__detail::__hyperg<long double>(__a, __b, __c, __x); }
 
+  /**
+   * Return the hypergeometric function @f$ {}_2F_1(a,b;c;x) @f$
+   * of real numeratorial parameters @c a and @c b,
+   * denominatorial parameter @c c, and argument @c x.
+   *
+   * The hypergeometric function is defined by
+   * @f[
+   *    {}_2F_1(a;c;x) = \sum_{n=0}^{\infty} \frac{(a)_n (b)_n x^n}{(c)_n n!}
+   * @f]
+   * where the Pochhammer symbol is @f$ (x)_k = (x)(x+1)...(x+k-1) @f$,
+   * @f$ (x)_0 = 1 @f$
+   *
+   * @param __a The first numeratorial parameter
+   * @param __b The second numeratorial parameter
+   * @param __c The denominatorial parameter
+   * @param __x The argument
+   */
   template<typename _Tpa, typename _Tpb, typename _Tpc, typename _Tp>
     inline typename __gnu_cxx::__promote_4<_Tpa, _Tpb, _Tpc, _Tp>::__type
     hyperg(_Tpa __a, _Tpb __b, _Tpc __c, _Tp __x)