diff --git a/Project.toml b/Project.toml
index e91ca6f..791b9dc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "2.1.1"
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 
 [compat]
 DataStructures = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17"
diff --git a/README.md b/README.md
index 14dfde9..d1f8baf 100644
--- a/README.md
+++ b/README.md
@@ -9,11 +9,13 @@ Documentation:
 
 This package provides support for one-dimensional numerical integration in Julia using adaptive
 Gauss-Kronrod quadrature.
-The code was originally part of Base Julia.
+The code was originally part of Base Julia.  It supports integration of arbitrary numeric types,
+including arbitrary precision (`BigFloat`), and even integration of arbitrary normed vector spaces
+(e.g. matrix-valued integrands).
 
 The package provides three functions: `quadgk`, `gauss`, and `kronrod`.
 `quadgk` performs the integration, `gauss` computes Gaussian quadrature points and weights for integrating
-over the interval [-1, 1], and `kronrod` computes Kronrod points, weights, and embedded Gaussian quadrature
+over the interval [a, b], and `kronrod` computes Kronrod points, weights, and embedded Gaussian quadrature
 weights for integrating over [-1, 1].   Typical usage looks like:
 ```jl
 using QuadGK
@@ -23,10 +25,36 @@ which computes the integral of exp(–x²) from x=0 to x=1 to a relative toleran
 
 For more information, see the [documentation](https://JuliaMath.github.io/QuadGK.jl/stable).
 
+## Gaussian quadrature and arbitrary weight functions
+
+If you are computing many similar integrals of smooth functions, you may not need an adaptive
+integration — with a little experimentation, you may be able to decide on an appropriate number
+`N` of integration points in advance, and re-use this for all of your integrals.    In this case
+you can use `x, w = gauss(N, a, b)` to find the quadrature points `x` and weights `w`, so that
+`sum(f.(x) .* w)` is an `N`-point approximation to `∫f(x)dx` from `a` to `b`.
+
+For computing many integrands of similar functions with *singularities* (and/or infinite intervals),
+`x, w = gauss(W, N, a, b)` function allows you to pass a *weight function* `W(x)` as the first argument,
+so that `sum(f.(x) .* w)` is an `N`-point approximation to `∫W(x)f(x)dx` from `a` to `b`.   In this way,
+you can put all of the singularities etcetera into `W` and precompute an accurate quadrature rule as
+long as the remaining `f(x)` terms are smooth.   For example,
+```jl
+using QuadGK
+x, w = gauss(x -> exp(-x) / sqrt(x), 10, 0, Inf, rtol=1e-10)
+```
+computes the points and weights for performing `∫exp(-x)f(x)/√x dx` integrals from `0` to `∞`, so that
+there is a `1/√x` singularity in the integrand at `x=0` and an implicit singularity from the `x=∞` endpoint.    For example, with `f(x) = sin(x)`, the exact answer is `0.5^0.25*sin(π/8)*√π ≈ 0.5703705559915792`.  Using the points and weights above with `sum(sin.(x) .* w)`, we obtain `0.5703705399484373`, which is correct to 7 digits using only 10 `f(x)` evaluations.  Obtaining similar
+accuracy for the same integral from `quadgk` requires nearly 300 function evaluations.   However, the
+`gauss` function itself computes many (`2N`) numerical integrals of your weight function (multiplied
+by polynomials), so this is only more efficient if your `f(x)` is very expensive or if you need
+to compute a large number of integrals with the same `W`.
+
+See the [`gauss` documentation](https://juliamath.github.io/QuadGK.jl/stable/#QuadGK.gauss) for more information.
+
 ## Similar packages
 
 The [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl) package provides
-non-adaptive Gaussian quadrature with a wider variety of weight functions — it is a good choice you need to go to very high orders N, e.g. to integrate rapidly oscillating functions, or use weight functions that incorporate some known singularity in your integrand.  QuadGK, on the other hand, keeps the order N of the quadrature rule fixed and improves accuracy by subdividing the integration domain, which can be better if fine resolution is required only in a part of your domain (e.g if your integrand has a sharp peak or singularity somewhere that is not known in advance).
+non-adaptive Gaussian quadrature variety of built-in weight functions — it is a good choice you need to go to very high orders N, e.g. to integrate rapidly oscillating functions, or use weight functions that incorporate some standard singularity in your integrand.  QuadGK, on the other hand, keeps the order N of the quadrature rule fixed and improves accuracy by subdividing the integration domain, which can be better if fine resolution is required only in a part of your domain (e.g if your integrand has a sharp peak or singularity somewhere that is not known in advance).
 
 For multidimensional integration, see the [HCubature.jl](https://github.com/stevengj/HCubature.jl), [Cubature.jl](https://github.com/stevengj/Cubature.jl), and
 [Cuba.jl](https://github.com/giordano/Cuba.jl) packages.
diff --git a/src/QuadGK.jl b/src/QuadGK.jl
index ec4764d..fd5a168 100644
--- a/src/QuadGK.jl
+++ b/src/QuadGK.jl
@@ -31,5 +31,6 @@ import Base.Order.Reverse
 include("gausskronrod.jl")
 include("evalrule.jl")
 include("adapt.jl")
+include("weightedgauss.jl")
 
 end # module QuadGK
diff --git a/src/gausskronrod.jl b/src/gausskronrod.jl
index c4f5b56..9d3c53a 100644
--- a/src/gausskronrod.jl
+++ b/src/gausskronrod.jl
@@ -82,10 +82,10 @@ function eigvec1(b,z::Number,m=length(b)+1)
 end
 
 """
-    gauss([T,] N)
+    gauss([T,] N, a=-1, b=1)
 
 Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to
-integrate functions on the interval `(-1, 1)`,  i.e. `sum(w .* f.(x))`
+integrate functions on the interval `(a, b)`,  i.e. `sum(w .* f.(x))`
 approximates the integral.  Uses the method described in Trefethen &
 Bau, Numerical Linear Algebra, to find the `N`-point Gaussian quadrature
 in O(`N`²) operations.
@@ -104,7 +104,17 @@ function gauss(::Type{T}, N::Integer) where T<:AbstractFloat
     return (x, w)
 end
 
-gauss(N::Integer) = gauss(Float64, N)
+gauss(N::Integer) = gauss(Float64, N) # integration on the standard interval (-1,1)
+
+# re-scaled to an arbitrary interval:
+gauss(N::Integer, a::Real, b::Real) = gauss(typeof(float(b-a)), N, a, b)
+function gauss(::Type{T}, N::Integer, a::Real, b::Real) where T<:AbstractFloat
+    x, w = gauss(T, N)
+    s = T(b-a)/2
+    x .= a .+ (x .+ 1) .* s
+    w .*= abs(s)
+    return (x, w)
+end
 
 """
     kronrod([T,] n)
diff --git a/src/weightedgauss.jl b/src/weightedgauss.jl
new file mode 100644
index 0000000..8b7c333
--- /dev/null
+++ b/src/weightedgauss.jl
@@ -0,0 +1,54 @@
+# Constructing Gaussian quadrature weights for an
+# arbitrary weight function and integration bounds.
+
+using Polynomials
+
+"""
+    gauss(W, N, a, b; rtol=sqrt(eps), quad=quadgk)
+
+Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to
+integrate functions on the interval `(a, b)` multiplied by the weight function
+`W(x)`.  That is, `sum(w .* f.(x))` approximates the integral `∫ W(x)f(x)dx`
+from `a` to `b`.
+
+This function performs `2N` numerical integrals of polynomials against `W(x)`
+using the integration function `quad` (defaults to `quadgk`) with relative tolerance `rtol`
+(which defaults to half of the precision `eps` of the endpoints).
+This is followed by an O(N²) calculations. So, using a large order `N` is expensive.
+
+If `W` has lots of singularities that make it hard to integrate numerically,
+you may need to decrease `rtol`.   You can also pass in a specialized quadrature routine
+via the `quad` keyword argument, which should accept arguments `quad(f,a,b,rtol=_,atol=_)`
+similar to `quadgk`.  (This is useful if your weight function has discontinuities, in which
+case you might want to break up the integration interval at the discontinuities.)
+"""
+function gauss(W, N, a::Real,b::Real; rtol::Real=sqrt(eps(typeof(float(b-a)))), quad=quadgk)
+    # Uses the Lanczos recurrence described in Trefethen & Bau,
+    # Numerical Linear Algebra, to find the `N`-point Gaussian quadrature
+    # using O(N) integrals and O(N²) operations.
+    α = zeros(N)
+    β = zeros(N+1)
+    T = typeof(float(b-a))
+    x = Poly(T[0, 1])
+    q₀ = Poly(T[0])
+    wint = first(quad(W, a, b, rtol=rtol))
+    atol = rtol*wint
+    q₁ = Poly(T[T(1)/sqrt(wint)])
+    for n = 1:N
+        v = x * q₁
+        q₁v = q₁ * v
+        α[n] = first(quad(x -> W(x) * q₁v(x), a, b, rtol=rtol, atol=atol))
+        n == N && break
+        v -= β[n]*q₀ + α[n]*q₁
+        v² = v*v
+        β[n+1] = sqrt(first(quad(x -> W(x) * v²(x), a, b, rtol=rtol, atol=atol)))
+        q₀ = q₁
+        q₁ = v / β[n+1]
+    end
+
+    # TODO: handle BigFloat etcetera — requires us to update eignewt() to
+    #       support nonzero diagonal entries.
+    E = eigen(SymTridiagonal(α, β[2:N]))
+
+    return (E.values, wint .* abs2.(E.vectors[1,:]))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d71731e..45ec69a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,4 +58,17 @@ end
     @test @inferred(quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)) isa Tuple{Float64,Float64}
     @test @inferred(quadgk(x -> 1, 0, 1im)) === (1.0im, 0.0)
     @test @inferred(quadgk(x -> sin(10x), 0,1))[1] ≈ (1 - cos(10))/10
+end
+
+@testset "gauss" begin
+    x,w = gauss(10, -1, 1)
+    x′,w′ = @inferred gauss(x->1, 10, -1, 1)
+    @test x ≈ x′ ≈ [-0.9739065285171717, -0.8650633666889845, -0.6794095682990244, -0.4333953941292472, -0.14887433898163124, 0.14887433898163124, 0.4333953941292472, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717]
+    @test w ≈ w′ ≈ [0.06667134430868811, 0.14945134915058064, 0.2190863625159821, 0.26926671930999635, 0.29552422471475276, 0.29552422471475276, 0.26926671930999635, 0.2190863625159821, 0.14945134915058064, 0.06667134430868811]
+
+    # Gauss–Hermite quadrature
+    xH,wH = @inferred gauss(x->exp(-x^2), 5, -Inf, Inf, rtol=1e-10)
+    @test xH ≈ [-2.020182870456085632929, -0.9585724646138185071128, 0, 0.9585724646138185071128, 2.020182870456085632929]
+    @test wH ≈ [0.01995324205904591320774, 0.3936193231522411598285, 0.9453087204829418812257, 0.393619323152241159829, 0.01995324205904591320774]
+    @test sum(xH.^4 .* wH) ≈ quadgk(x -> x^4 * exp(-x^2), -Inf, Inf)[1]
 end
\ No newline at end of file