MathUtils Module Source file:MathUtils

Dependencies

precision Defined in constants
MpiUtils
constants Defined in constants
MiscUtils
config Defined in config
iso_fortran_env

Subroutines

airy_fast(x, ai, aip)

Fast real Airy Ai(x) and Ai'(x) Optimized for < 5e-8 absolute error on Ai and < 7e-8 on Ai'. The Ai branch cuts and polynomial coefficients are shared verbatim with airy_ai_fast; keep the two in step when refitting. Relative Airy errors can be large near Airy zeros, so absolute error is the intended accuracy diagnostic for this routine. Branch joins are endpoint-matched for smoothness. Measured Ai jumps at the hot-interval joins are below 1.0e-12, and Ai' jumps are below 2.0e-14 in double precision. Outside the hot interval the fallbacks are deliberately simplified for speed at the existing +/-6.5 cut points: the positive fallback keeps one correction term in the decaying asymptotic form, while the negative fallback uses first-degree polynomial corrections to the oscillatory asymptotic form. They are not intended to be moved much closer to the origin without refitting. Asymptotic-branch accuracy, measured against scipy.special.airy on [-50,-6.5) and (6.5,25.77], has max absolute Ai errors about 2.3e-8 on the negative side and 7.9e-10 on the positive side, and max absolute Ai' errors about 5.9e-8 and 2.4e-9. The positive decaying branch has worst relative errors about 2.9e-4 for Ai and 3.3e-4 for Ai' before the x > 25.77 zero cutoff; beyond the cutoff relative error is 100%, but absolute error is below 1e-38 at entry.
- real(dl) intent(in) :: x
- real(dl) intent(out) :: ai
- real(dl) intent(out) :: aip
AiryAiFastArray(n, ai, x)
- integer intent(in) :: n
- real(dl) intent(out) :: ai(n)
- real(dl) intent(in) :: x(n)
AiryFastArray(n, ai, aip, x)
- integer intent(in) :: n
- real(dl) intent(out) :: ai(n)
- real(dl) intent(out) :: aip(n)
- real(dl) intent(in) :: x(n)
brentq(obj, f, ax, bx, tol, xzero, fzero, iflag, fax, fbx)

> Find a zero of the function \( f(x) \) in the given interval \( [a_x,b_x] \) to within a tolerance \( 4 \epsilon |x| + tol \), where \( \epsilon \) is the relative machine precision defined as the smallest representable number such that \( 1.0 + \epsilon > 1.0 \). It is assumed that \( f(a_x) \) and \( f(b_x) \) have opposite signs. #References R. P. Brent, "[An algorithm with guaranteed convergence for finding a zero of a function](http://maths-people.anu.edu.au/~brent/pd/rpb005.pdf)", The Computer Journal, Vol 14, No. 4., 1971. R. P. Brent, "Algorithms for minimization without derivatives", http://maths-people.anu.edu.au/~brent/pub/pub011.html Prentice-Hall, Inc., 1973. # See also 1. [zeroin.f](http://www.netlib.org/go/zeroin.f) from Netlib
- class(*) :: obj
- procedure(obj_function) :: f f(obj, x), unwrapping obj with select type
- real(dl) intent(in) :: ax left endpoint of initial interval
- real(dl) intent(in) :: bx right endpoint of initial interval
- real(dl) intent(in) :: tol desired length of the interval of uncertainty
- real(dl) intent(out) :: xzero abscissa approximating a zero of `f` in the interval `ax`,`bx`
- real(dl) intent(out) :: fzero value of `f` at the root (`f(xzero)`)
- integer intent(out) :: iflag status flag (`-1` = error, `0` = root found)
- real(dl) intent(in), optional :: fax if `f(ax)` is already known, it can be input here
- real(dl) intent(in), optional :: fbx if `f(bx)` is already known, it can be input here
Gauss_Legendre(n, x, w)

Get n-point Gauss-Legendre nodes x and weights w on [-1, 1].
- integer intent(in) :: n
- real(dl) intent(out) :: x(n)
- real(dl) intent(out) :: w(n)
GetThreeJs(l2in, l3in, m2in, m3in, thrcof)

Recursive evaluation of 3j symbols. Does minimal error checking on input parameters. Generates the set of 3j-coeffs (l1,l2,l3\\ m1,m2,m3) for all allowed l1 from l1min = max(abs(l2-l3),abs(m1)) to l1max = l2 + l3, with m1 = -(m2 + m3). The resulting 3j-coeffs are stored as thrcof(l1-l1min+1). For numerical stability the recursion proceeds simultaneously forwards from l1min and backwards from l1max, with the two branches matched over three overlapping points where the forward recursion stops being stable.
- integer intent(in) :: l2in
- integer intent(in) :: l3in
- integer intent(in) :: m2in
- integer intent(in) :: m3in
- real(dl) intent(out) :: thrcof(*)
integrate_3j(lmax, lmax_w, n, W, dopol, M)

$ use omp_lib, only: omp_get_thread_num, omp_get_max_threads Get coupling matrix, eg for pesudo-CL
- integer intent(in) :: lmax
- integer intent(in) :: lmax_w
- integer intent(in) :: n
- real(dl) intent(in) :: W(0:lmax_w, n)
- logical intent(in) :: dopol
- real(dl) intent(out) :: M(0:lmax, 0:lmax, merge(4, n, dopol)
Legendre_Table(lmax, npoints, x, P, dP)

Legendre polynomials P_l(x_i) and derivatives dP_l(x_i)/dx for all 0 <= l <= lmax at each of the npoints x values, requiring |x| < 1. P and dP are (0:lmax, npoints) arrays (C-ordered (npoints, lmax+1) from python).
- integer intent(in) :: lmax
- integer intent(in) :: npoints
- real(dl) intent(in) :: x(npoints)
- real(dl) intent(out) :: P(0:lmax, npoints)
- real(dl) intent(out) :: dP(0:lmax, npoints)

Functions

real(dl)
airy_ai_fast(x)

Fast real Airy Ai(x), optimized for < 5e-8 absolute error. Ai-only version of airy_fast. Uses the same branch cuts, Ai polynomial coefficients, and simplified asymptotic fallbacks. The coefficients are duplicated verbatim rather than shared, to keep both routines branch-free of any present(aip) test on the hot path; a refit must update airy_fast too. camb.tests.mathutils_test checks the two agree to round-off. Maximum absolute Ai error is about 4.9e-08 on [-6.5,6.5]. On [-50,-6.5) and (6.5,25.77], the corresponding asymptotic-branch max absolute Ai errors are about 2.3e-08 and 7.9e-10. Relative errors can be large near zeros; absolute error is the intended accuracy diagnostic.
- real(dl) intent(in) :: x
real(dl)
airy_p1evl(x, coef)
- real(dl) intent(in) :: x
- real(dl) intent(in) :: coef(:)
real(dl)
airy_polevl(x, coef)
- real(dl) intent(in) :: x
- real(dl) intent(in) :: coef(:)
real(dl)
GetChiSquared(n, Y, c_inv)

get dot_product(matmul(C_inv,Y), Y) efficiently assuming c_inv symmetric
- integer intent(in) :: n
- real(dl) intent(in) :: Y(n)
- real(dl) intent(in) :: c_inv(n, n)
real(dl)
Integrate_Romberg(obj, f, a, b, tol, maxit, minsteps, abs_tol)

Rombint returns the integral from a to b of f(obj,x) using Romberg integration. The method converges provided that f is continuous in (a,b). f must be real(dl). The first argument is a class instance. tol indicates the desired relative accuracy in the integral. Modified by AL to specify max iterations and minimum number of steps (min steps useful to stop wrong results on periodic or sharp functions)
- class(*) :: obj
- procedure(obj_function) :: f f(obj, x), unwrapping obj with select type
- real(dl) intent(in) :: a
- real(dl) intent(in) :: b
- real(dl) intent(in) :: tol
- integer intent(in), optional :: maxit
- integer intent(in), optional :: minsteps
- logical intent(in), optional :: abs_tol
real(dl)
Newton_Raphson2(xxl, xxh, funcs, param, param2)
- real(dl) intent(in) :: xxl root bracket 1
- real(dl) intent(in) :: xxh root bracket 2
- external :: funcs subroutine for non-linear equation
- real(dl) intent(in) :: param parameters for function
- real(dl) intent(in) :: param2 parameters for function
real(dl)
obj_function(obj, x)

Callback taking an arbitrary class instance; the callback does select type on obj.
- class(*) :: obj
- real(dl) intent(in) :: x