Add citations to docstrings

Author: abhro

docs/src/functions_list.md

@@ -151,3 +151,7 @@ ellipe
 eta
 zeta
 ```
+
+## Citations
+```@bibliography
+```

src/beta_inc.jl

@@ -190,7 +190,7 @@ External links:
 See also: [`beta_inc`](@ref)
 
 # Implementation
-`BFRAC(A,B,X,Y,LAMBDA,EPS)` from Didonato and Morris (1982)
+`BFRAC(A,B,X,Y,LAMBDA,EPS)` from [Didonato and Morris (1992)](@cite didonato_1992)
 """
 function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, lambda::Float64, epps::Float64)
     @assert a > 1.0
@@ -260,7 +260,7 @@ External links:
 See also: [`beta_inc`](@ref)
 
 # Implementation
-`BASYM(A,B,LAMBDA,EPS)` from Didonato and Morris (1982)
+`BASYM(A,B,LAMBDA,EPS)` from [Didonato and Morris (1992)](@cite didonato_1992)
 """
 function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64)
     @assert a >= 15.0
@@ -443,7 +443,7 @@ External links:
 See also: [`beta_inc`](@ref)
 
 # Implementation
-`FPSER(A,B,X,EPS)` from Didonato and Morris (1982)
+`FPSER(A,B,X,EPS)` from [Didonato and Morris (1992)](@cite didonato_1992)
 """
 function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64)
     @assert b < epps*min(1.0, a)
@@ -489,7 +489,7 @@ External links:
 See also: [`beta_inc`](@ref)
 
 # Implementation
-`APSER(A,B,X,EPS)` from Didonato and Morris (1982)
+`APSER(A,B,X,EPS)` from [Didonato and Morris (1992)](@cite didonato_1992)
 """
 function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64)
     @assert a   <= epps*min(1.0, b)
@@ -534,7 +534,7 @@ External links:
 See also: [`beta_inc`](@ref)
 
 # Implementation
-`BPSER(A,B,X,EPS)` from Didonato and Morris (1982)
+`BPSER(A,B,X,EPS)` from [Didonato and Morris (1992)](@cite didonato_1992)
 """
 function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64)
     @assert b <= 1.0 || b*x <= 0.7

src/betanc.jl

@@ -149,8 +149,8 @@ Compute the CDF of the noncentral beta distribution given by
 ```math
 I_{x}(a,b; \lambda) = \sum_{j=0}^{\infty} q(\lambda/2,j) I_{x}(a+j,b;0)
 ```
-For ``\lambda < 54`` : algorithm suggested by Lenth(1987) in `ncbeta_tail(a,b,lambda,x)`.
-Else for ``\lambda \geq 54``: modification in Chattamvelli(1997) in
+For ``\lambda < 54``: algorithm suggested by [Lenth (1987)](@cite lenth_1987) in `ncbeta_tail(a,b,lambda,x)`.
+Else for ``\lambda \geq 54``: modification in [Chattamvelli (1997)](@cite chattamvelli_1997) in
 `ncbeta_poisson(a,b,lambda,x)` by using both forward and backward recurrences.
 """
 function ncbeta(a::Float64, b::Float64, lambda::Float64, x::Float64)

src/ellip.jl

@@ -27,19 +27,10 @@ See also: [`ellipe(m)`](@ref SpecialFunctions.ellipe).
        ``\alpha`` by ``k = \sin \alpha``.
 
 # Implementation
-Using piecewise approximation polynomial as given in
-> 'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions',
-> Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr,
-> DOI 10.1007/s10569-009-9228-z,
-> <https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf>
-
-For ``m<0``, followed by
-> Fukushima, Toshio. (2014).
-> 'Precise, compact, and fast computation of complete elliptic integrals by piecewise
-> minimax rational function approximation'.
-> Journal of Computational and Applied Mathematics. 282.
-> DOI 10.13140/2.1.1946.6245.,
-> <https://www.researchgate.net/publication/267330394>
+Using piecewise approximation polynomial as given in [fukushima_2009](@citet).
+
+For ``m < 0``, followed by [fukushima_2015](@citet).
+
 As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
 """
 ellipk(m::Real) = _ellipk(float(m))
@@ -207,19 +198,10 @@ See also: [`ellipk(m)`](@ref SpecialFunctions.ellipk).
        ``\alpha`` by ``k=\sin \alpha``.
 
 # Implementation
-Using piecewise approximation polynomial as given in
-> 'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions',
-> Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr,
-> DOI 10.1007/s10569-009-9228-z,
-> <https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf>
-
-For ``m<0``, followed by
-> Fukushima, Toshio. (2014).
-> 'Precise, compact, and fast computation of complete elliptic integrals by piecewise
-> minimax rational function approximation'.
-> Journal of Computational and Applied Mathematics. 282.
-> DOI 10.13140/2.1.1946.6245.,
-> <https://www.researchgate.net/publication/267330394>
+Using piecewise approximation polynomial as given in [fukushima_2015](@citet)
+
+For ``m<0``, followed by [fukushima_2015](@citet).
+
 As suggested in this paper, the domain is restricted to ``(-\infty,1]``.
 """
 ellipe(m::Real) = _ellipe(float(m))

src/erf.jl

@@ -244,13 +244,8 @@ External links:
 See also: [`erf(x)`](@ref erf).
 
 # Implementation
-Using the rational approximants tabulated in:
-> J. M. Blair, C. A. Edwards, and J. H. Johnson,
-> "Rational Chebyshev approximations for the inverse of the error function",
-> Math. Comp. 30, pp. 827--830 (1976).
-> <https://doi.org/10.1090/S0025-5718-1976-0421040-7>,
-> <http://www.jstor.org/stable/2005402>
-combined with Newton iterations for `BigFloat`.
+Using the rational approximants tabulated in [Blair (1976)](@cite blair_1976) combined with
+Newton iterations for `BigFloat`.
 """
 erfinv(x::Real) = _erfinv(float(x))
 
@@ -431,12 +426,7 @@ External links:
 See also: [`erfc(x)`](@ref erfc).
 
 # Implementation
-Using the rational approximants tabulated in:
-> J. M. Blair, C. A. Edwards, and J. H. Johnson,
-> "Rational Chebyshev approximations for the inverse of the error function",
-> Math. Comp. 30, pp. 827--830 (1976).
-> <https://doi.org/10.1090/S0025-5718-1976-0421040-7>,
-> <http://www.jstor.org/stable/2005402>
+Using the rational approximants tabulated in [Blair (1976)](@citet blair_1976)
 combined with Newton iterations for `BigFloat`.
 """
 erfcinv(x::Real) = _erfcinv(float(x))

src/gamma_inc.jl

@@ -504,7 +504,8 @@ and ``P(a,x)/x^a`` is given by:
 (1 - J) / \Gamma(a+1)
 ```
 resulting from term-by-term integration of `gamma_inc(a,x,ind)`.
-This is used when `a < 1` and `x < 1.1` - Refer Eqn (9) in the paper.
+This is used when `a < 1` and `x < 1.1` - Refer Eqn (9) in the
+[paper by DiDonato & Morris (1986)](@cite didonato_1986).
 
 See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
 """
@@ -553,8 +554,8 @@ T(a,\lambda) = \sum_{0}^{N} c_{k}(z)a^{-k}
 This is a higher accuracy approximation of Temme expansion, which deals with
 the region near `a ≈ x` with `a` large.
 Refer Appendix F in the paper for the extensive set of coefficients calculated
-using Brent's multiple precision arithmetic(set at 50 digits) in
-> Brent, R. P. A Fortran multiple-precision arithmetic package, ACM Trans. Math. Softw. 4(1978), 57-70, doi: 10.1145/355769.355775.
+using Brent's multiple precision arithmetic (set at 50 digits) in
+[Brent (1978)](@cite brent_1978).
 
 External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8)
 
@@ -700,7 +701,7 @@ end
 
 Compute using Finite Sums for ``Q(a,x)`` when `a >= 1 && 2a` is integer.
 Used when `a <= x <= x0` and `a = n/2`.
-Refer Eqn (14) in the paper.
+Refer Eqn (14) in the [paper by DiDonato and Morris (1986)](@cite didonato_1986).
 
 See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
 """
@@ -820,10 +821,8 @@ function gamma_inc_inv_alarge(a::Float64, minpq::Float64, pcase::Bool)
     fp = -sqrt(inv2π*a)*exp(-0.5*a*eta*eta)/gammax(a)
     return (x0, fp)
 end
-# Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato, Alfred H Morris.
-# Published in Journal: ACM Transactions on Mathematical Software (TOMS)
-# Volume 12 Issue 4, Dec. 1986 Pages 377-393
-# doi>10.1145/22721.23109
+# Reference: # DiDonato & Morris (1986), doi: 10.1145/22721.23109,
+# citation key: didonato_1986
 
 @doc raw"""
     gamma_inc(a,x,IND=0)
@@ -841,8 +840,10 @@ In terms of these, the lower incomplete gamma function is
 ``\gamma(a,x) = P(a,x) \Gamma(a)`` and the upper incomplete
 gamma function is ``\Gamma(a,x) = Q(a,x) \Gamma(a)``.
 
-`IND ∈ [0,1,2]` sets accuracy: `IND=0` means 14 significant digits accuracy,
-`IND=1` means 6 significant digit, and `IND=2` means only 3 digit accuracy.
+`IND ∈ [0,1,2]` sets accuracy:
+- `IND=0` means 14 significant digits accuracy,
+- `IND=1` means 6 significant digit, and
+- `IND=2` means only 3 digit accuracy.
 
 External links:
 [DLMF 8.2.4](https://dlmf.nist.gov/8.2.4),

src/sincosint.jl

@@ -245,12 +245,7 @@ External links:
 See also: [`cosint(x)`](@ref SpecialFunctions.cosint).
 
 # Implementation
-Using the rational approximants tabulated in:
-> A.J. MacLeod,
-> "Rational approximations, software and test methods for sine and cosine integrals",
-> Numer. Algor. 12, pp. 259--272 (1996).
-> <https://doi.org/10.1007/BF02142806>,
-> <https://link.springer.com/article/10.1007/BF02142806>.
+Using the rational approximants tabulated in [MacLeod (1996)](@cite macleod_1996).
 
 Note: the second zero of ``\text{Ci}(x)`` has a typo that is fixed:
 ``r_1 = 3.38418 0422\mathbf{8} 51186 42639 78511 46402`` in the article, but is in fact:
@@ -278,12 +273,7 @@ External links:
 See also: [`sinint(x)`](@ref SpecialFunctions.sinint).
 
 # Implementation
-Using the rational approximants tabulated in:
-> A.J. MacLeod,
-> "Rational approximations, software and test methods for sine and cosine integrals",
-> Numer. Algor. 12, pp. 259--272 (1996).
-> <https://doi.org/10.1007/BF02142806>,
-> <https://link.springer.com/article/10.1007/BF02142806>.
+Using the rational approximants tabulated in [MacLeod (1996)](@cite macleod_1996).
 
 Note: the second zero of ``\text{Ci}(x)`` has a typo that is fixed:
 ``r_1 = 3.38418 0422\mathbf{8} 51186 42639 78511 46402`` in the article, but is in fact: