Add more docstrings, and generate document by Documenter.jl (#87)

* add docstrings and doctest for gausslegendre * fix jldoctest * fix jldoctest for gausslegendre * add docstring and jldoctest for gausshermite * add docstring and jldoctest for gausslaguerre * fix typos and form of README.md * add docstring and jldoctest for gausschebyshev * add docstring and jldoctest for gaussjacobi * add docstring and jldoctest for gaussradau * add docstring and jldoctest for gausslobatto * first commit for Documenter.jl * update gitignore for Documenter.jl * fix wording in docstring * update .travis.yml for Documenter * update authors in make.jl * update wording in index.md * update benchmark results * add pages for besselroots and benchmark * add badges for Documenter.jl * update docstring for besselroots * update travis badge (fix #82) * fix typos * update wording for clarity * update symbols to use unicode
JuliaApproximation · Jan 11, 2021 · bd8857c · bd8857c · dlfivefifty · Jan 11, 2021
1 parent dca47a7
commit bd8857c
Show file tree

Hide file tree

Showing 21 changed files with 872 additions and 421 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,2 +1,6 @@
 .DS_Store
 /Manifest.toml
+/dev/
+/docs/build/
+/docs/site/
+/docs/Manifest.toml
diff --git a/.travis.yml b/.travis.yml
@@ -7,9 +7,19 @@ julia:
   - 1.0
   - 1.5
   - nightly
-matrix:
-  allow_failures:
-    - julia: nightly  
 notifications:
   email: false
+jobs:
+  allow_failures:
+    - julia: nightly
+  fast_finish: true
+  include:
+    - stage: Documentation
+      julia: nightly
+      script: julia --project=docs -e '
+          using Pkg;
+          Pkg.develop(PackageSpec(path=pwd()));
+          Pkg.instantiate();
+          include("docs/make.jl");'
+      after_success: skip
 codecov: true
diff --git a/README.md b/README.md
@@ -1,47 +1,54 @@
 FastGaussQuadrature.jl
 =========
-[![Build Status](https://travis-ci.org/JuliaApproximation/FastGaussQuadrature.jl.svg?branch=master)](https://travis-ci.org/JuliaApproximation/FastGaussQuadrature.jl) [![codecov](https://codecov.io/gh/JuliaApproximation/FastGaussQuadrature.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastGaussQuadrature.jl)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastGaussQuadrature.jl/stable)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/FastGaussQuadrature.jl/dev)
+[![Build Status](https://travis-ci.com/JuliaApproximation/FastGaussQuadrature.jl.svg?branch=master)](https://travis-ci.com/JuliaApproximation/FastGaussQuadrature.jl)
+[![codecov](https://codecov.io/gh/JuliaApproximation/FastGaussQuadrature.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastGaussQuadrature.jl)
 
-A Julia package to compute `n`-point Gauss quadrature nodes and weights to 16-digit accuracy and in `O(n)` time. So far the package includes `gausschebyshev()`, `gausslegendre()`, `gaussjacobi()`, `gaussradau()`, `gausslobatto()`, `gausslaguerre()`, and `gausshermite()`. This package is heavily influenced by <a href="http://www.chebfun.org">Chebfun</a>.
+A Julia package to compute `n`-point Gauss quadrature nodes and weights to 16-digit accuracy and in `O(n)` time.
+So far the package includes `gausschebyshev()`, `gausslegendre()`, `gaussjacobi()`, `gaussradau()`, `gausslobatto()`, `gausslaguerre()`, and `gausshermite()`.
+This package is heavily influenced by [Chebfun](http://www.chebfun.org).
 
-An introduction to Gauss quadrature can be found <a href="http://en.wikipedia.org/wiki/Gaussian_quadrature">here</a>. For a quirky account on the history of computing Gauss-Legendre quadrature, see <a href="http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf">[6]</a>.
+An introduction to Gauss quadrature can be found [here](http://en.wikipedia.org/wiki/Gaussian_quadrature).
+For a quirky account on the history of computing Gauss-Legendre quadrature, see [[6]](http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf).
 
 ## Our Aims
 
 * 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.
+Here we compute `100000` nodes and weights of the Gauss rules.
+Try a million or ten million.
 
-```
-@time gausschebyshev( 100000 );
-0.002681 seconds (9 allocations: 1.526 MB, 228.45% gc time)
+```julia
+julia> @time gausschebyshev( 100000 );
+  0.001788 seconds (4 allocations: 1.526 MiB)
 
-@time gausslegendre( 100000 ); 
-0.007110 seconds (17 allocations: 2.671 MB)
+julia> @time gausslegendre( 100000 );
+  0.002976 seconds (10 allocations: 2.289 MiB)
 
-@time gaussjacobi( 100000, .9, -.1 );
-1.782347 seconds (20.84 k allocations: 1.611 GB, 22.89% gc time)
+julia> @time gaussjacobi( 100000, .9, -.1 );
+  0.894373 seconds (3.59 k allocations: 1.255 GiB, 36.38% gc time)
 
-@time gaussradau( 100000 );
-1.849520 seconds (741.84 k allocations: 1.625 GB, 22.59% gc time)
+julia> @time gaussradau( 100000 );
+  0.684122 seconds (3.59 k allocations: 1.256 GiB, 21.71% gc time)
 
-@time gausslobatto( 100000 );
-1.905083 seconds (819.73 k allocations: 1.626 GB, 23.45% gc time)
+julia> @time gausslobatto( 100000 );
+  0.748166 seconds (3.57 k allocations: 1.256 GiB, 27.78% gc time)
 
-@time gausslaguerre( 100000 )
-.891567 seconds (115.19 M allocations: 3.540 GB, 3.05% gc time)
+julia> @time gausslaguerre( 100000 );
+  0.156867 seconds (7 allocations: 2.292 MiB)
 
-@time gausshermite( 100000 );
-0.249756 seconds (201.22 k allocations: 131.643 MB, 4.92% gc time)
+julia> @time gausshermite( 100000 );
+  0.175055 seconds (386 allocations: 67.916 MiB, 9.18% gc time)
 ```
 
-The paper <a href="http://epubs.siam.org/doi/abs/10.1137/140954969">[1]</a> computed a billion Gauss-Legendre nodes. So here we will do a billion + 1.
-```
-@time gausslegendre( 1000000001 );
-131.392154 seconds (17 allocations: 26.077 GB, 1.17% 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.)
 
@@ -56,80 +63,71 @@ There are four kinds of Gauss-Chebyshev quadrature rules, corresponding to four
 
 4. 4th kind, weight function `w(x) = sqrt((1-x)/(1+x))`
 
-They are all have explicit simple formulas for the nodes and weights <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&pg=PA180#v=onepage&q&f=false">[4]</a>.
+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 `Pn(x)=0`. Evaluate `Pn` and `Pn'` by 3-term recurrence. Weights are related to `Pn'`.
-
-* For `n>60`: Use asymptotic expansions for the Legendre nodes and weights <a href="http://epubs.siam.org/doi/abs/10.1137/140954969">[1]</a>.  
+* 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 `Pn(x)=0`. Evaluate `Pn` and `Pn'` by three-term recurrence.
+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 `Pn(x)=0`. Evaluate `Pn` and `Pn'` 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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.)
-
-* For `max(a,b)>5`: Use the Golub-Welsch algorithm requiring `O(n^2)` operations. 
+* 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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.
+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 <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">[3]</a>.
+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 `Ln` and `Ln'` by using Taylor series expansions near roots generated by solving the second-order differential equation that `Ln` satisfies, see <a href="http://epubs.siam.org/doi/pdf/10.1137/06067016X">[2]</a>.
-
-* 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.
+* 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 `Hn(x)=0`. Evaluate `Hn` and `Hn'` by three-term recurrence.
+* 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).
 
-* For `n>=200`: Use Newton's method to solve `Hn(x)=0`. Evaluate `Hn` and `Hn'` by a uniform asymptotic expansion, see <a href="http://arxiv.org/abs/1410.5286">[7]</a>.
-*
-The paper <a href="http://arxiv.org/abs/1410.5286">[7]</a> 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.
+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
-
-
-```
-@time nodes, weights = gausslegendre( 100000 );
-0.007890 seconds (19 allocations: 2.671 MB)
+```julia
+julia> @time nodes, weights = gausslegendre( 100000 );
+  0.002192 seconds (10 allocations: 2.289 MiB)
 
 # integrates f(x) = x^2 from -1 to 1
-@time dot( weights, nodes.^2 )
-0.004264 seconds (7 allocations: 781.484 KB)
-0.666666666666666
+julia> @time dot( weights, nodes.^2 )
+  0.000184 seconds (7 allocations: 781.422 KiB)
+0.6666666666666665
 ```
 
 ## References:
-[1] I. Bogaert, <a href="http://epubs.siam.org/doi/abs/10.1137/140954969">"Iteration-free computation of Gauss-Legendre quadrature nodes and weights"</a>, SIAM J. Sci. Comput., 36(3), A1008-A1026, 2014.
-       
-[2] A. Glaser, X. Liu, and V. Rokhlin. <a href="http://epubs.siam.org/doi/pdf/10.1137/06067016X">"A fast algorithm for the calculation of the roots of special functions."</a> SIAM J. Sci. Comput., 29 (2007), 1420-1438.
+[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, <a href="http://epubs.siam.org/doi/abs/10.1137/120889873">"Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature
-       nodes and weights"</a>, SIAM J. Sci. Comput., 2012.
+[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, <a href="http://books.google.com/books?id=8FHf0P3to0UC&lpg=PP1&dq=Mason%20and%20Handscomb&pg=PP1#v=onepage&q=Mason%20and%20Handscomb&f=false">"Chebyshev Polynomials"</a>, CRC Press, 2002.
+[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, <a href="http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf"> The race for high order Gauss-Legendre quadrature</a>, in SIAM News, March 2015.  
+[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, <a href="http://arxiv.org/abs/1410.5286">"Fast computation of Gauss quadrature nodes and weights on the whole real line"</a>, to appear in IMA Numer. Anal., 2014.
+[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
@@ -0,0 +1,3 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
diff --git a/docs/make.jl b/docs/make.jl
@@ -0,0 +1,23 @@
+using Documenter
+using FastGaussQuadrature
+
+makedocs(;
+    modules = [FastGaussQuadrature],
+    format = Documenter.HTML(
+        canonical = "https://JuliaApproximation.github.io/FastGaussQuadrature.jl/stable/",
+        assets = ["assets/favicon.ico"],
+    ),
+    pages = [
+        "Home" => "index.md",
+        "Gaussian Quadrature" => "gaussquadrature.md",
+        "Benchmark" => "benchmark.md",
+        "Roots of Bessel function" => "besselroots.md",
+        "References" => "reference.md",
+    ],
+    repo = "https://github.com/JuliaApproximation/FastGaussQuadrature.jl/blob/{commit}{path}#L{line}",
+    sitename = "FastGaussQuadrature.jl",
+    authors = "Alex Townsend, Sheehan Olver, Peter Opsomer, and contributors.",
+    assets = String[],
+)
+
+deploydocs(; repo = "github.com/JuliaApproximation/FastGaussQuadrature.jl")
diff --git a/docs/src/benchmark.md b/docs/src/benchmark.md
@@ -0,0 +1,35 @@
+# Benchmark
+
+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.)
diff --git a/docs/src/besselroots.md b/docs/src/besselroots.md
@@ -0,0 +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`.
+
+```@docs
+besselroots(ν::Real, n::Integer)
+```
+
+This method `besselroots` is used to calculate `gaussjacobi` and `gausslaguerre`.
diff --git a/docs/src/gaussquadrature.md b/docs/src/gaussquadrature.md
@@ -0,0 +1,55 @@
+# Gaussian Quadrature
+
+
+## Gauss-Legendre quadrature
+
+```@docs
+gausslegendre(n::Integer)
+```
+
+
+## Gauss-Hermite quadrature
+
+```@docs
+gausshermite(n::Integer)
+```
+
+
+## Gauss-Laguerre quadrature
+
+```@docs
+gausslaguerre(n::Integer)
+```
+
+```@docs
+gausslaguerre(n::Integer, α::Real)
+```
+
+
+## Gauss-Chebyshev quadrature
+
+```@docs
+gausschebyshev(n::Integer, kind::Integer)
+```
+
+
+## Gauss-Jacobi quadrature
+
+```@docs
+gaussjacobi(n::Integer, α::Real, β::Real)
+```
+
+
+## Gauss-Radau quadrature
+
+```@docs
+gaussradau(n::Integer)
+```
+
+
+## Gauss-Lobatto quadrature
+
+```@docs
+gausslobatto(n::Integer)
+```
+
diff --git a/docs/src/index.md b/docs/src/index.md
@@ -0,0 +1,34 @@
+# FastGaussQuadrature.jl
+
+## Abstract
+
+FastGaussQuadrature.jl is a Julia package to compute `n`-point Gauss quadrature nodes and weights to 16-digit accuracy and in `O(n)` time.
+So far the package includes `gausschebyshev()`, `gausslegendre()`, `gaussjacobi()`, `gaussradau()`, `gausslobatto()`, `gausslaguerre()`, and `gausshermite()`.
+This package is heavily influenced by [Chebfun](http://www.chebfun.org).
+
+An introduction to Gauss quadrature can be found [here](http://en.wikipedia.org/wiki/Gaussian_quadrature).
+For a quirky account on the history of computing Gauss-Legendre quadrature, see [[6]](http://pi.math.cornell.edu/~ajt/papers/QuadratureEssay.pdf).
+
+## Our Aims
+
+* The fastest Julia code for Gauss quadrature nodes and weights (without tabulation).
+* Change the perception that Gauss quadrature rules are expensive to compute.
+
+## First example
+To check an integral
+```math
+\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
+```