Numerical computing tools and functions in Swift.

Fills in a few missing basic functions.

• expm1mx(_:) - Calculate e^x - 1 - x with precision even when x is small

• log1pmx(_:) - Calculate log(1 + x) - x with precision even when x is small

• xmsin(_:) - Calculate x - sin(x) with precision even when x is small

A collection of functions with many applications in stats/ml and the sciences.

• Error Function

  • invErf(_:) - Inverse of the Error function (erf)

  • invErfC(_:) - Inverse of the complement of the Error function (erfc)

• Gamma Function

  • gammaReg(a:x:) - The regularized incomplete gamma function. This return a Probability value that can store values close to 0 or 1 with precision.

  • pGamma(a:x:) - The regularized lower incomplete gamma function Pₐ(x) = γ(a,x) / Γ(a). This is the lower tail of gammaReg(a:x:).

  • qGamma(a:x:) - The regularized upper incomplete gamma function Qₐ(x) = Γ(a,x) / Γ(a). This is the upper tail of gammaReg(a:x:).

  • pGammaDeriv(a:x:) - The derivative of P(a,x), P'(a,x)

  • invGammaReg(a:pq:) - The inverse of the regularized incomplete gamma function such that gamma_reg(a,x) = p, where p is a Probability value.

  • invPGamma(a:p:) - The inverse of P such that P(a,x) = p

  • invQGamma(a:q:) - The inverse of Q such that Q(a,x) = q

• Beta Function

  • beta(a:b:) - The beta function, Beta(a,b) = Gamma(a) * Gamma(b) / Gamma(a + b)

  • lbeta(a:b:) - The log of the beta function

  • betaReg(x:a:b) - The regularized incomplete beta function, I_x(a,b). This return a Probability value that can store values close to 0 or 1 with precision.

  • betaRegDeriv(x:a:b:) - The derivative of I, Iʹ_x(a,b)

  • invBetaReg(p:a:b:) - The inverse of I such that I_x(a,b) = p, where p is a Probability value.

• Marcum Q Function

  • marcum(mu:x:y:) - The Marcum Q function, Qᵤ(x,y) = e^(-x) Σᵢ₌₀... xⁱ / i! Qᵤ₊ᵢ(y). Output is a Probability value that can store values close to 0 or 1 with precision.

  • marcumDeriv(mu:x:y:) - The derivative of the Marcum Q function with respect to y.

  • invMarcum(mu:x:p:) - The inverse of the Marcum Q function with respect to y such that Qᵤ(x,y) = p, where p is a Probability value.

Functions to find the root of a function of interest.

• root(guess:xmin?:xmax?:tolerance:method:f) - Finds a root of f without any derivatives. It first brackets the root and then finds it within that interval. Multiple methods to find the root are available, including: Secant, Brent's (default), Dekker's, Ridders', and TOMS Algo 748.

• root(guess:xmin?:xmax?:tolerance:f:f1) - Finds a root of f using its first derivative. This is the Newton-Raphson method.

• root(guess:xmin?:xmax?:tolerance:f:f1:f2) - Finds a root of f using its first and second derivatives. This uses Halley's method.

• root(guess:xmin?:xmax?:tolerance:f:f1:f2f1) - Same as above but takes the ratio of the second derivative to the first. For use in cases where the ratio is cheaper to calculate.

Numerically integrate a function on a closed interval.

• integrate(range:maxIter:method:f:) - Integrates the function f on the specified range using the specified method.

  • trapezoidal - Trapezoidal Rule: a simple method that is is especially efficient when integrating a periodic function over its period or when integrating a peak function. In other cases it will usually not be the most efficient method.

  • romberg - (default) Romberg's Method: this method layers Richardson extrapolation on top of the Trapezoidal Rule to achieve more accuracy for the same number of function evaluations in many cases.

• More flexible stopping criteria for sequences - extensions on Sequence that can prove powerful for defining converence and other stopping criteria for numerical methods

  • first(where:) - a variant of the Standard Library function where the closure also has access to the preceding element.

  • until(minIter:maxIter:_) - finds the first element after minIter and before maxIter satisfying the predicate, or the last element it processed otherwise. It returns an IterativeResult value that tells how many elements it saw and how it exited. Available with the closure seeing only the last element or also seeing the preceding element.

• Approximate Equality - methods for determing if two FloatingPoint numbers are close to each other in a way useful for numerical computations. This introduces an EqualityTarget type which let's you specify if the number you are comparing to is known to be zero or only might be zero, and an EqualityTolerance type which lets you specify tolerance both relatively and absolutely. If you know you are comparing to zero there is an option to provide a value relative to which tolerance will be calculated. Several preset tolerances are provided as static members to represent best practices of different strictness.

  • isApprox(_:tolerance:) - extension on FloatingPoint that takes an EqualityTarget and an EqualityTolerance. Default tolerance is matching about half the digits supported by your type.

  • isWithinULP(of:n:) - extension on BinaryFloatingPoint that lets you specify tolerance as the number of discrete floating point values away another number is. This one has no special handling for zero and should not be used there.

• Summation - Provides some standard summation methods

  • sumNaive() - Simple sequential summation. This method is subject to accumulation of rounding and truncation errors. Implemented as extension on Sequence where Element is FloatingPoint.

  • sumPairwise() - Pairwise summation. As fast as naive sum but accumulates error much more slowly. Implemented as extension on Collection where Element is FloatingPoint and Index is Int.

  • sumKahan() - Kahan's compensated sum. Estimates the error after adding each term and tries to correct for it. More accurate in many cases than pairwise but slower. Implemented as extension on Sequence where Element is FloatingPoint.

  • sumKBN() - Kahan-Babuška-Neumaier sum. Improves on Kahan by compensating for rounding error in either the sum or the addend. The same number of operations as Kahan sum and more accurate so KBN should be preferred. Implemented as extension on Sequence where Element is FloatingPoint.

• Series - Some tools around evaluating convergent series and products.

  • series(indices:initialSum:initialState:maxIter:tolerance:update:) - Evaluates the potentially infinite series truncated either at its maxIter term or when it converges within tolerance. Terms are defined by the closure update which has access to the index and optionally some State as of the previous term.

  • product(indices:initialProduct:initialState:maxIter:tolerance:update:) - Evaluates the potentially infinite product truncated either at its maxIter term or when it converges within tolerance. Terms are defined by the closure update which has access to the index and optionally some State as of the previous term.

• Polynomials - Evaluation of polynomials using Horner's method. Also includes ratios of polynomials.

  • polynomial(coeffs:z:) - evaluate the polynomial with coefficients coeffs at point z

  • polynomialRatio(num:denom:z:) - evaluate the ratio of polynomials defined by num and denom at point z

• Chebyshev polynomials - Evaluation of Chebyshev polynomials, including rescaling the interval.

  • chebyshev(coeffs:z:interval:maxTerms:) - evaluates the Chebyshev polynomial defined by coeffs at point z. Optionally an interval other than [-1,1] can be specified and/or a maximum number of terms.

• Continued Fractions - Evaluation of continued fractions using the modified Lentz method.

  • continuedFraction(b0:a:b:maxIter:) - evaluate a continued fraction of the form b₀ + a₁ / (b₁ +) a₂ / (b₂ +) ... The terms aᵢ are defined by the closure a and the terms bᵢ by the closure b for i>0. The initial term b₀ is specified by b0. A maximum number of levels may be set by maxIter.

Extensive accuracy measurement and testing to make sure accuracy doesn't change from the expected levels. Accuracy is measured in terms of Log Relative Error.

Easy is 10k random doubles from [0,1] and their negatives shuffled into an array. Hard is 10k random doubles drawn from [0,1] multiplied by random powers of ten between -10 and 10. Peters is the small pathological case [1.0,1e100,1.0,-1e100].

(number of function evaluations in parentheses. for bracketing methods first number in parentheses is how many function iterations to bracket. * indicates method didn't converge. † indicates converged to a different root.)

(number of terms in parentheses. * indicates method didn't converge.)

(number of terms in parentheses. * indicates method didn't converge.)

(number of terms in parentheses. * indicates method didn't converge.)

(number of function evaluations in parentheses. * indicates method didn't converge.)