From 527d54fd67899b4ff8cc578cb8a0ffb2247c3319 Mon Sep 17 00:00:00 2001
From: Alexander Voigt <Alexander.Voigt@physik.rwth-aachen.de>
Date: Wed, 4 Jul 2018 08:50:26 +0200
Subject: [PATCH] mark threshold functions and dilog as pure

to allow for improved optimization in HSSUSY
---
 src/dilog.cpp                    |  12 +-
 src/dilog.hpp                    |  11 +-
 src/threshold_loop_functions.cpp | 208 +++++++++++++++----------------
 src/threshold_loop_functions.hpp | 120 +++++++++---------
 4 files changed, 177 insertions(+), 174 deletions(-)

diff --git a/src/dilog.cpp b/src/dilog.cpp
index a4a71a804..60a2a15ad 100644
--- a/src/dilog.cpp
+++ b/src/dilog.cpp
@@ -25,7 +25,7 @@ namespace gm2calc {
 
 namespace {
    template <typename T>
-   T sqr(T x) { return x*x; }
+   T sqr(T x) noexcept { return x*x; }
 } // namespace
 
 /**
@@ -34,7 +34,7 @@ namespace {
  * @note Implementation translated by R.Brun from CERNLIB DILOG function C332
  * @return \f$\mathrm{Li}_2(z)\f$
  */
-double dilog(double x) {
+double dilog(double x) noexcept {
    const double PI = M_PI;
    const double HF  = 0.5;
    const double PI2 = PI*PI;
@@ -110,7 +110,7 @@ double dilog(double x) {
  *
  * @note not full long double precision yet!
  */
-long double dilog(long double x) {
+long double dilog(long double x) noexcept {
    return dilog(static_cast<double>(x));
 }
 
@@ -120,7 +120,7 @@ long double dilog(long double x) {
  * @note Implementation translated from SPheno to C++
  * @return \f$\mathrm{Li}_2(z)\f$
  */
-std::complex<long double> dilog(const std::complex<long double>& z) {
+std::complex<long double> dilog(const std::complex<long double>& z) noexcept {
    using flexiblesusy::fast_log;
    const long double PI = 3.1415926535897932384626433832795l;
    std::complex<long double> cy, cz;
@@ -215,7 +215,7 @@ std::complex<long double> dilog(const std::complex<long double>& z) {
  * @note Implementation translated from SPheno to C++
  * @return \f$\mathrm{Li}_2(z)\f$
  */
-std::complex<double> dilog(const std::complex<double>& z) {
+std::complex<double> dilog(const std::complex<double>& z) noexcept {
    const auto rz = static_cast<long double>(std::real(z));
    const auto iz = static_cast<long double>(std::imag(z));
    const auto re = dilog(std::complex<long double>(rz, iz));
@@ -227,7 +227,7 @@ std::complex<double> dilog(const std::complex<double>& z) {
  * @param x real angle
  * @return \f$\mathrm{Cl}_2(\theta)\f$
  */
-double clausen_2(double x)
+double clausen_2(double x) noexcept
 {
    using std::exp;
    using gm2calc::dilog;
diff --git a/src/dilog.hpp b/src/dilog.hpp
index 4d1b4d2ad..3782f34ec 100644
--- a/src/dilog.hpp
+++ b/src/dilog.hpp
@@ -20,23 +20,24 @@
 #define DILOG_H
 
 #include <complex>
+#include "cextensions.hpp"
 
 namespace gm2calc {
 
 /// real dilogarithm
-double dilog(double);
+double dilog(double) noexcept ATTR(const);
 
 /// real dilogarithm
-long double dilog(long double);
+long double dilog(long double) noexcept ATTR(const);
 
 /// complex dilogarithm
-std::complex<double> dilog(const std::complex<double>&);
+std::complex<double> dilog(const std::complex<double>&) noexcept ATTR(const);
 
 /// complex dilogarithm
-std::complex<long double> dilog(const std::complex<long double>&);
+std::complex<long double> dilog(const std::complex<long double>&) noexcept ATTR(const);
 
 /// Clausen function Cl_2(x)
-double clausen_2(double);
+double clausen_2(double) noexcept ATTR(const);
 
 } // namespace gm2calc
 
diff --git a/src/threshold_loop_functions.cpp b/src/threshold_loop_functions.cpp
index 5920468d6..923c378fa 100644
--- a/src/threshold_loop_functions.cpp
+++ b/src/threshold_loop_functions.cpp
@@ -31,30 +31,30 @@ namespace threshold_loop_functions {
 namespace {
    const double Pi = 3.1415926535897932384626433832795;
 
-   template <typename T> T sqr(T x) { return x*x; }
-   template <typename T> T cube(T x) { return x*x*x; }
-   template <typename T> T quad(T x) { return x*x*x*x; }
-   template <typename T> T pow5(T x) { return x*x*x*x*x; }
-   template <typename T> T pow6(T x) { return x*x*x*x*x*x; }
-   template <typename T> T pow7(T x) { return x*x*x*x*x*x*x; }
-   template <typename T> T pow8(T x) { return x*x*x*x*x*x*x*x; }
-   template <typename T> T pow9(T x) { return x*x*x*x*x*x*x*x*x; }
-   template <typename T> T power10(T x) { return x*x*x*x*x*x*x*x*x*x; }
+   template <typename T> T sqr(T x) noexcept { return x*x; }
+   template <typename T> T cube(T x) noexcept { return x*x*x; }
+   template <typename T> T quad(T x) noexcept { return x*x*x*x; }
+   template <typename T> T pow5(T x) noexcept { return x*x*x*x*x; }
+   template <typename T> T pow6(T x) noexcept { return x*x*x*x*x*x; }
+   template <typename T> T pow7(T x) noexcept { return x*x*x*x*x*x*x; }
+   template <typename T> T pow8(T x) noexcept { return x*x*x*x*x*x*x*x; }
+   template <typename T> T pow9(T x) noexcept { return x*x*x*x*x*x*x*x*x; }
+   template <typename T> T power10(T x) noexcept { return x*x*x*x*x*x*x*x*x*x; }
 
    template <typename T>
-   bool is_zero(T a, T prec = std::numeric_limits<T>::epsilon())
+   bool is_zero(T a, T prec = std::numeric_limits<T>::epsilon()) noexcept
    {
       return std::fabs(a) < prec;
    }
 
    template <typename T>
-   bool is_equal(T a, T b, T prec = std::numeric_limits<T>::epsilon())
+   bool is_equal(T a, T b, T prec = std::numeric_limits<T>::epsilon()) noexcept
    {
       return is_zero(a - b, prec);
    }
 
    template <typename T>
-   bool is_equal_rel(T a, T b, T prec = std::numeric_limits<T>::epsilon())
+   bool is_equal_rel(T a, T b, T prec = std::numeric_limits<T>::epsilon()) noexcept
    {
       if (is_equal(a, b, std::numeric_limits<T>::epsilon()))
          return true;
@@ -67,7 +67,7 @@ namespace {
 
 } // anonymous namespace
 
-double F1(double x)
+double F1(double x) noexcept
 {
    const double x2 = sqr(x);
 
@@ -80,7 +80,7 @@ double F1(double x)
    return x*std::log(x2)/(x2-1);
 }
 
-double F2(double x)
+double F2(double x) noexcept
 {
    const double x2 = sqr(x);
 
@@ -93,7 +93,7 @@ double F2(double x)
    return 6*x2*(2-2*x2+(1+x2)*std::log(x2))/cube(x2-1);
 }
 
-double F3(double x)
+double F3(double x) noexcept
 {
    const double x2 = sqr(x);
 
@@ -106,7 +106,7 @@ double F3(double x)
    return 2*x*(5*(1-x2)+(1+4*x2)*std::log(x2))/(3*sqr(x2-1));
 }
 
-double F4(double x)
+double F4(double x) noexcept
 {
    const double x2 = sqr(x);
 
@@ -119,7 +119,7 @@ double F4(double x)
    return 2*x*(x2-1-std::log(x2))/sqr(x2-1);
 }
 
-double F5(double x)
+double F5(double x) noexcept
 {
    const double x2 = sqr(x);
    const double x4 = quad(x);
@@ -136,7 +136,7 @@ double F5(double x)
    return 3*x*(1-x4+2*x2*std::log(x2))/cube(1-x2);
 }
 
-double F6(double x)
+double F6(double x) noexcept
 {
    const double x2 = sqr(x);
 
@@ -152,7 +152,7 @@ double F6(double x)
    return (x2-3)/(4*(1-x2)) + x2*(x2-2)/(2*sqr(1.-x2))*std::log(x2);
 }
 
-double F7(double x)
+double F7(double x) noexcept
 {
    const double x2 = sqr(x);
    const double x4 = quad(x);
@@ -175,14 +175,14 @@ double F7(double x)
 }
 
 /// F8(x1,x2) in the limit x1 -> 1 and x2 -> 1
-static double F8_1_1(double x1, double x2)
+static double F8_1_1(double x1, double x2) noexcept
 {
    return 1. + 1.333333333333333*(-1 + x1) +
       (1.333333333333333 - 0.6666666666666661*(-1 + x1))*(-1 + x2);
 }
 
 /// F8(x1,x2) in the limit x1 -> 1
-static double F8_1_x2(double x1, double x2)
+static double F8_1_x2(double x1, double x2) noexcept
 {
    const double lx22 = std::log(sqr(x2));
 
@@ -210,7 +210,7 @@ static double F8_1_x2(double x1, double x2)
 }
 
 /// F8(x1,x2) in the limit x1 -> 0
-static double F8_0_x2(double x1, double x2)
+static double F8_0_x2(double x1, double x2) noexcept
 {
    const double lx22 = std::log(sqr(x2));
 
@@ -219,7 +219,7 @@ static double F8_0_x2(double x1, double x2)
 }
 
 // F8(x1,x2) in the limit x1 -> x2
-static double F8_x1_x2(double x1, double x2)
+static double F8_x1_x2(double x1, double x2) noexcept
 {
    const double lx22 = std::log(sqr(x2));
 
@@ -241,7 +241,7 @@ static double F8_x1_x2(double x1, double x2)
     (x2*pow5(-1. + sqr(x2)));
 }
 
-double F8(double x1, double x2)
+double F8(double x1, double x2) noexcept
 {
    if (is_equal(x1, 0.) && is_equal(x2, 0.))
       return -2.;
@@ -288,7 +288,7 @@ double F8(double x1, double x2)
 }
 
 /// F9(x1,x2) in the limit x1 -> 1 and x2 -> 1
-static double F9_1_1(double x1, double x2)
+static double F9_1_1(double x1, double x2) noexcept
 {
    return 8.223809523809523 - 12.863492063492064*x2
       + 10.580952380952382*sqr(x2) - 4.609523809523809*cube(x2)
@@ -311,7 +311,7 @@ static double F9_1_1(double x1, double x2)
 }
 
 /// F9(x1,x2) in the limit x1 -> 1
-static double F9_1_x2(double x1, double x2)
+static double F9_1_x2(double x1, double x2) noexcept
 {
    const double lx22 = std::log(sqr(x2));
 
@@ -335,13 +335,13 @@ static double F9_1_x2(double x1, double x2)
 }
 
 /// F9(x1,x2) in the limit x1 -> 0
-static double F9_0_x2(double, double x2)
+static double F9_0_x2(double, double x2) noexcept
 {
    return (2.*std::log(sqr(x2)))/(-1. + sqr(x2));
 }
 
 /// F9(x1,x2) in the limit x1 -> x2
-static double F9_x1_x2(double x1, double x2)
+static double F9_x1_x2(double x1, double x2) noexcept
 {
    const double lx22 = std::log(sqr(x2));
 
@@ -363,7 +363,7 @@ static double F9_x1_x2(double x1, double x2)
       (sqr(x2)*quad(-1. + sqr(x2)));
 }
 
-double F9(double x1, double x2)
+double F9(double x1, double x2) noexcept
 {
    if (is_equal(x1, 1., 0.01) && is_equal(x2, 1., 0.01))
       return F9_1_1(x1, x2);
@@ -400,17 +400,17 @@ double F9(double x1, double x2)
    return 2./(x12-x22)*(x12/(x12-1.)*std::log(x12)-x22/(x22-1.)*std::log(x22));
 }
 
-double f(double r)
+double f(double r) noexcept
 {
    return F5(r);
 }
 
-double g(double r)
+double g(double r) noexcept
 {
    return F7(r);
 }
 
-double f1(double r)
+double f1(double r) noexcept
 {
    if (is_equal(r, 0., 0.01))
       return 18./7.*sqr(r);
@@ -431,7 +431,7 @@ double f1(double r)
       + (6*(r2-5)*quad(r)*std::log(r2))/(7*cube(r2-1));
 }
 
-double f2(double r)
+double f2(double r) noexcept
 {
    if (is_equal(r, 0., 0.01))
       return 22./9.*sqr(r);
@@ -452,7 +452,7 @@ double f2(double r)
       + (2*(5*r2-17)*quad(r)*std::log(r2))/(9*cube(r2-1));
 }
 
-double f3(double r)
+double f3(double r) noexcept
 {
    if (is_equal(r, 0., 0.001))
       return 4./3.;
@@ -472,7 +472,7 @@ double f3(double r)
       + (2*(r4-7*r2-6)*r2*std::log(r2))/(3*cube(r2-1));
 }
 
-double f4(double r)
+double f4(double r) noexcept
 {
    if (is_equal(r, 0., 0.001))
       return 12./7.;
@@ -491,7 +491,7 @@ double f4(double r)
 }
 
 /// f5(r1,r2) in the limit r1 -> 1 and r2 -> 1
-static double f5_1_1(double r1, double r2)
+static double f5_1_1(double r1, double r2) noexcept
 {
    return 0.772943722943723
       - 0.5524891774891774*r2
@@ -520,7 +520,7 @@ static double f5_1_1(double r1, double r2)
 }
 
 /// f5(r1,r2) in the limit r1 -> 1
-static double f5_1_r2(double r1, double r2)
+static double f5_1_r2(double r1, double r2) noexcept
 {
    const double lr22 = std::log(sqr(r2));
 
@@ -560,7 +560,7 @@ static double f5_1_r2(double r1, double r2)
 }
 
 /// f5(r1,r2) in the limit r1 -> 0
-static double f5_0_r2(double r1, double r2)
+static double f5_0_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -572,13 +572,13 @@ static double f5_0_r2(double r1, double r2)
 }
 
 /// f5(r1,r2) in the limit r1 -> 0 and r2 -> 1
-static double f5_0_1(double, double r2)
+static double f5_0_1(double, double r2) noexcept
 {
    return 0.75*(1 + (-1 + r2)/3. + sqr(-1 + r2)/6.);
 }
 
 /// f5(r1,r2) in the limit r1 -> r2
-static double f5_r1_r2(double r1, double r2)
+static double f5_r1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -605,7 +605,7 @@ static double f5_r1_r2(double r1, double r2)
     (6.*r2*pow6(-1 + r22));
 }
 
-double f5(double r1, double r2)
+double f5(double r1, double r2) noexcept
 {
    if (is_equal(r1, 0., 0.0001) && is_equal(r2, 0., 0.0001))
       return 0.75;
@@ -651,7 +651,7 @@ double f5(double r1, double r2)
 }
 
 /// f6(r1,r2) in the limit r1 -> 1 and r2 -> 1
-static double f6_1_1(double r1, double r2)
+static double f6_1_1(double r1, double r2) noexcept
 {
    return 1 + (4*(-1 + r1))/7. - (2*sqr(-1 + r1))/35.
       - cube(-1 + r1)/70. + (9*quad(-1 + r1))/490.
@@ -670,7 +670,7 @@ static double f6_1_1(double r1, double r2)
 }
 
 /// f6(r1,r2) in the limit r1 -> 1
-static double f6_1_r2(double r1, double r2)
+static double f6_1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
 
@@ -699,7 +699,7 @@ static double f6_1_r2(double r1, double r2)
 }
 
 /// f6(r1,r2) in the limit r1 -> 0
-static double f6_0_r2(double r1, double r2)
+static double f6_0_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -710,7 +710,7 @@ static double f6_0_r2(double r1, double r2)
 }
 
 // f6(r1,r2) in the limit r1 -> r2
-static double f6_r1_r2(double r1, double r2)
+static double f6_r1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -736,13 +736,13 @@ static double f6_r1_r2(double r1, double r2)
 }
 
 /// f6(r1,r2) in the limit r1 -> 0 and r2 -> 1
-static double f6_0_1(double, double r2)
+static double f6_0_1(double, double r2) noexcept
 {
    return 6./7.*(0.5 + (2*(-1 + r2))/3. - cube(-1 + r2)/15.
                  + quad(-1 + r2)/20.);
 }
 
-double f6(double r1, double r2)
+double f6(double r1, double r2) noexcept
 {
    if (is_equal(r1, 0., 0.0001) && is_equal(r2, 0., 0.0001))
       return 0.;
@@ -788,7 +788,7 @@ double f6(double r1, double r2)
 }
 
 /// f7(r1,r2) in the limit r1 -> 1 and r2 -> 1
-static double f7_1_1(double r1, double r2)
+static double f7_1_1(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
 
@@ -806,7 +806,7 @@ static double f7_1_1(double r1, double r2)
 }
 
 /// f7(r1,r2) in the limit r1 -> 1
-static double f7_1_r2(double r1, double r2)
+static double f7_1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -832,7 +832,7 @@ static double f7_1_r2(double r1, double r2)
 }
 
 /// f7(r1,r2) in the limit r1 -> 0
-static double f7_0_r2(double r1, double r2)
+static double f7_0_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -843,14 +843,14 @@ static double f7_0_r2(double r1, double r2)
 }
 
 /// f7(r1,r2) in the limit r1 -> 0 and r2 -> 1
-static double f7_0_1(double, double r2)
+static double f7_0_1(double, double r2) noexcept
 {
    return 6.*(0.5 + (1 - r2)/3. + sqr(-1 + r2)/6. - cube(-1 + r2)/15.
               + quad(-1 + r2)/60.);
 }
 
 /// f7(r1,r2) in the limit r1 -> r2
-static double f7_r1_r2(double r1, double r2)
+static double f7_r1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -876,7 +876,7 @@ static double f7_r1_r2(double r1, double r2)
     (6.*r2*pow6(-1 + r22));
 }
 
-double f7(double r1, double r2)
+double f7(double r1, double r2) noexcept
 {
    if (is_equal(r1, 0., 0.0001) && is_equal(r2, 0., 0.0001))
       return 6.;
@@ -922,7 +922,7 @@ double f7(double r1, double r2)
 }
 
 /// f8(r1,r2) in the limit r1 -> 1 and r2 -> 1
-static double f8_1_1(double r1, double r2)
+static double f8_1_1(double r1, double r2) noexcept
 {
    return 1 - sqr(-1 + r1)/10. + (3*cube(-1 + r1))/40.
       - (3*quad(-1 + r1))/70.
@@ -941,7 +941,7 @@ static double f8_1_1(double r1, double r2)
 }
 
 /// f8(r1,r2) in the limit r1 -> 1
-static double f8_1_r2(double r1, double r2)
+static double f8_1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -969,7 +969,7 @@ static double f8_1_r2(double r1, double r2)
 }
 
 /// f8(r1,r2) in the limit r1 -> 0
-static double f8_0_r2(double r1, double r2)
+static double f8_0_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -980,14 +980,14 @@ static double f8_0_r2(double r1, double r2)
 }
 
 /// f8(r1,r2) in the limit r1 -> 0 and r2 -> 1
-static double f8_0_1(double, double r2)
+static double f8_0_1(double, double r2) noexcept
 {
    return 1.5*(0.5 + (-1 + r2)/6. - sqr(-1 + r2)/6. + cube(-1 + r2)/10.
                - quad(-1 + r2)/20.);
 }
 
 /// f8(r1,r2) in the limit r1 -> r2
-static double f8_r1_r2(double r1, double r2)
+static double f8_r1_r2(double r1, double r2) noexcept
 {
    const double r22 = sqr(r2);
    const double lr22 = std::log(r22);
@@ -1010,7 +1010,7 @@ static double f8_r1_r2(double r1, double r2)
     (6.*pow6(-1 + r22));
 }
 
-double f8(double r1, double r2)
+double f8(double r1, double r2) noexcept
 {
    if (is_equal(r1, 0., 0.0001) && is_equal(r2, 0., 0.0001))
       return 0.;
@@ -1056,7 +1056,7 @@ double f8(double r1, double r2)
 }
 
 /// First derivative of F1
-double D1F1(double x)
+double D1F1(double x) noexcept
 {
    using std::log;
 
@@ -1067,7 +1067,7 @@ double D1F1(double x)
 }
 
 /// First derivative of F2
-double D1F2(double x)
+double D1F2(double x) noexcept
 {
    using std::log;
 
@@ -1079,7 +1079,7 @@ double D1F2(double x)
 }
 
 /// First derivative of F3
-double D1F3(double x)
+double D1F3(double x) noexcept
 {
    using std::log;
 
@@ -1091,7 +1091,7 @@ double D1F3(double x)
 }
 
 /// First derivative of F4
-double D1F4(double x)
+double D1F4(double x) noexcept
 {
    using std::log;
 
@@ -1104,7 +1104,7 @@ double D1F4(double x)
 }
 
 /// First derivative of F5
-double D1F5(double x)
+double D1F5(double x) noexcept
 {
    using std::log;
 
@@ -1116,7 +1116,7 @@ double D1F5(double x)
 }
 
 /// First derivative of F6
-double D1F6(double x)
+double D1F6(double x) noexcept
 {
    using std::log;
 
@@ -1127,7 +1127,7 @@ double D1F6(double x)
 }
 
 /// First derivative of F7
-double D1F7(double x)
+double D1F7(double x) noexcept
 {
    using std::log;
 
@@ -1139,7 +1139,7 @@ double D1F7(double x)
 }
 
 /// First derivative of f
-double D1f(double x)
+double D1f(double x) noexcept
 {
    using std::log;
 
@@ -1151,7 +1151,7 @@ double D1f(double x)
 }
 
 /// First derivative of g
-double D1g(double x)
+double D1g(double x) noexcept
 {
    using std::log;
 
@@ -1163,7 +1163,7 @@ double D1g(double x)
 }
 
 /// First derivative of f1
-double D1f1(double x)
+double D1f1(double x) noexcept
 {
    using std::log;
 
@@ -1175,7 +1175,7 @@ double D1f1(double x)
 }
 
 /// First derivative of f2
-double D1f2(double x)
+double D1f2(double x) noexcept
 {
    using std::log;
 
@@ -1187,7 +1187,7 @@ double D1f2(double x)
 }
 
 /// First derivative of f3
-double D1f3(double x)
+double D1f3(double x) noexcept
 {
    using std::log;
 
@@ -1199,7 +1199,7 @@ double D1f3(double x)
 }
 
 /// First derivative of f4
-double D1f4(double x)
+double D1f4(double x) noexcept
 {
    using std::log;
 
@@ -1211,7 +1211,7 @@ double D1f4(double x)
 }
 
 /// First derivative of f5 w.r.t. 1st argument
-double D10f5(double x, double y)
+double D10f5(double x, double y) noexcept
 {
    using std::log;
 
@@ -1253,7 +1253,7 @@ double D10f5(double x, double y)
 }
 
 /// First derivative of f5 w.r.t. 2nd argument
-double D01f5(double x, double y)
+double D01f5(double x, double y) noexcept
 {
    using std::log;
 
@@ -1295,7 +1295,7 @@ double D01f5(double x, double y)
 }
 
 /// First derivative of f6 w.r.t. 1st argument
-double D10f6(double x, double y)
+double D10f6(double x, double y) noexcept
 {
    using std::log;
 
@@ -1335,7 +1335,7 @@ double D10f6(double x, double y)
 }
 
 /// First derivative of f6 w.r.t. 2nd argument
-double D01f6(double x, double y)
+double D01f6(double x, double y) noexcept
 {
    using std::log;
 
@@ -1374,7 +1374,7 @@ double D01f6(double x, double y)
 }
 
 /// First derivative of f7 w.r.t. 1st argument
-double D10f7(double x, double y)
+double D10f7(double x, double y) noexcept
 {
    using std::log;
 
@@ -1413,7 +1413,7 @@ double D10f7(double x, double y)
 }
 
 /// First derivative of f7 w.r.t. 2nd argument
-double D01f7(double x, double y)
+double D01f7(double x, double y) noexcept
 {
    using std::log;
 
@@ -1452,7 +1452,7 @@ double D01f7(double x, double y)
 }
 
 /// First derivative of f8 w.r.t. 1st argument
-double D10f8(double x, double y)
+double D10f8(double x, double y) noexcept
 {
    using std::log;
 
@@ -1491,7 +1491,7 @@ double D10f8(double x, double y)
 }
 
 /// First derivative of f8 w.r.t. 2nd argument
-double D01f8(double x, double y)
+double D01f8(double x, double y) noexcept
 {
    using std::log;
 
@@ -1530,7 +1530,7 @@ double D01f8(double x, double y)
 }
 
 /// Second derivative of F1
-double D2F1(double x)
+double D2F1(double x) noexcept
 {
    using std::log;
 
@@ -1543,7 +1543,7 @@ double D2F1(double x)
 }
 
 /// Second derivative of F2
-double D2F2(double x)
+double D2F2(double x) noexcept
 {
    using std::log;
 
@@ -1555,7 +1555,7 @@ double D2F2(double x)
 }
 
 /// Second derivative of F3
-double D2F3(double x)
+double D2F3(double x) noexcept
 {
    using std::log;
 
@@ -1568,7 +1568,7 @@ double D2F3(double x)
 }
 
 /// Second derivative of F4
-double D2F4(double x)
+double D2F4(double x) noexcept
 {
    using std::log;
 
@@ -1580,7 +1580,7 @@ double D2F4(double x)
 }
 
 /// Second derivative of F5
-double D2F5(double x)
+double D2F5(double x) noexcept
 {
    using std::log;
 
@@ -1592,7 +1592,7 @@ double D2F5(double x)
 }
 
 /// Second derivative of F6
-double D2F6(double x)
+double D2F6(double x) noexcept
 {
    using std::log;
 
@@ -1604,7 +1604,7 @@ double D2F6(double x)
 }
 
 /// Second derivative of F7
-double D2F7(double x)
+double D2F7(double x) noexcept
 {
    using std::log;
 
@@ -1616,14 +1616,14 @@ double D2F7(double x)
 }
 
 /// Iabc(a,a,a)
-static double Iaaa(double a, double b, double c)
+static double Iaaa(double a, double b, double c) noexcept
 {
    return (151.*quad(a) + 13.*sqr(b)*sqr(c) - 128.*cube(a)*(b + c) - 40.*a*b*c*(b + c)
            + sqr(a)*(37.*sqr(b) + 128.*b*c + 37.*sqr(c))) / (60.*pow6(a));
 }
 
 /// Iabc(a,a,c)
-static double Iaac(double a, double b, double c)
+static double Iaac(double a, double b, double c) noexcept
 {
    return ((sqr(a) - sqr(c))
            * (17.*pow6(a) - 16.*pow5(a)*b - 40.*cube(a)*b*sqr(c)
@@ -1635,7 +1635,7 @@ static double Iaac(double a, double b, double c)
 }
 
 /// Iabc(a,a,0)
-double Iaa0(double a, double b)
+double Iaa0(double a, double b) noexcept
 {
    return (17.*sqr(a) - 16.*a*b + 5.*sqr(b)) / (6.*quad(a));
 }
@@ -1646,7 +1646,7 @@ double I0bc(double b, double c)
    return log(sqr(b/c))/(sqr(b) - sqr(c));
 }
 
-double Iabc(double a, double b, double c) {
+double Iabc(double a, double b, double c) noexcept {
    if ((is_zero(a) && is_zero(b) && is_zero(c)) ||
        (is_zero(a) && is_zero(b)) ||
        (is_zero(a) && is_zero(c)) ||
@@ -1690,20 +1690,20 @@ double Iabc(double a, double b, double c) {
 }
 
 /// Delta function from hep-ph/0907.47682v1
-double delta_xyz(double x, double y, double z)
+double delta_xyz(double x, double y, double z) noexcept
 {
    return sqr(x)+sqr(y)+sqr(z)-2*(x*y+x*z+y*z);
 }
 
 namespace {
    /// lambda^2(u,v)
-   double lambda_2(double u, double v)
+   double lambda_2(double u, double v) noexcept
    {
       return sqr(1 - u - v) - 4*u*v;
    }
 
    /// u < 1 && v < 1, lambda^2(u,v) > 0
-   double phi_pos(double u, double v)
+   double phi_pos(double u, double v) noexcept
    {
       using std::log;
       using gm2calc::dilog;
@@ -1717,7 +1717,7 @@ namespace {
    }
 
    /// lambda^2(u,v) < 0
-   double phi_neg(double u, double v)
+   double phi_neg(double u, double v) noexcept
    {
       using std::acos;
       using std::sqrt;
@@ -1735,7 +1735,7 @@ namespace {
     * The following identities hold:
     * Phi(u,v) = Phi(v,u) = Phi(1/u,v/u)/u = Phi(1/v,u/v)/v
     */
-   double phi_uv(double u, double v)
+   double phi_uv(double u, double v) noexcept
    {
       const auto lambda = lambda_2(u,v);
 
@@ -1764,7 +1764,7 @@ namespace {
  *
  * @return \f$\Phi(x,y,z)\f$
  */
-double phi_xyz(double x, double y, double z)
+double phi_xyz(double x, double y, double z) noexcept
 {
    const auto u = x/z, v = y/z;
    return phi_uv(u,v);
@@ -1779,7 +1779,7 @@ double phi_xyz(double x, double y, double z)
  *
  * @return \f$B_0(p=0,m_1,m_2,Q)\f$
  */
-double B0(double m1, double m2, double scale)
+double B0(double m1, double m2, double scale) noexcept
 {
    return passarino_veltman::ReB0(0, m1*m1, m2*m2, scale*scale);
 }
@@ -1792,7 +1792,7 @@ double B0(double m1, double m2, double scale)
  *
  * @return \f$B_0'(p=0,m_1,m_2)\f$
  */
-double DB0(double m1, double m2)
+double DB0(double m1, double m2) noexcept
 {
    const double m12 = sqr(m1);
    const double m14 = sqr(m12);
@@ -1820,12 +1820,12 @@ double DB0(double m1, double m2)
  *
  * @return \f$C_0(p=0,m_1,m_2,m_3)\f$
  */
-double C0(double m1, double m2, double m3)
+double C0(double m1, double m2, double m3) noexcept
 {
    return softsusy::c0(m1, m2, m3);
 }
 
-double D0(double m1, double m2, double m3, double m4)
+double D0(double m1, double m2, double m3, double m4) noexcept
 {
    return softsusy::d0(m1,m2,m3,m4);
 }
@@ -1840,7 +1840,7 @@ double D0(double m1, double m2, double m3, double m4)
  *
  * @return \f$\tilde{D}_2(m_1,m_2,m_3,m_4)\f$
  */
-double D2t(double m1, double m2, double m3, double m4)
+double D2t(double m1, double m2, double m3, double m4) noexcept
 {
    return C0(m2, m3, m4) + m1*m1 * D0(m1, m2, m3, m4);
 }
@@ -1856,7 +1856,7 @@ double D2t(double m1, double m2, double m3, double m4)
  *
  * @return \f$\tilde{D}_4(m_1,m_2,m_3,m_4,Q)\f$
  */
-double D4t(double m1, double m2, double m3, double m4, double scale)
+double D4t(double m1, double m2, double m3, double m4, double scale) noexcept
 {
    return B0(m3, m4, scale) + (m1*m1 + m2*m2) * C0(m2, m3, m4)
       + quad(m1) * D0(m1, m2, m3, m4);
@@ -1871,7 +1871,7 @@ double D4t(double m1, double m2, double m3, double m4, double scale)
  *
  * @return \f$W(m_1,m_2,Q)\f$
  */
-double W(double m1, double m2, double scale)
+double W(double m1, double m2, double scale) noexcept
 {
    const double m12 = sqr(m1);
    const double m14 = sqr(m12);
diff --git a/src/threshold_loop_functions.hpp b/src/threshold_loop_functions.hpp
index 6b7df8716..b2ea1acdf 100644
--- a/src/threshold_loop_functions.hpp
+++ b/src/threshold_loop_functions.hpp
@@ -41,96 +41,98 @@
 #define TCD1f0 threshold_loop_functions::D1f
 #define TCD1g0 threshold_loop_functions::D1g
 
+#include "cextensions.hpp"
+
 namespace flexiblesusy {
 namespace threshold_loop_functions {
 
 // loop functions from arXiv:1407.4081
 
-double F1(double);
-double F2(double);
-double F3(double);
-double F4(double);
-double F5(double);
-double F6(double);
-double F7(double);
-double F8(double, double);
-double F9(double, double);
-
-double f(double);
-double g(double);
-
-double f1(double);
-double f2(double);
-double f3(double);
-double f4(double);
-double f5(double, double);
-double f6(double, double);
-double f7(double, double);
-double f8(double, double);
+double F1(double) noexcept ATTR(const);
+double F2(double) noexcept ATTR(const);
+double F3(double) noexcept ATTR(const);
+double F4(double) noexcept ATTR(const);
+double F5(double) noexcept ATTR(const);
+double F6(double) noexcept ATTR(const);
+double F7(double) noexcept ATTR(const);
+double F8(double, double) noexcept ATTR(const);
+double F9(double, double) noexcept ATTR(const);
+
+double f(double) noexcept ATTR(const);
+double g(double) noexcept ATTR(const);
+
+double f1(double) noexcept ATTR(const);
+double f2(double) noexcept ATTR(const);
+double f3(double) noexcept ATTR(const);
+double f4(double) noexcept ATTR(const);
+double f5(double, double) noexcept ATTR(const);
+double f6(double, double) noexcept ATTR(const);
+double f7(double, double) noexcept ATTR(const);
+double f8(double, double) noexcept ATTR(const);
 
 // first derivatives
 
-double D1F1(double);
-double D1F2(double);
-double D1F3(double);
-double D1F4(double);
-double D1F5(double);
-double D1F6(double);
-double D1F7(double);
-double D1f(double);
-double D1g(double);
-double D1f1(double);
-double D1f2(double);
-double D1f3(double);
-double D1f4(double);
-double D10f5(double, double);
-double D01f5(double, double);
-double D10f6(double, double);
-double D01f6(double, double);
-double D10f7(double, double);
-double D01f7(double, double);
-double D10f8(double, double);
-double D01f8(double, double);
+double D1F1(double) noexcept ATTR(const);
+double D1F2(double) noexcept ATTR(const);
+double D1F3(double) noexcept ATTR(const);
+double D1F4(double) noexcept ATTR(const);
+double D1F5(double) noexcept ATTR(const);
+double D1F6(double) noexcept ATTR(const);
+double D1F7(double) noexcept ATTR(const);
+double D1f(double) noexcept ATTR(const);
+double D1g(double) noexcept ATTR(const);
+double D1f1(double) noexcept ATTR(const);
+double D1f2(double) noexcept ATTR(const);
+double D1f3(double) noexcept ATTR(const);
+double D1f4(double) noexcept ATTR(const);
+double D10f5(double, double) noexcept ATTR(const);
+double D01f5(double, double) noexcept ATTR(const);
+double D10f6(double, double) noexcept ATTR(const);
+double D01f6(double, double) noexcept ATTR(const);
+double D10f7(double, double) noexcept ATTR(const);
+double D01f7(double, double) noexcept ATTR(const);
+double D10f8(double, double) noexcept ATTR(const);
+double D01f8(double, double) noexcept ATTR(const);
 
 // second derivatives
 
-double D2F1(double);
-double D2F2(double);
-double D2F3(double);
-double D2F4(double);
-double D2F5(double);
-double D2F6(double);
-double D2F7(double);
+double D2F1(double) noexcept ATTR(const);
+double D2F2(double) noexcept ATTR(const);
+double D2F3(double) noexcept ATTR(const);
+double D2F4(double) noexcept ATTR(const);
+double D2F5(double) noexcept ATTR(const);
+double D2F6(double) noexcept ATTR(const);
+double D2F7(double) noexcept ATTR(const);
 
 /// \f$I_{abc}(a,b,c)\f$ (arguments are interpreted as unsquared)
-double Iabc(double, double, double);
+double Iabc(double, double, double) noexcept ATTR(const);
 
 /// \f$Delta_{xyz}(x,y,z)\f$ (arguments are interpreted as squared masses)
-double delta_xyz(double, double, double);
+double delta_xyz(double, double, double) noexcept ATTR(const);
 
 /// \f$phi_{xyz}(x,y,z)\f$ (arguments are interpreted as squared masses)
-double phi_xyz(double, double, double);
+double phi_xyz(double, double, double) noexcept ATTR(const);
 
 /// \f$B_0(p=0,m_1,m_2,Q)\f$ (arguments are interpreted as unsquared)
-double B0(double, double, double);
+double B0(double, double, double) noexcept ATTR(const);
 
 /// \f$B_0'(p=0,m_1,m_2)\f$ (arguments are interpreted as unsquared)
-double DB0(double, double);
+double DB0(double, double) noexcept ATTR(const);
 
 /// \f$C_0(p=0,m_1,m_2,m_3)\f$ (arguments are interpreted as unsquared)
-double C0(double, double, double);
+double C0(double, double, double) noexcept ATTR(const);
 
 /// \f$D_0(p=0,m_1,m_2,m_3,m_4)\f$ (arguments are interpreted as unsquared)
-double D0(double, double, double, double);
+double D0(double, double, double, double) noexcept ATTR(const);
 
 /// \f$\tilde{D}_2(m_1,m_2,m_3,m_4)\f$ (arguments are interpreted as unsquared)
-double D2t(double, double, double, double);
+double D2t(double, double, double, double) noexcept ATTR(const);
 
 /// \f$\tilde{D}_4(m_1,m_2,m_3,m_4,Q)\f$ (arguments are interpreted as unsquared)
-double D4t(double, double, double, double, double);
+double D4t(double, double, double, double, double) noexcept ATTR(const);
 
 /// \f$Q(m_1,m_2,Q)\f$ (arguments are interpreted as unsquared)
-double W(double, double, double);
+double W(double, double, double) noexcept ATTR(const);
 
 } // namespace threshold_loop_functions
 } // namespace flexiblesusy