Update documents (#96)

* rename besselroots to approx_besselroots * add warning for compatibility of besselroots * use DocMeta.setdocmeta! instead of setup = * add test for besselroots * update benchmark.md to generate results automatically * fix typo in besselroots.jl * move some texts from README.md to docs/src/*.md * update ci.yml to trigger ci with release * update ci.yml for jldoctest * update ci.yml again for jldoctest * update autogeneration of benchmark * add more (a)repl block * bump version to v0.4.7
JuliaApproximation · Feb 16, 2021 · 7740d0a · 7740d0a · hyrodium · Feb 17, 2021
1 parent 770bdd6
commit 7740d0a
Show file tree

Hide file tree

Showing 20 changed files with 109 additions and 173 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -1,7 +1,12 @@
 name: CI
 on:
-  - push
-  - pull_request
+  pull_request:
+    branches:
+      - master
+  push:
+    branches:
+      - master
+    tags: '*'
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
@@ -56,8 +61,9 @@ jobs:
             Pkg.instantiate()'
       - run: |
           julia --project=docs -e '
-            using Documenter: doctest
+            using Documenter
             using FastGaussQuadrature
+            DocMeta.setdocmeta!(FastGaussQuadrature, :DocTestSetup, :(using LinearAlgebra, SpecialFunctions, FastGaussQuadrature))
             doctest(FastGaussQuadrature)'
       - run: julia --project=docs docs/make.jl
         env:

diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "FastGaussQuadrature"
 uuid = "442a2c76-b920-505d-bb47-c5924d526838"
-version = "0.4.6"
+version = "0.4.7"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

diff --git a/README.md b/README.md
@@ -17,93 +17,6 @@ For a quirky account on the history of computing Gauss-Legendre quadrature, see
 * The fastest Julia code for Gauss quadrature nodes and weights (without tabulation).
 * Change the perception that Gauss quadrature rules are expensive to compute.
 
-## Examples
-Here we compute `100000` nodes and weights of the Gauss rules.
-Try a million or ten million.
-
-```julia
-julia> @time gausschebyshev( 100000 );
-  0.001788 seconds (4 allocations: 1.526 MiB)
-
-julia> @time gausslegendre( 100000 );
-  0.002976 seconds (10 allocations: 2.289 MiB)
-
-julia> @time gaussjacobi( 100000, .9, -.1 );
-  0.894373 seconds (3.59 k allocations: 1.255 GiB, 36.38% gc time)
-
-julia> @time gaussradau( 100000 );
-  0.684122 seconds (3.59 k allocations: 1.256 GiB, 21.71% gc time)
-
-julia> @time gausslobatto( 100000 );
-  0.748166 seconds (3.57 k allocations: 1.256 GiB, 27.78% gc time)
-
-julia> @time gausslaguerre( 100000 );
-  0.156867 seconds (7 allocations: 2.292 MiB)
-
-julia> @time gausshermite( 100000 );
-  0.175055 seconds (386 allocations: 67.916 MiB, 9.18% gc time)
-```
-
-The paper [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969) computed a billion Gauss-Legendre nodes.
-So here we will do a billion + 1.
-```julia
-julia> @time gausslegendre( 1000000001 );
- 24.441304 seconds (10 allocations: 22.352 GiB, 2.08% gc time)
-```
-(The nodes near the endpoints coalesce in 16-digits of precision.)
-
-## The algorithm for Gauss-Chebyshev
-There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions:
-
-1. 1st kind, weight function `w(x) = 1/sqrt(1-x^2)`
-
-2. 2nd kind, weight function `w(x) = sqrt(1-x^2)`
-
-3. 3rd kind, weight function `w(x) = sqrt((1+x)/(1-x))`
-
-4. 4th kind, weight function `w(x) = sqrt((1-x)/(1+x))`
-
-They are all have explicit simple formulas for the nodes and weights [[4]](https://books.google.co.jp/books?id=8FHf0P3to0UC).
-
-## The algorithm for Gauss-Legendre
-Gauss quadrature for the weight function `w(x) = 1`.
-
-* For `n ≤ 5`: Use an analytic expression.
-* For `n ≤ 60`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate Legendre polynomials `Pₙ` and their derivatives `Pₙ'` by 3-term recurrence. Weights are related to `Pₙ'`.
-* For `n > 60`: Use asymptotic expansions for the Legendre nodes and weights [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969).
-
-## The algorithm for Gauss-Jacobi
-Gauss quadrature for the weight functions `w(x) = (1-x)^a(1+x)^b`, `a,b > -1`.
-
-* For `n ≤ 100`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate `Pₙ` and `Pₙ'` by three-term recurrence.
-* For `n > 100`: Use Newton's method to solve `Pₙ(x)=0`. Evaluate `Pₙ` and `Pₙ'` by an asymptotic expansion (in the interior of `[-1,1]`) and the three-term recurrence `O(n^-2)` close to the endpoints. (This is a small modification to the algorithm described in [[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).)
-* For `max(a,b) > 5`: Use the Golub-Welsch algorithm requiring `O(n^2)` operations.
-
-## The algorithm for Gauss-Radau
-Gauss quadrature for the weight function `w(x)=1`, except the endpoint `-1` is included as a quadrature node.
-
-The Gauss-Radau nodes and weights can be computed via the `(0,1)` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).
-
-## The algorithm for Gauss-Lobatto
-Gauss quadrature for the weight function `w(x)=1`, except the endpoints `-1` and `1` are included as nodes.
-
-The Gauss-Lobatto nodes and weights can be computed via the `(1,1)` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).
-
-## The algorithm for Gauss-Laguerre
-Gauss quadrature for the weight function `w(x) = exp(-x)` on `[0,Inf)`
-
-* For `n < 128`: Use the Golub-Welsch algorithm.
-* For `method=GLR`: Use the Glaser-Lui-Rohklin algorithm. Evaluate Laguerre polynomials `Lₙ` and their derivatives `Lₙ'` by using Taylor series expansions near roots generated by solving the second-order differential equation that `Lₙ` satisfies, see [[2]](http://epubs.siam.org/doi/pdf/10.1137/06067016X).
-* For `n ≥ 128`: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight `w(x) = x^alpha exp(-q_m x^m)` and this is O(sqrt(n)) when allowed to stop when the weights are below the smallest positive floating point number.
-
-## The algorithm for Gauss-Hermite
-Gauss quadrature for the weight function `w(x) = exp(-x^2)` on the real line.
-
-* For `n < 200`: Use Newton's method to solve `Hₙ(x)=0`. Evaluate Hermite polynomials `Hₙ` and their derivatives `Hₙ'` by three-term recurrence.
-* For `n ≥ 200`: Use Newton's method to solve `Hₙ(x)=0`. Evaluate `Hₙ` and `Hₙ'` by a uniform asymptotic expansion, see [[7]](http://arxiv.org/abs/1410.5286).
-
-The paper [[7]](http://arxiv.org/abs/1410.5286) also derives an `O(n)` algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form `exp(-V(x))`, where `V(x)` is a real polynomial.
-
 ## Example usage
 ```julia
 julia> @time nodes, weights = gausslegendre( 100000 );
@@ -114,20 +27,3 @@ julia> @time dot( weights, nodes.^2 )
   0.000184 seconds (7 allocations: 781.422 KiB)
 0.6666666666666665
 ```
-
-## References:
-[1] I. Bogaert, ["Iteration-free computation of Gauss-Legendre quadrature nodes and weights"](http://epubs.siam.org/doi/abs/10.1137/140954969), SIAM J. Sci. Comput., 36(3), A1008-A1026, 2014.
-
-[2] A. Glaser, X. Liu, and V. Rokhlin. ["A fast algorithm for the calculation of the roots of special functions."](http://epubs.siam.org/doi/pdf/10.1137/06067016X) SIAM J. Sci. Comput., 29 (2007), 1420-1438.
-
-[3] N. Hale and A. Townsend, ["Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights"](http://epubs.siam.org/doi/abs/10.1137/120889873), SIAM J. Sci. Comput., 2012.
-
-[4] J. C. Mason and D. C. Handscomb, ["Chebyshev Polynomials"](https://books.google.co.jp/books?id=8FHf0P3to0UC), CRC Press, 2002.
-
-[5] P. Opsomer, (in preparation).
-
-[6] A. Townsend, [The race for high order Gauss-Legendre quadrature](http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf), in SIAM News, March 2015.
-
-[7] A. Townsend, T. Trogdon, and S. Olver, ["Fast computation of Gauss quadrature nodes and weights on the whole real line"](http://arxiv.org/abs/1410.5286), to appear in IMA Numer. Anal., 2014.
-
-[8] M. Vanlessen, "Strong asymptotics of Laguerre-Type orthogonal polynomials and applications in Random Matrix Theory", Constr. Approx., 25:125-175, 2007.
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,6 +1,9 @@
 using Documenter
 using FastGaussQuadrature
 
+# Setup for doctests in docstrings
+DocMeta.setdocmeta!(FastGaussQuadrature, :DocTestSetup, :(using LinearAlgebra, SpecialFunctions, FastGaussQuadrature))
+
 makedocs(;
     modules = [FastGaussQuadrature],
     format = Documenter.HTML(

diff --git a/docs/src/benchmark.md b/docs/src/benchmark.md
@@ -3,33 +3,23 @@
 Here we compute `100000` nodes and weights of the Gauss rules.
 Try a million or ten million.
 
-```julia
-julia> @time gausschebyshev( 100000 );
-  0.001788 seconds (4 allocations: 1.526 MiB)
-
-julia> @time gausslegendre( 100000 );
-  0.002976 seconds (10 allocations: 2.289 MiB)
-
-julia> @time gaussjacobi( 100000, .9, -.1 );
-  0.894373 seconds (3.59 k allocations: 1.255 GiB, 36.38% gc time)
-
-julia> @time gaussradau( 100000 );
-  0.684122 seconds (3.59 k allocations: 1.256 GiB, 21.71% gc time)
-
-julia> @time gausslobatto( 100000 );
-  0.748166 seconds (3.57 k allocations: 1.256 GiB, 27.78% gc time)
-
-julia> @time gausslaguerre( 100000 );
-  0.156867 seconds (7 allocations: 2.292 MiB)
-
-julia> @time gausshermite( 100000 );
-  0.175055 seconds (386 allocations: 67.916 MiB, 9.18% gc time)
+```@repl
+using LinearAlgebra, BenchmarkTools, FastGaussQuadrature
+@btime gausschebyshev(100000);
+@btime gausslegendre(100000);
+@btime gaussjacobi(100000, 0.9, -0.1);
+@btime gaussradau(100000);
+@btime gausslobatto(100000);
+@btime gausslaguerre(100000);
+@btime gausshermite(100000);
 ```
 
 The paper [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969) computed a billion Gauss-Legendre nodes.
 So here we will do a billion + 1.
+
 ```julia
-julia> @time gausslegendre( 1000000001 );
- 24.441304 seconds (10 allocations: 22.352 GiB, 2.08% gc time)
+julia> @btime gausslegendre(1000000001);
+  23.363 s (10 allocations: 22.35 GiB)
 ```
+
 (The nodes near the endpoints coalesce in 16-digits of precision.)
diff --git a/docs/src/besselroots.md b/docs/src/besselroots.md
@@ -1,9 +1,9 @@
 # Roots of Bessel function
 
-Since [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl) doesn't have a method to calculate roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function), we implemented `besselroots`.
+Since [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl) doesn't have a method to calculate roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function), we implemented `approx_besselroots`.
 
 ```@docs
-besselroots(ν::Real, n::Integer)
+approx_besselroots(ν::Real, n::Integer)
 ```
 
-This method `besselroots` is used to calculate `gaussjacobi` and `gausslaguerre`.
+This method `approx_besselroots` is used to calculate `gaussjacobi` and `gausslaguerre`.
diff --git a/docs/src/gaussquadrature.md b/docs/src/gaussquadrature.md
@@ -2,20 +2,36 @@
 
 
 ## Gauss-Legendre quadrature
+Gauss quadrature for the weight function ``w(x) = 1``.
+
+* For ``n ≤ 5``: Use an analytic expression.
+* For ``n ≤ 60``: Use Newton's method to solve ``P_n(x)=0``. Evaluate Legendre polynomials ``P_n`` and their derivatives ``P_n'`` by 3-term recurrence. Weights are related to ``P_n'``.
+* For ``n > 60``: Use asymptotic expansions for the Legendre nodes and weights [[1]](http://epubs.siam.org/doi/abs/10.1137/140954969).
 
 ```@docs
 gausslegendre(n::Integer)
 ```
 
 
 ## Gauss-Hermite quadrature
+Gauss quadrature for the weight function ``w(x) = \exp(-x^2)`` on the real line.
+
+* For ``n < 200``: Use Newton's method to solve ``H_n(x)=0``. Evaluate Hermite polynomials ``H_n`` and their derivatives ``H_n'`` by three-term recurrence.
+* For ``n ≥ 200``: Use Newton's method to solve ``H_n(x)=0``. Evaluate ``H_n`` and ``H_n'`` by a uniform asymptotic expansion, see [[7]](http://arxiv.org/abs/1410.5286).
+
+The paper [[7]](http://arxiv.org/abs/1410.5286) also derives an ``O(n)`` algorithm for generalized Gauss-Hermite nodes and weights associated to weight functions of the form ``\exp(-V(x))``, where ``V(x)`` is a real polynomial.
 
 ```@docs
 gausshermite(n::Integer)
 ```
 
 
 ## Gauss-Laguerre quadrature
+Gauss quadrature for the weight function ``w(x) = \exp(-x)`` on ``[0,+\infty)``
+
+* For ``n < 128``: Use the Golub-Welsch algorithm.
+* For `method=GLR`: Use the Glaser-Lui-Rohklin algorithm. Evaluate Laguerre polynomials ``L_n`` and their derivatives ``L_n'`` by using Taylor series expansions near roots generated by solving the second-order differential equation that ``L_n`` satisfies, see [[2]](http://epubs.siam.org/doi/pdf/10.1137/06067016X).
+* For ``n ≥ 128``: Use a Newton procedure on Riemann-Hilbert asymptotics of Laguerre polynomials, see [5], based on [8]. There are some heuristics to decide which expression to use, it allows a general weight ``w(x) = x^\alpha \exp(-q_m x^m)`` and this is ``O(\sqrt{n})`` when allowed to stop when the weights are below the smallest positive floating point number.
 
 ```@docs
 gausslaguerre(n::Integer)
@@ -27,27 +43,46 @@ gausslaguerre(n::Integer, α::Real)
 
 
 ## Gauss-Chebyshev quadrature
+There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four weight functions:
+
+* 1st kind, weight function ``w(x) = 1/\sqrt{1-x^2}``
+* 2nd kind, weight function ``w(x) = \sqrt{1-x^2}``
+* 3rd kind, weight function ``w(x) = \sqrt{(1+x)/(1-x)}``
+* 4th kind, weight function ``w(x) = \sqrt{(1-x)/(1+x)}``
+
+They are all have explicit simple formulas for the nodes and weights [[4]](https://books.google.co.jp/books?id=8FHf0P3to0UC).
 
 ```@docs
 gausschebyshev(n::Integer, kind::Integer)
 ```
 
 
 ## Gauss-Jacobi quadrature
+Gauss quadrature for the weight functions ``w(x) = (1-x)^\alpha(1+x)^\beta``, ``\alpha,\beta > -1``.
+
+* For ``n ≤ 100``: Use Newton's method to solve ``P_n(x)=0``. Evaluate ``P_n`` and ``P_n'`` by three-term recurrence.
+* For ``n > 100``: Use Newton's method to solve ``P_n(x)=0``. Evaluate ``P_n`` and ``P_n'`` by an asymptotic expansion (in the interior of ``[-1,1]``) and the three-term recurrence ``O(n^{-2})`` close to the endpoints. (This is a small modification to the algorithm described in [[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).)
+* For ``max(a,b) > 5``: Use the Golub-Welsch algorithm requiring ``O(n^2)`` operations.
 
 ```@docs
 gaussjacobi(n::Integer, α::Real, β::Real)
 ```
 
 
 ## Gauss-Radau quadrature
+Gauss quadrature for the weight function ``w(x)=1``, except the endpoint ``-1`` is included as a quadrature node.
+
+The Gauss-Radau nodes and weights can be computed via the ``(0,1)`` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).
 
 ```@docs
 gaussradau(n::Integer)
 ```
 
 
 ## Gauss-Lobatto quadrature
+Gauss quadrature for the weight function ``w(x)=1``, except the endpoints ``-1`` and ``1`` are included as nodes.
+
+The Gauss-Lobatto nodes and weights can be computed via the ``(1,1)`` Gauss-Jacobi nodes and weights[[3]](http://epubs.siam.org/doi/abs/10.1137/120889873).
 
 ```@docs
 gausslobatto(n::Integer)

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -20,15 +20,11 @@ To check an integral
 \int_{-1}^{1} x^4 dx = \frac{2}{5}
 ```
 by numerically, try following code.
-```julia
-julia> using FastGaussQuadrature, LinearAlgebra
 
-julia> x, w = gausslegendre(3);
-
-julia> f(x) = x^4;
-
-julia> I = dot(w, f.(x));
-
-julia> I ≈ 2/5
-true
+```@repl
+using FastGaussQuadrature, LinearAlgebra
+x, w = gausslegendre(3)
+f(x) = x^4
+I = dot(w, f.(x))
+I ≈ 2/5
 ```
diff --git a/src/FastGaussQuadrature.jl b/src/FastGaussQuadrature.jl
@@ -11,6 +11,7 @@ export gausshermite
 export gaussjacobi
 export gausslobatto
 export gaussradau
+export approx_besselroots
 export besselroots
 
 import SpecialFunctions: besselj, airyai, airyaiprime

diff --git a/src/besselroots.jl b/src/besselroots.jl
@@ -1,30 +1,31 @@
 @doc raw"""
-    gausslegendre(ν::Real, n::Integer) -> Vector{Float64}
+    approx_besselroots(ν::Real, n::Integer) -> Vector{Float64}
 
 Return the first ``n`` roots of [Bessel function](https://en.wikipedia.org/wiki/Bessel_function).
+Note that this function is only 12-digits of precision.
 
 ```math
 J_{\nu}(x) = \sum_{m=0}^{\infty}\frac{(-1)^j}{\Gamma(\nu+j+1)j!} \left(\frac{x}{2}\right)^{2j+\nu}
 ```
 
 # Examples
-```jldoctest; setup = :(using FastGaussQuadrature, SpecialFunctions)
+```jldoctest
 julia> ν = 0.3;
 
-julia> roots = besselroots(ν, 10);
+julia> roots = approx_besselroots(ν, 10);
 
 julia> zeros = (x -> besselj(ν, x)).(roots);
 
 julia> all(zeros .< 1e-12)
 true
 ```
 """
-function besselroots(ν::Real, n::Integer)
+function approx_besselroots(ν::Real, n::Integer)
     # FIXME (related issue #22 and #80)
-    return besselroots(Float64(ν), n)
+    return approx_besselroots(Float64(ν), n)
 end
 
-function besselroots(ν::Float64, n::Integer)
+function approx_besselroots(ν::Float64, n::Integer)
 # DEVELOPERS NOTES:
 #   ν = 0 --> Full Float64 precision for n ≤ 20 (Wolfram Alpha), and very
 #     accurate approximations for n > 20 (McMahon's expansion)
@@ -219,3 +220,8 @@ function besselJ1(m)
     end
     return Jk2
 end
+
+function besselroots(ν::Real, n::Integer)
+    @warn "`besselroots` has been renamed to `approx_besselroots`, and `besselroots` will be removed in the next major release."
+    return approx_besselroots(ν, n)
+end
diff --git a/src/constants.jl b/src/constants.jl
@@ -71,8 +71,8 @@ const Piessens_C = @SMatrix [
 Values of Bessel function ``J_1`` on first ten roots of Bessel function `J_0`.
 
 # Examples
-```jldoctest; setup = :(using FastGaussQuadrature, SpecialFunctions)
-julia> roots = besselroots(0,10);
+```jldoctest
+julia> roots = approx_besselroots(0,10);
 
 julia> (besselj1.(roots)).^2 ≈ FastGaussQuadrature.besselJ1_10
 true