FML
Typedefs | Functions | Variables
FML::MATH Namespace Reference

This namespace deals with special math functions and general math algorithms needed in the library. More...

Typedefs

using Spline = FML::INTERPOLATION::SPLINE::Spline
 
using ODESolver = FML::SOLVERS::ODESOLVER::ODESolver
 
using ODEFunction = FML::SOLVERS::ODESOLVER::ODEFunction
 
using DVector = std::vector< double >
 

Functions

std::optional< double > find_root_newton (const std::function< double(double)> &function, const std::function< double(double)> &derivative, double x0, double epsilon=1e-7)
 Find roots \( f(x) = 0 \) using Newtons method. Requires C++17. More...
 
double find_root_bisection (const std::function< double(double)> &function, std::pair< double, double > xrange, double epsilon=1e-7)
 Find roots \( f(x) = 0 \) using bisection. More...
 
double find_root_bisection (const Spline &y, double y_value, std::pair< double, double > xrange={0.0, 0.0}, double epsilon=1e-7)
 Find roots \( f(x) = 0 \) using bisection with f given as a spline. Requires GSL. More...
 
DVector Hyperspherical_j_ell_array (int lmax, const double nu, const double chi, double K)
 Hyperspherical bessel functions using recursion formula. More...
 
DVector j_ell_array (int lmax, const double x)
 Spherical bessel functions using recursion formula. More...
 
double j_ell_gsl (const int ell, const double arg)
 
double j_ell_cxx (const int ell, const double arg)
 
double j_ell (const int ell, const double arg)
 Spherical bessel function \( j_\ell(x) \) from CXX or GSL with fix for very small or large arguments. More...
 
DVector linspace (double xmin, double xmax, int num)
 Python linspace. Generate a lineary spaced array. More...
 
double Airy_Ai (double z)
 Airy function \( {\rm Ai}(x) \) (requires GSL). Found by solving and splining \( y'' - xy = 0 \). More...
 
std::pair< double, bool > GeneralizedLentzMethod (std::function< double(int)> &a, std::function< double(int)> &b, double epsilon, int maxsteps)
 Evaluate continued fraction: (b0 + a1/(b1 + a2 /( ... ))). More...
 
double Gegenbauer (int n, double x, double alpha)
 
DVector GegenbauerArray (int nmax, double x, double alpha)
 
double HyperSphericalBesselWKB (int ell, double nu, double chi, double K)
 WKB approximation for The hyper spherical bessel functions (for a curved Universe) For a flat Universe call with nu = 1.0, chi = k*(eta0-eta) and K = 0.0. More...
 
std::vector< double > legendre_ell_of_mu_vector (double mu, int ell_max)
 Legendre polynomials. More...
 

Variables

double init_airy_function_at_startup = Airy_Ai(0.0)
 

Detailed Description

This namespace deals with special math functions and general math algorithms needed in the library.

Typedef Documentation

◆ DVector

using FML::MATH::DVector = typedef std::vector<double>

Definition at line 39 of file Math.h.

◆ ODEFunction

using FML::MATH::ODEFunction = typedef FML::SOLVERS::ODESOLVER::ODEFunction

Definition at line 37 of file Math.h.

◆ ODESolver

using FML::MATH::ODESolver = typedef FML::SOLVERS::ODESOLVER::ODESolver

Definition at line 36 of file Math.h.

◆ Spline

using FML::MATH::Spline = typedef FML::INTERPOLATION::SPLINE::Spline

Definition at line 35 of file Math.h.

Function Documentation

◆ Airy_Ai()

double FML::MATH::Airy_Ai ( double  z)

Airy function \( {\rm Ai}(x) \) (requires GSL). Found by solving and splining \( y'' - xy = 0 \).

Definition at line 375 of file Math.cpp.

375 {
376 static AiryFunction * a = AiryFunction::getInstance();
377 return a->Ai(z);
378 }
example.z
float z
Definition: example.py:56

◆ find_root_bisection() [1/2]

double FML::MATH::find_root_bisection ( const Spline &  y,
double  y_value,
std::pair< double, double >  xrange,
double  epsilon 
)

Find roots \( f(x) = 0 \) using bisection with f given as a spline. Requires GSL.

Definition at line 89 of file Math.cpp.

89 {
90 const int nsearch = 20;
91
92 // The search range. If not provided use the full range from the spline
93 if (xrange.first == 0.0 && xrange.second == 0.0) {
94 xrange = y.get_xrange();
95 }
96
97 // Arrange search range so that x_low < x_high)
98 double x_low = std::min(xrange.first, xrange.second);
99 double x_high = std::max(xrange.first, xrange.second);
100
101 // We measure function values relative to y_value
102 double y_high = y(x_high) - y_value;
103 double y_low = y(x_low) - y_value;
104
105 // If the endpoints are both larger than zero do
106 // a coarse grid search to see if we can find a good value
107 // If not then throw an error
108 if (y_low * y_high >= 0.0) {
109 // Make a search grid
110 auto x_array = linspace(x_low, x_high, nsearch);
111
112 // Are the points positive or negative?
113 int sign = y_low > 0.0 ? 1 : -1;
114
115 for (int i = 0; i < nsearch; i++) {
116 double ycur = y(x_array[i]) - y_value;
117 if (sign * ycur < 0.0) {
118 y_low = ycur;
119 x_low = x_array[i];
120 break;
121 }
122 if (i == nsearch - 1) {
123 throw std::runtime_error("Error binary_search_for_value. Could not find a good interval to "
124 "start search");
125 }
126 }
127 }
128 xrange.first = x_low;
129
130 auto function = [&](double x) -> double { return (y(x) - y_value); };
131 return find_root_bisection(function, xrange, epsilon);
132 }
FML::MATH::linspace
DVector linspace(double xmin, double xmax, int num)
Python linspace. Generate a lineary spaced array.
Definition: Math.cpp:269
FML::MATH::find_root_bisection
double find_root_bisection(const Spline &y, double y_value, std::pair< double, double > xrange, double epsilon)
Find roots using bisection with f given as a spline. Requires GSL.
Definition: Math.cpp:89
FML::y
return y
Definition: Global.h:265
test.x
x
Definition: test.py:18

◆ find_root_bisection() [2/2]

double FML::MATH::find_root_bisection ( const std::function< double(double)> &  function,
std::pair< double, double >  xrange,
double  epsilon 
)

Find roots \( f(x) = 0 \) using bisection.

Definition at line 51 of file Math.cpp.

53 {
54 double x_low = std::min(xrange.first, xrange.second);
55 double x_high = std::max(xrange.first, xrange.second);
56 double y_low = function(x_low);
57 double y_high = function(x_high);
58 const int sign = y_high > y_low ? 1 : -1;
59
60 // Accuracy we want
61 const double delta_x = (x_high - x_low) * epsilon;
62
63 // Sanity check
64 if (y_low * y_high > 0.0) {
65 throw std::runtime_error("Error binary_search_for_value. f(xlow) and f(xhigh) "
66 "have the same sign!");
67 }
68
69 // Do the binary search
70 int count = 0;
71 while (x_high - x_low > delta_x) {
72 double x_mid = (x_high + x_low) / 2.0;
73 double y_mid = function(x_mid) * sign;
74 if (y_mid > 0.0) {
75 x_high = x_mid;
76 } else {
77 x_low = x_mid;
78 }
79 count++;
80 }
81
82 return x_low;
83 }

◆ find_root_newton()

std::optional< double > FML::MATH::find_root_newton ( const std::function< double(double)> &  function,
const std::function< double(double)> &  derivative,
double  x0,
double  epsilon 
)

Find roots \( f(x) = 0 \) using Newtons method. Requires C++17.

Definition at line 27 of file Math.cpp.

30 {
31 const int niter_max = 100;
32 double x = x0;
33 int iter = 0;
34 while (iter++ < niter_max) {
35 double dydx = derivative(x);
36 double y = function(x);
37 if (std::abs(y) < epsilon)
38 return x;
39 if (dydx == 0.0)
40 return {};
41 x = x - y / dydx;
42 }
43 return {};
44 }

◆ Gegenbauer()

double FML::MATH::Gegenbauer ( int  n,
double  x,
double  alpha 
)

Definition at line 455 of file Math.cpp.

455 {
456 double C0 = 1.0;
457 if (n == 0)
458 return C0;
459 double C1 = 2.0 * alpha * x;
460 if (n == 1)
461 return C1;
462 double Cn = 0.0;
463 for (int i = 2; i <= n; i++) {
464 Cn = 1.0 / i * (2 * x * (i + alpha - 1) * C1 - (i + 2 * alpha - 2) * C0);
465 C0 = C1;
466 C1 = Cn;
467 }
468 return Cn;
469 }

◆ GegenbauerArray()

DVector FML::MATH::GegenbauerArray ( int  nmax,
double  x,
double  alpha 
)

Definition at line 472 of file Math.cpp.

472 {
473 if (nmax == 0)
474 return {1.0};
475 DVector Cn(nmax + 1, 1.0);
476 Cn[1] = 2.0 * alpha * x;
477 for (int i = 2; i <= nmax; i++) {
478 Cn[i] = 1.0 / i * (2 * x * (i + alpha - 1) * Cn[i - 1] - (i + 2 * alpha - 2) * Cn[i - 2]);
479 }
480 return Cn;
481 }
DVector
FML::INTERPOLATION::SPLINE::DVector DVector
Definition: Simulation.h:39

◆ GeneralizedLentzMethod()

std::pair< double, bool > FML::MATH::GeneralizedLentzMethod ( std::function< double(int)> &  a,
std::function< double(int)> &  b,
double  epsilon,
int  maxsteps 
)

Evaluate continued fraction: (b0 + a1/(b1 + a2 /( ... ))).

For example useful for computing ratios of spherical bessel functions.

Parameters
[in]aThe function a[i]
[in]bThe function b[i]
[in]epsilonConvergence criterion
[in]maxstepsThe maximum steps before we deem it not to converge
Returns
The result and if it has converged or not

Definition at line 413 of file Math.cpp.

416 {
417
418 const double tiny = 1e-30;
419
420 double b0 = b(0);
421 if (std::fabs(b0) < tiny)
422 b0 = tiny;
423 double f = b0;
424 double C = b0;
425 double D = 0.0;
426 int j = 1;
427 for (;;) {
428 const double aa = a(j);
429 const double bb = b(j);
430 D = bb + aa * D;
431 if (std::fabs(D) < tiny)
432 D = tiny;
433 C = bb + aa / C;
434 if (std::fabs(C) < tiny)
435 C = tiny;
436 D = 1 / D;
437 const double delta = C * D;
438 f *= delta;
439 j++;
440
441 // Check for convergence
442 if (std::fabs(delta - 1) < epsilon)
443 break;
444
445 // Did not converge
446 if (j > maxsteps) {
447 std::cout << "Did not converge\n";
448 return {f, false};
449 }
450 }
451 return {f, true};
452 }

◆ Hyperspherical_j_ell_array()

DVector FML::MATH::Hyperspherical_j_ell_array ( int  lmax,
const double  nu,
const double  chi,
double  K 
)

Hyperspherical bessel functions using recursion formula.

Definition at line 140 of file Math.cpp.

140 {
141 // In this method chi is sqrt|K|(eta0-eta) and K is 0,-1,+1
142 const double chinu = chi * nu;
143 const double chinu2 = chinu * chinu;
144
145 // Geometry factors (cotc is equivalent of sinc, i.e. cot(x)*x)
146 double cotcKchi = 1.0;
147 double sincKchi = 1.0;
148 if (K == 0.0) {
149 cotcKchi = 1.0;
150 sincKchi = 1.0;
151 } else if (K < 0.0) {
152 cotcKchi = chi > 1e-8 ? chi / std::tanh(chi) : 1.0;
153 sincKchi = chi > 1e-8 ? std::sinh(chi) / chi : 1.0;
154 } else {
155 cotcKchi = chi > 1e-8 ? chi / std::tan(chi) : 1.0;
156 sincKchi = chi > 1e-8 ? std::sin(chi) / chi : 1.0;
157 }
158
159 // Start the recursion at a large enough lmax such that j_ell/j_ell-1 ~ 0
160 // and bring it down to lmax where we start to store the values
161 int lstart = std::max(lmax, int(lmax < 10 ? 5 * chinu : (lmax < 100 ? 1.6 * chinu : 1.2 * chinu)));
162 double h = 0.0;
163 double sqrtnu2chi2old = std::sqrt(chinu2 - chi * chi * K * (lstart + 1) * (lstart + 1));
164 for (int k = lstart; k >= lmax + 1; k--) {
165 double sqrtnu2chi2 = std::sqrt(chinu2 - chi * chi * K * k * k);
166 h = sqrtnu2chi2old / ((2 * k + 1) * cotcKchi - sqrtnu2chi2 * h);
167 sqrtnu2chi2old = sqrtnu2chi2;
168 }
169
170 // Recursion relation for j_(n+1) / jn
171 DVector res(lmax + 1, 0.0);
172 for (int k = lmax; k >= 1; k--) {
173 double sqrtnu2chi2 = std::sqrt(chinu2 - chi * chi * K * k * k);
174 h = sqrtnu2chi2old / ((2 * k + 1) * cotcKchi - sqrtnu2chi2 * h);
175 res[k] = h;
176 sqrtnu2chi2old = sqrtnu2chi2;
177 }
178
179 // Transform ratios into j_ell
180 res[0] = chinu == 0.0 ? (chi == 0.0 ? 1.0 : 1.0 / sincKchi) : std::sin(chinu) / (chinu * sincKchi);
181 for (int ell = 1; ell <= lmax; ell++) {
182 res[ell] *= res[ell - 1];
183 }
184
185 return res;
186 }
h
const double h
Definition: Main.cpp:66
FML::TRIANGULATION::K
CGAL::Exact_predicates_inexact_constructions_kernel K
Definition: PeriodicDelaunay.h:50

◆ HyperSphericalBesselWKB()

double FML::MATH::HyperSphericalBesselWKB ( int  ell,
double  nu,
double  chi,
double  K 
)

WKB approximation for The hyper spherical bessel functions (for a curved Universe) For a flat Universe call with nu = 1.0, chi = k*(eta0-eta) and K = 0.0.

Requires GSL.

Definition at line 484 of file Math.cpp.

484 {
485 if (chi == 0.0)
486 return ell == 0 ? 1.0 : 0.0;
487
488 //================================================================
489 // WKB approximation for the spherical bessel function
490 // Depends on the Airy function
491 // This is accurate enough for most purposes and as fast as possible
492 // Here chi = sqrt(|OmegaK|) H0 (eta0 - eta(x))
493 // nu = k / [sqrt(|OmegaK|) H0]
494 // K is the sign of the curvature i.e. 0,-1,+1
495 // For the flat case we can call this with (ell, 1.0, k*eta, 0.0)
496 //================================================================
497
498 // The geometry factor
499 double sincKchi = 1.0;
500 if (K == 0.0) {
501 sincKchi = 1.0;
502 } else if (K < 0.0) {
503 sincKchi = chi > 1e-8 ? std::sinh(chi) / chi : 1.0;
504 } else {
505 sincKchi = chi > 1e-8 ? std::sin(chi) / chi : 1.0;
506 }
507
508 // For small curvatures then alpha -> infty
509 const double sqrtellellp1 = std::sqrt(ell * (ell + 1));
510 const double alpha = nu / sqrtellellp1, a2 = alpha * alpha;
511 const double w = alpha * chi * sincKchi, w2 = w * w;
512 const double sign = w < 1.0 ? 1.0 : -1.0;
513 double S;
514 const double sqrtw2m1 = std::sqrt(sign * (1 - w2));
515 if (sign < 0.0) {
516 if (K == 0.0) {
517 S = (std::atan(1.0 / sqrtw2m1) + sqrtw2m1 - M_PI / 2.);
518 } else if (K < 0.0) {
519 const double sqrtw2pa2 = std::sqrt(w2 + a2);
520 S = (alpha * std::log((sqrtw2m1 + sqrtw2pa2) / std::sqrt(1 + a2)) +
521 std::atan(sqrtw2pa2 / sqrtw2m1) / alpha - M_PI / 2.);
522 } else {
523 const double v = std::sqrt(a2 - w2) / sqrtw2m1 / alpha;
524 S = (std::atan(v) + alpha * std::atan(1.0 / (v * alpha)) - M_PI / 2.);
525 }
526 } else {
527 if (K == 0.0) {
528 S = std::atanh(sqrtw2m1) - sqrtw2m1;
529 } else if (K < 0.0) {
530 const double u = sqrtw2m1 * alpha / std::sqrt(a2 + w2);
531 S = (atanh(u) - alpha * std::atan(u / alpha));
532 } else {
533 S = (atanh(sqrtw2m1 * alpha / std::sqrt(a2 - w2)) -
534 alpha * std::log((std::sqrt(a2 - w2) + sqrtw2m1) / std::sqrt(a2 - 1)));
535 }
536 }
537
538 const double Zpow1over6 = std::pow(1.5 * std::fabs(S) * sqrtellellp1, 1. / 6.);
539 double result = std::sqrt(M_PI / (w * sqrtw2m1)) * Zpow1over6 / sqrtellellp1;
540 result *= Airy_Ai(sign * Zpow1over6 * Zpow1over6 * Zpow1over6 * Zpow1over6);
541 return result;
542 }
FML::MATH::Airy_Ai
double Airy_Ai(double z)
Airy function (requires GSL). Found by solving and splining .
Definition: Math.cpp:375

◆ j_ell()

double FML::MATH::j_ell ( const int  ell,
const double  arg 
)

Spherical bessel function \( j_\ell(x) \) from CXX or GSL with fix for very small or large arguments.

Definition at line 248 of file Math.cpp.

248 {
249 // In this regime the function is ~1e-6 times the maximum value so put it to zero
250 // to avoid issues with the library functions failing to compute it
251 if (ell >= 10 && arg < (1.0 - 2.6 / std::sqrt(ell)) * ell)
252 return 0.0;
253
254#if ((defined(_MSVC_LANG) && _MSVC_LANG >= 201703L) || (__cplusplus >= 201703L))
255 return j_ell_cxx(ell, arg);
256#else
257#ifdef USE_GSL
258 return j_ell_gsl(ell, arg);
259#endif
260 throw std::runtime_error("No j_ell implementation availiable");
261#endif
262 }
FML::MATH::j_ell_gsl
double j_ell_gsl(const int ell, const double arg)
Definition: Math.cpp:216
FML::MATH::j_ell_cxx
double j_ell_cxx(const int ell, const double arg)
Definition: Math.cpp:233

◆ j_ell_array()

DVector FML::MATH::j_ell_array ( int  lmax,
const double  x 
)

Spherical bessel functions using recursion formula.

Definition at line 189 of file Math.cpp.

189 {
190 // Start the recursion at a large enough lmax such that j_ell/j_ell-1 ~ 0
191 // and bring it down to lmax where we start to store the values
192 int lstart = std::max(lmax, int(lmax < 10 ? 5 * x : (lmax < 100 ? 1.6 * x : 1.2 * x)));
193 double h = 0.0;
194 for (int k = lstart; k >= lmax + 1; k--) {
195 h = x / (2 * k + 1 - x * h);
196 }
197
198 // Recursion relation for j_(n+1) / jn
199 DVector res(lmax + 1, 0.0);
200 for (int k = lmax; k >= 1; k--) {
201 h = x / (2 * k + 1 - x * h);
202 res[k] = h;
203 }
204
205 // Transform ratios into j_ell
206 res[0] = x == 0.0 ? 1.0 : std::sin(x) / x;
207 for (int ell = 1; ell <= lmax; ell++) {
208 res[ell] *= res[ell - 1];
209 }
210
211 return res;
212 }

◆ j_ell_cxx()

double FML::MATH::j_ell_cxx ( const int  ell,
const double  arg 
)

Definition at line 233 of file Math.cpp.

233 {
234 // For x=0 we have 0/0
235 if (arg == 0.0)
236 return ell == 0 ? 1.0 : 0.0;
237
238 // The implementation fails for the largest arguments so simply put to zero
239 // to avoid any issues with this
240 const double arg_max_libcxx = 14000.0;
241 if (arg > arg_max_libcxx)
242 return 0.0;
243 return std::sph_bessel(ell, arg); // NOLINT
244 }

◆ j_ell_gsl()

double FML::MATH::j_ell_gsl ( const int  ell,
const double  arg 
)

Definition at line 216 of file Math.cpp.

216 {
217 // For x=0 we have 0/0
218 if (arg == 0.0)
219 return ell == 0 ? 1.0 : 0.0;
220
221 // The implementation fails for the largest arguments so simply put to zero
222 // to avoid any issues with this
223 const double arg_max_gsl = 9000.0;
224 if (ell < 500 && arg > arg_max_gsl)
225 return 0.0;
226
227 return gsl_sf_bessel_jl(ell, arg);
228 }

◆ legendre_ell_of_mu_vector()

std::vector< double > FML::MATH::legendre_ell_of_mu_vector ( double  mu,
int  ell_max 
)

Legendre polynomials.

Computes an array with \( P_0(\mu), P_1(\mu), ..., P_{ellmax}(\mu) \) using a recursion relation

Definition at line 546 of file Math.cpp.

546 {
547 std::vector<double> P_ell_array(ell_max + 1);
548 if (ell_max >= 0)
549 P_ell_array[0] = 1.0;
550 if (ell_max >= 1)
551 P_ell_array[1] = mu;
552 for (int ell = 2; ell <= ell_max; ell++)
553 P_ell_array[ell] =
554 ((2 * ell - 1.) * mu * P_ell_array[ell - 1] - (ell - 1.) * P_ell_array[ell - 2]) / ell;
555 return P_ell_array;
556 }

◆ linspace()

DVector FML::MATH::linspace ( double  xmin,
double  xmax,
int  num 
)

Python linspace. Generate a lineary spaced array.

Definition at line 269 of file Math.cpp.

269 {
270 DVector res(num);
271 double delta_x = (xmax - xmin) / double(num - 1);
272 for (int i = 0; i < num - 1; i++) {
273 res[i] = xmin + delta_x * i;
274 }
275 // Just to make sure its exactly xmax at the endpoint
276 res[num - 1] = xmax;
277 return res;
278 }
FML::MATH::DVector
std::vector< double > DVector
Definition: Math.h:39

Variable Documentation

◆ init_airy_function_at_startup

double FML::MATH::init_airy_function_at_startup = Airy_Ai(0.0)

Definition at line 383 of file Math.cpp.

Generated by doxygen 1.9.3