Reference

Quadriceps

Positive-weight cubature rules in d ≥ 1 dimensions for two weight functions:

• GH (ghpos): the Gaussian weight, by default the standard normal density N(0, I_d);
• Le (lepos): the uniform weight, by default the uniform density on [0,1]^d.

A rule of degree p integrates every polynomial of total degree ≤ p exactly. All weights are strictly positive. The rules shipped with the package are the smallest ones known to its author; see the README and the documentation for where each one comes from.

The two exported functions mirror gausshermite and gausslegendre from FastGaussQuadrature.jl:

X, w = ghpos(3, 4)        # as exact as the 4×4×4 Gauss–Hermite grid (degree 7), with 27 nodes
X, w = ghpos(3; p = 7)    # the same rule, requested by degree
X, w = lepos(2, 5)        # degree 9 on the unit square: 17 nodes instead of 25

The second argument q is the number of nodes of the one-dimensional Gauss rule whose exactness is wanted, as in gausshermite(q); the degree is p = 2q - 1.

Unexported but public: Quadriceps.available, Quadriceps.extended, Quadriceps.nnodes, Quadriceps.ruleinfo, Quadriceps.exactness_error.

ghpos([T = Float64,] d, q; normalize = true, pragmatic = false) -> X, w
ghpos([T = Float64,] d; p, normalize = true, pragmatic = false) -> X, w

Positive-weight cubature rule for the Gaussian weight in d dimensions: the smallest one the package has.

The first method follows gausshermite(q; normalize) of FastGaussQuadrature.jl. There, q is the number of nodes of the one-dimensional Gauss rule, which is exact to degree 2q - 1. Here, ghpos(d, q) returns a rule with the same exactness in d dimensions, degree p = 2q - 1, that replaces the q^d-node product grid; for d = 1 it is gausshermite(q) itself. The second method takes the degree directly, as the keyword p; ghpos(d, q) is ghpos(d; p = 2q - 1).

X is an n × d matrix with one node per row, w the vector of n weights, all strictly positive, and

\[\sum_{i=1}^n w_i f(X_{i,1}, …, X_{i,d}) = \int f(x)\, ω(x)\, dx\]

for every polynomial f of total degree ≤ p (up to rounding; see Quadriceps.ruleinfo for the measured error of each stored rule).

Keyword arguments

• p: the degree of exactness, p ≥ 0 (second method only, required). Rules are stored at odd degrees; a request is served by the smallest stored rule of degree ≥ p, so an even p gets the rule for p + 1.
• normalize = true: the weight is the standard normal density ω(x) = (2π)^(-d/2) exp(-|x|²/2), so the weights sum to 1 and the rule computes E f(Z) for Z ~ N(0, I_d). With normalize = false the weight is ω(x) = exp(-|x|²) and the weights sum to π^(d/2). The keyword means what it means in gausshermite, but the default differs: gausshermite defaults to normalize = false, and only its normalize = true is a rule for the normal density.
• pragmatic = false: what to do when no stored rule covers the request. With false, throw an ArgumentError. With true, return the cheapest (fewest nodes) tensor product of lower-dimensional rules instead: stored rules and one-dimensional Gauss–Hermite rules, combined over the split of d that minimizes the number of nodes. Such a product is a valid positive-weight rule of the requested degree; it is just not small. With pragmatic = true the product is also returned in the rare case that it has strictly fewer nodes than the stored rule, so the result is always the cheapest the package can build.

Number type

The optional first argument T is the element type of X and w. Beyond Float64 the data depend on how many significant bits T holds (precision(T); a type without that method counts as wide):

• up to 113 bits — Float128 of Quadmath.jl, Double64 of DoubleFloats.jl, BigFloat under setprecision(BigFloat, 113) — are served from data/rules128.bin: every node and weight is the 80-digit rule correctly rounded to IEEE binary128 (113 bits, about 34 digits), and for Float128 the conversion is exact;
• more than 113 bits — BigFloat at its default 256 bits, Float64x4 of MultiFloats.jl, … — are served from data/rules80.bin: the 80-digit rules of the Zenodo deposit (record 10.5281/zenodo.22881864) themselves, held in a 40-byte binary format. The result carries 80 correct digits, whatever the precision of T.

Both files ship with the package; nothing is fetched from anywhere. Rounded on to Float64 either is the Float64 rule (to the last bit, except that a coordinate of size 1e-30 standing for zero can come out one unit in the last place off). One-dimensional Gauss factors are computed to the accuracy of the data they are combined with. Quadriceps.extended lists the cells with the measured error of each rule in both formats: every stored rule is one of them, so T never fails on a cell the Float64 call answers.

using Quadmath
X, w = ghpos(Float128, 3, 4)             # the 27-node rule in quadruple precision
X, w = ghpos(BigFloat, 3, 4)             # the same rule to 80 digits
Float64.(X) ≈ ghpos(3, 4)[1]             # true: the same rule

Details

• d = 1 returns the Gauss–Hermite rule with q (or p ÷ 2 + 1) nodes, as an n × 1 matrix; it is never an error.
• The result is a fresh copy; mutating it does not affect later calls. Loaded rules are cached, so repeated calls are cheap.

Examples

X, w = ghpos(3, 4)                       # degree 7: 27 nodes instead of 4³ = 64
sum(w .* X[:, 1] .^ 2 .* X[:, 2] .^ 4)   # E[Z₁² Z₂⁴] = 3
ghpos(3; p = 7) == ghpos(3, 4)           # true

ghpos(3, 21)                             # ArgumentError: nothing stored at that degree
X, w = ghpos(3, 21; pragmatic = true)    # tensor product of lower-dimensional rules
X, w = ghpos(7, 5; pragmatic = true)     # d = 7 as (d = 2) × (d = 5)

See also lepos, Quadriceps.nnodes, Quadriceps.available.

lepos([T = Float64,] d, q; normalize = true, pragmatic = false) -> X, w
lepos([T = Float64,] d; p, normalize = true, pragmatic = false) -> X, w

Positive-weight cubature rule for the uniform weight on a d-dimensional cube: the smallest one the package has.

The first method follows gausslegendre(q) of FastGaussQuadrature.jl: q is the number of nodes of the one-dimensional Gauss rule, and lepos(d, q) returns a rule with the same exactness in d dimensions, degree p = 2q - 1, that replaces the q^d-node product grid; for d = 1 it is gausslegendre(q) itself. The second method takes the degree directly, as the keyword p; lepos(d, q) is lepos(d; p = 2q - 1).

X is an n × d matrix with one node per row, w the vector of n weights, all strictly positive, and sum(w[i] * f(X[i, :])) equals the integral of f against the weight for every polynomial f of total degree ≤ p (up to rounding; see Quadriceps.ruleinfo).

Keyword arguments

• p: the degree of exactness, p ≥ 0 (second method only, required); an even p gets the rule for p + 1.
• normalize = true: the weight is the uniform density on [0,1]^d, so the weights sum to 1 and the rule computes E f(U) for U uniform on the unit cube. With normalize = false the rule is for ∫ f(x) dx over [-1,1]^d, the convention of gausslegendre, and the weights sum to 2^d.
• pragmatic = false: as for ghpos. With false, a request that no stored rule covers is an ArgumentError; with true, the cheapest tensor product of lower-dimensional rules (stored rules and one-dimensional Gauss–Legendre rules) is returned instead.

The optional number type T and the details listed under ghpos apply here too. All nodes of every stored Le rule lie inside the cube.

Examples

X, w = lepos(2, 5)                       # degree 9: 17 nodes on [0,1]² instead of 25
sum(w .* X[:, 1] .^ 3 .* X[:, 2] .^ 2)   # 1/4 · 1/3

X, w = lepos(2, 5; normalize = false)    # the same rule on [-1,1]², weights sum to 4
X, w = lepos(6, 6; pragmatic = true)     # d = 6 from a product of stored rules

See also ghpos, Quadriceps.nnodes, Quadriceps.available.

Quadriceps.available(family)

All stored rules of family (:gh or :le), as a vector of Quadriceps.RuleInfo sorted by dimension and degree. Tensor products, which pragmatic = true builds on demand, are not listed.

[(r.d, r.p, r.n) for r in Quadriceps.available(:gh) if r.d == 3]

Quadriceps.extended(family)

The cells of family (:gh or :le) that are stored beyond double precision — every stored cell — and so available as ghpos(T, d, q) or lepos(T, d, q) for a number type T wider than Float64: a vector of named tuples (family, d, p, n, relerr128, relerr80), sorted by dimension and degree. Both errors are the largest relative monomial error of the rule, measured in much wider arithmetic: relerr128 of the rule rounded to IEEE binary128 (Float128; the unit is the machine epsilon of that format, 2^-112 ≈ 1.93e-34, and a few epsilons is the floor of the format, as a few 1e-16 is for the relerr of the Float64 rules), relerr80 of the 80-digit rule that a type wider than 113 bits receives, as the deposit measured it.

Quadriceps.nnodes(family, d, q; pragmatic = false) -> Integer
Quadriceps.nnodes(family, d; p, pragmatic = false) -> Integer

Number of nodes of the rule that ghpos (family = :gh) or lepos (family = :le) returns for the same arguments, without building it. Throws the same ArgumentError when there is no rule. The count is an Int when it fits and a BigInt otherwise (tensor products in high dimensions).

Quadriceps.ruleinfo(family, d, q; pragmatic = false)
Quadriceps.ruleinfo(family, d; p, pragmatic = false)

Describe the rule that ghpos or lepos returns for the same arguments: a vector with one entry per tensor factor, in the order of the columns of X. An entry is either the Quadriceps.RuleInfo of a stored rule (dimension, degree, node count, measured error, origin) or, for a one-dimensional Gauss factor, the named tuple (family, d = 1, p, n, origin = "Gauss"). A request answered by a single stored rule gives a one-element vector.

Use it to find out whom to cite: the origin field names the published source of every rule that is not the package author's own.

Quadriceps.RuleInfo

Catalog entry of one stored rule. Fields:

• family — :gh or :le;
• d, p, n — dimension, degree of exactness, number of nodes;
• moller — Möller's lower bound on n for this (d, p), or -1 where it is not tabulated; a rule with n == moller is proven minimal;
• relerr — largest relative monomial error of the stored Float64 rule over all monomials of total degree ≤ p, measured when the package data were built (see Quadriceps.exactness_error);
• minweight — smallest weight (always > 0);
• interior — true when every node lies inside the integration domain (always true for GH);
• origin — who the rule belongs to: "own", or "derived: …", "same-rule: …", "transcribed: …" followed by the published source;
• source_id — the same information as the small integer stored in data/rules.bin (0: own; 3: derived from Diallo and Worku; 10 and up: a published rule; see the format page).

Quadriceps.exactness_error(X, w, p, family) -> err

Largest relative monomial error of the rule (X, w) over all monomials of total degree ≤ p, for the normalized weight of family (:gh: N(0, I_d); :le: uniform on [0,1]^d):

\[\max_{|a| ≤ p} \frac{|\sum_i w_i x_i^a - \mathrm{E}\,x^a|}{\max(\sum_i |w_i|\,|x_i^a|,\ 1)} .\]

The denominator is the scale of the sum being computed, so a value near eps() means the rule is exact to rounding. The element type of X and w sets the arithmetic: pass BigFloat arrays to measure beyond double precision. Rules returned with normalize = false must be checked in the normalized frame.

X, w = ghpos(4, 5)                           # degree 2·5 - 1 = 9
Quadriceps.exactness_error(X, w, 9, :gh)     # ≈ 4e-16