4. SecDecUtil

SecDecUtil is a standalone autotools-c++ package, that collects common helper classes and functions needed by the c++ code generated using loop_package or make_package. Everything defined by the SecDecUtil is put into the c++ namepace secdecutil.

4.1. Amplitude

A collection of utilities for evaluating amplitudes (sums of integrals multiplied by coefficients).

4.1.1. WeightedIntegral

A class template containing an integral, I, and the coefficient of the integral, C. A WeightedIntegral is interpreted as the product C*I and can be used to represent individual terms in an amplitude.

template<typename integral_t, typename coefficient_t>
struct WeightedIntegral
std::shared_ptr<integral_t> integral;

A shared pointer to the integral.

coefficient_t coefficient;

The coefficient which will be multiplied on to the integral.

std::string display_name = "WINTEGRAL";

A string used to indicate the name of the current weighted integral.

WeightedIntegral(const std::shared_ptr<integral_t> &integral, const coefficient_t &coefficient = coefficient_t(1))

The arithmetic operators (+, -, *, /) are overloaded for WeightedIntegral types.

4.1.2. WeightedIntegralHandler

A class template for integrating a sum of WeightedIntegral types.

template<typename integrand_return_t, typename real_t, typename coefficient_t, template<typename...> class container_t> class WeightedIntegralHandler
bool verbose

Controls the verbosity of the output of the amplitude.

real_t min_decrease_factor

If the next refinement iteration is expected to make the total time taken for the code to run longer than wall_clock_limit then the number of points to be requested in the next iteration will be reduced by at least min_decrease_factor.

real_t decrease_to_percentage

If remaining_time * decrease_to_percentage > time_for_next_iteration then the number of points requested in the next refinement iteration will be reduced. Here: remaining_time = wall_clock_limit - elapsed_time and time_for_next_iteration is the estimated time required for the next refinement iteration. Note: if this condition is met this means that the expected precision will not match the desired precision.

real_t wall_clock_limit

If the current elapsed time has passed wall_clock limit and a refinement iteration finishes then a new refinement iteration will not be started. Instead, the code will return the current result and exit.

size_t number_of_threads

The number of threads used to compute integrals concurrently. Note: The integrals themselves may also be computed with multiple threads irrespective of this option.

size_t reset_cuda_after

The cuda driver does not automatically remove unnecessary functions from the device memory such that the device may run out of memory after some time. This option controls after how many integrals cudaDeviceReset() is called to clear the memory. With the default 0, cudaDeviceReset() is never called. This option is ignored if compiled without cuda.

const container_t<std::vector<term_t>> &expression

The sum of terms to be integrated.

real_t epsrel

The desired relative accuracy for the numerical evaluation of the weighted sum of the sectors.

real_t epsabs

The desired absolute accuracy for the numerical evaluation of the weighted sum of the sectors.

unsigned long long int maxeval

The maximal number of integrand evaluations for each sector.

unsigned long long int mineval

The minimal number of integrand evaluations for each sector.

real_t maxincreasefac

The maximum factor by which the number of integrand evaluations will be increased in a single refinement iteration.

real_t min_epsrel

The minimum relative accuracy required for each individual sector.

real_t min_epsabs

The minimum absolute accuracy required for each individual sector.

real_t max_epsrel

The maximum relative accuracy assumed possible for each individual sector. Any sector known to this precision will not be refined further. Note: if this condition is met this means that the expected precision will not match the desired precision.

real_t max_epsabs

The maximum absolute accuracy assumed possible for each individual sector. Any sector known to this precision will not be refined further. Note: if this condition is met this means that the expected precision will not match the desired precision.

ErrorMode errormode

With enum ErrorMode : int { abs=0, all, largest, real, imag};

Defines how epsrel and epsabs are defined for complex values. With the choice largest, the relative uncertainty is defined as max( |Re(error)|, |Im(error)|)/max( |Re(result)|, |Im(result)|). Choosing all will apply epsrel and epsabs to both the real and imaginary part separately.

4.2. Series

A class template for containing (optionally truncated) Laurent series. Multivariate series can be represented as series of series.

This class overloads the arithmetic operators (+, -, *, /) and the comparator operators (==, !=). A string representation can be obtained using the << operator. The at(i) and [i] operators return the coefficient of the ith power of the expansion parameter. Otherwise elements can be accessed identically to std::vector.

template<typename T>
class Series
std::string expansion_parameter

A string representing the expansion parameter of the series (default x)

int get_order_min() const

Returns the lowest order in the series.

int get_order_max() const

Returns the highest order in the series.

bool get_truncated_above() const

Checks whether the series is truncated from above.

bool has_term(int order) const

Checks whether the series has a term at order order.

Series(int order_min, int order_max, std::vector<T> content, bool truncated_above = true, const std::string expansion_parameter = "x")

Example:

#include <iostream>
#include <secdecutil/series.hpp>

int main()
{
    secdecutil::Series<int> exact(-2,1,{1,2,3,4},false,"eps");
    secdecutil::Series<int> truncated(-2,1,{1,2,3,4},true,"eps");
    secdecutil::Series<secdecutil::Series<int>> multivariate(1,2,
                                                             {
                                                                 {-2,-1,{1,2},false,"alpha"},
                                                                 {-2,-1,{3,4},false,"alpha"},
                                                             },false,"eps"
                                                             );

    std::cout << "exact:        " << exact << std::endl;
    std::cout << "truncated:    " << truncated << std::endl;
    std::cout << "multivariate: " << multivariate << std::endl << std::endl;

    std::cout << "exact + 1:         " << exact + 1 << std::endl;
    std::cout << "exact * exact:     " << exact * exact << std::endl;
    std::cout << "exact * truncated: " << exact * truncated << std::endl;
    std::cout << "exact.at(-2):      " << exact.at(-2) << std::endl;
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -std=c++14 example.cpp -o example -lm && ./example

Output:

exact:         + (1)*eps^-2 + (2)*eps^-1 + (3) + (4)*eps
truncated:     + (1)*eps^-2 + (2)*eps^-1 + (3) + (4)*eps + O(eps^2)
multivariate:  + ( + (1)*alpha^-2 + (2)*alpha^-1)*eps + ( + (3)*alpha^-2 + (4)*alpha^-1)*eps^2

exact + 1:          + (1)*eps^-2 + (2)*eps^-1 + (4) + (4)*eps
exact * exact:      + (1)*eps^-4 + (4)*eps^-3 + (10)*eps^-2 + (20)*eps^-1 + (25) + (24)*eps + (16)*eps^2
exact * truncated:  + (1)*eps^-4 + (4)*eps^-3 + (10)*eps^-2 + (20)*eps^-1 + O(eps^0)
exact.at(-2):      1

4.3. Deep Apply

A general concept to apply a std::function to a nested data structure. If the applied std::function is not void then deep_apply() returns a nested data structure of the return values. Currently secdecutil implements this for std::vector and Series.

This concept allows, for example, the elements of a nested series to be edited without knowing the depth of the nested structure.

template<typename Tout, typename Tin, template<typename...> class Tnest>
Tnest<Tout> deep_apply(Tnest<Tin> &nest, std::function<Tout(Tin)> &func)

Example (complex conjugate a Series):

#include <iostream>
#include <complex>
#include <secdecutil/series.hpp>
#include <secdecutil/deep_apply.hpp>

int main()
{
    std::function<std::complex<double>(std::complex<double>)> conjugate =
    [] (std::complex<double> element)
    {
        return std::conj(element);
    };

    secdecutil::Series<std::complex<double>> u(-1,0,{{1,2},{3,4}},false,"eps");
    secdecutil::Series<secdecutil::Series<std::complex<double>>> m(1,1,{{1,1,{{1,2}},false,"alpha"},},false,"eps");

    std::cout << "u: " << u << std::endl;
    std::cout << "m: " << m << std::endl << std::endl;

    std::cout << "conjugated u:   " << secdecutil::deep_apply(u, conjugate) << std::endl;
    std::cout << "conjugated m: " << secdecutil::deep_apply(m, conjugate) << std::endl;
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -std=c++14 example.cpp -o example -lm && ./example

Output:

u:  + ((1,2))*eps^-1 + ((3,4))
m:  + ( + ((1,2))*alpha)*eps

conjugated u:    + ((1,-2))*eps^-1 + ((3,-4))
conjugated m:  + ( + ((1,-2))*alpha)*eps

4.4. Uncertainties

A class template which implements uncertainty propagation for uncorrelated random variables by overloads of the +, -, * and partially /. Division by UncorrelatedDeviation is not implemented as it is not always defined. It has special overloads for std::complex<T>.

Note

Division by UncorrelatedDeviation is not implemented as this operation is not always well defined. Specifically, it is ill defined in the case that the errors are Gaussian distributed as the expectation value,

\mathrm{E}\left[\frac{1}{X}\right] = \int_{-\infty}^{\infty} \frac{1}{X} p(X)\ \mathrm{d}X,

where

p(X) = \frac{1}{\sqrt{2 \pi \sigma^2 }} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right),

is undefined in the Riemann or Lebesgue sense. The rule \delta(a/b) = |a/b| \sqrt{ (\delta a/a)^2 + (\delta b/b)^2 } can not be derived from the first principles of probability theory.

The rules implemented for real valued error propagation are:

\delta(a+b) = \sqrt{(\delta a)^2 + (\delta b)^2},

\delta(a-b) = \sqrt{(\delta a)^2 + (\delta b)^2},

\delta(ab) = \sqrt{ (\delta a)^2 b^2 + (\delta b)^2 a^2 + (\delta a)^2 (\delta b)^2 }.

For complex numbers the above rules are implemented for the real and imaginary parts individually.

template<typename T>
class UncorrelatedDeviation
T value

The expectation value.

T uncertainty

The standard deviation.

Example:

#include <iostream>
#include <complex>
#include <secdecutil/uncertainties.hpp>

int main()
{
    secdecutil::UncorrelatedDeviation<double> r(1.,0.5);
    secdecutil::UncorrelatedDeviation<std::complex<double>> c({2.,3.},{0.6,0.7});

    std::cout << "r: " << r << std::endl;
    std::cout << "c: " << c << std::endl << std::endl;

    std::cout << "r.value:       " << r.value << std::endl;
    std::cout << "r.uncertainty: " << r.uncertainty << std::endl;
    std::cout << "r + c:         " << r + c << std::endl;
    std::cout << "r * c:         " << r * c << std::endl;
    std::cout << "r / 3.0:       " << r / 3. << std::endl;
    // std::cout << "1. / r:     " << 1. / r << std::endl; // ERROR
    // std::cout << "c / r:      " << c / r << std::endl;  // ERROR
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -std=c++14 example.cpp -o example -lm && ./example

Output:

r: 1 +/- 0.5
c: (2,3) +/- (0.6,0.7)

r.value:       1
r.uncertainty: 0.5
r + c:         (3,3) +/- (0.781025,0.7)
r * c:         (2,3) +/- (1.20416,1.69189)
r / 3.0:       0.333333 +/- 0.166667

4.5. Integrand Container

A class template for containing integrands. It stores the number of integration variables and the integrand as a std::function.

This class overloads the arithmetic operators (+, -, *, /) and the call operator (()).

template<typename T, typename ...Args>
class IntegrandContainer
int number_of_integration_variables

The number of integration variables that the integrand depends on.

std::function<T(Args...)> integrand

The integrand function. The call operator forwards to this function.

Example (add two IntegrandContainer and evaluate one point):

#include <iostream>
#include <secdecutil/integrand_container.hpp>

int main()
{
    using input_t = const double * const;
    using return_t = double;

    const std::function<return_t(input_t, secdecutil::ResultInfo*)> f1 = [] (input_t x, secdecutil::ResultInfo* result_info) { return 2*x[0]; };
    secdecutil::IntegrandContainer<return_t,input_t> c1(1,f1);

    const std::function<return_t(input_t, secdecutil::ResultInfo*)> f2 = [] (input_t x, secdecutil::ResultInfo* result_info) { return x[0]*x[1]; };
    secdecutil::IntegrandContainer<return_t,input_t> c2(2,f2);

    secdecutil::IntegrandContainer<return_t,input_t> c3 = c1 + c2;
    
    const double point[]{1.0,2.0};
    const double parameters[]{};
    secdecutil::ResultInfo* result_info;

    std::cout << "c1.number_of_integration_variables: " << c1.number_of_integration_variables << std::endl;
    std::cout << "c2.number_of_integration_variables: " << c2.number_of_integration_variables << std::endl << std::endl;
    std::cout << "c3.number_of_integration_variables: " << c3.number_of_integration_variables << std::endl;
    std::cout << "c3.integrand(point, parameters, result_info): " << c3.integrand_with_parameters(point, parameters, result_info) << std::endl;
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -std=c++14 example.cpp -o example -lm && ./example

Output:

c1.number_of_integration_variables: 1
c2.number_of_integration_variables: 2

c3.number_of_integration_variables: 2
c3.integrand(point, parameters, result_info): 4

4.6. Integrator

A base class template from which integrator implementations inherit. It defines the minimal API available for all integrators.

template<typename return_t, typename input_t, typename container_t = secdecutil::IntegrandContainer<return_t, input_t const*const>>
class Integrator
bool together

(Only available if return_t is a std::complex type) If true after each call of the function both the real and imaginary parts are passed to the underlying integrator. If false after each call of the function only the real or imaginary part is passed to the underlying integrator. For some adaptive integrators considering the real and imaginary part of a complex function separately can improve the sampling. Default: false.

UncorrelatedDeviation<return_t> integrate(const IntegrandContainer<return_t, input_t const*const>&)

Integrates the IntegrandContainer and returns the value and uncertainty as an UncorrelatedDeviation.

An integrator that chooses another integrator based on the dimension of the integrand.

template<typename return_t, typename input_t>
class MultiIntegrator
Integrator<return_t, input_t> &low_dimensional_integrator

Reference to the integrator to be used if the integrand has a lower dimension than critical_dim.

Integrator<return_t, input_t> &high_dimensional_integrator

Reference to the integrator to be used if the integrand has dimension critical_dim or higher.

int critical_dim

The dimension below which the low_dimensional_integrator is used.

4.6.1. CQuad

For one dimensional integrals, we wrap the cquad integrator form the GNU scientifc library (gsl).

CQuad takes the following options:

  • epsrel - The desired relative accuracy for the numerical evaluation. Default: 0.01.

  • epsabs - The desired absolute accuracy for the numerical evaluation. Default: 1e-7.

  • n - The size of the workspace. This value can only be set in the constructor. Changing this attribute of an instance is not possible. Default: 100.

  • verbose - Whether or not to print status information. Default: false.

  • zero_border - The minimal value an integration variable can take. Default: 0.0. (new in version 1.3)

4.6.2. Qmc

The quasi-monte carlo integrator as described in [PSD18]. Using a quasi-monte integrator to compute sector decomposed integrals was pioneered in [LWY+15].

template<typename return_t, ::integrators::U maxdim, template<typename, typename, ::integrators::U> class transform_t, typename container_t = secdecutil::IntegrandContainer<return_t, typename remove_complex<return_t>::type const*const>, template<typename, typename, ::integrators::U> class fitfunction_t = void_template>
class Qmc : Integrator<return_t, return_t, container_t>, public ::integrators::Qmc<return_t, return_t, maxdim, transform_t, fitfunction_t>

Derived from secdecutil::Integrator and ::integrators::Qmc - the underlying standalone implementation of the Qmc.

The most important fields and template argments of Qmc are:

  • minn - The minimal number of points in the Qmc lattice. Will be augmented to the next larger available n.

  • minm - The minimal number of random shifts.

  • maxeval - The maximal number of integrand evaluations.

  • epsrel - The desired relative accuracy for the numerical evaluation.

  • epsabs - The desired absolute accuracy for the numerical evaluation.

  • maxdim - The highest dimension the Qmc instance can be used for.

  • transform_t - The periodizing transform to apply prior to integration.

  • fitfunction_t - The fit function transform to apply for adaptive integration.

  • verbosity - Controls the amount of status messages during integration. Can be 0, 1, 2, or 3.

  • devices - A std::set of devices to run on. -1 denotes the CPU, positive integers refer to GPUs.

Refer to the documentation of the standalone Qmc for the default values and additional information.

An integral transform has to be chosen by setting the template argument transform_t. Available transforms are e.g. Korobov<r0,r1> and Sidi<r0>, please refer to the underlying Qmc implementation for a complete list. The fit function for adaptive integration can be set by the fitfunction_t, e.g. PolySingular. If not set, the default of the underlying Qmc implementation is used.

Examples how to use the Qmc on the CPU and on both, CPU and GPU are shown below.

4.6.3. Cuba

Currently we wrap the following Cuba integrators:

  • Vegas

  • Suave

  • Divonne

  • Cuhre

The Cuba integrators all implement:

  • epsrel - The desired relative accuracy for the numerical evaluation. Default: 0.01.

  • epsabs - The desired absolute accuracy for the numerical evaluation. Default: 1e-7.

  • flags - Sets the Cuba verbosity flags. The flags=2 means that the Cuba input parameters and the result after each iteration are written to the log file of the numerical integration. Default: 0.

  • seed - The seed used to generate random numbers for the numerical integration with Cuba. Default: 0.

  • mineval - The number of evaluations which should at least be done before the numerical integrator returns a result. Default: 0.

  • maxeval - The maximal number of evaluations to be performed by the numerical integrator. Default: 1000000.

  • zero_border - The minimal value an integration variable can take. Default: 0.0. (new in version 1.3)

The available integrator specific parameters and their default values are:

Vegas

Suave

Divonne

Cuhre

nstart (10000)

nnew (1000)

key1 (2000)

key (0)

nincrease (5000)

nmin (10)

key2 (1)

nbatch (500)

flatness (25.0)

key3 (1)

maxpass (4)

border (0.0)

maxchisq (1.0)

mindeviation (0.15)

For the description of these more specific parameters we refer to the Cuba manual.

4.6.4. Examples

4.6.4.1. Integrate Real Function with Cuba Vegas

Example:

#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/cuba.hpp>

int main()
{
    using input_t = const double * const;
    using return_t = double;

    secdecutil::cuba::Vegas<return_t> integrator;
    integrator.epsrel = 1e-4;
    integrator.maxeval = 1e7;

    secdecutil::IntegrandContainer<return_t,input_t> c(2, [] (input_t x, secdecutil::ResultInfo* result_info) { return x[0]*x[1]; });
    secdecutil::UncorrelatedDeviation<return_t> result = integrator.integrate(c);

    std::cout << "result: " << result << std::endl;
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lcuba -lm && ./example

Output:

result: 0.250004 +/- 2.43875e-05

4.6.4.2. Integrate Complex Function with Cuba Vegas

Example:

#include <iostream>
#include <complex>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/cuba.hpp>

int main()
{
    using input_t = const double * const;
    using return_t = std::complex<double>;

    secdecutil::cuba::Vegas<return_t> integrator;
    const std::function<return_t(input_t, secdecutil::ResultInfo*)> f = [] (input_t x, secdecutil::ResultInfo* result_info) { return return_t{x[0],x[1]}; };
    secdecutil::IntegrandContainer<return_t,input_t> c(2,f);

    integrator.together = false; // integrate real and imaginary part separately (default)
    secdecutil::UncorrelatedDeviation<return_t> result_separate = integrator.integrate(c);

    integrator.together = true; // integrate real and imaginary part simultaneously
    secdecutil::UncorrelatedDeviation<return_t> result_together = integrator.integrate(c);

    std::cout << "result_separate: " << result_separate << std::endl;
    std::cout << "result_together: " << result_together << std::endl;
    
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lcuba -lm && ./example

Output:

result_separate: (0.499937,0.499937) +/- (0.00288675,0.00288648)
result_together: (0.499937,0.499937) +/- (0.00288675,0.00288648)

4.6.4.3. Integrate Real Function with Cuba Vegas or CQuad

Example:

#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/integrator.hpp>
#include <secdecutil/integrators/cuba.hpp>
#include <secdecutil/integrators/cquad.hpp>

int main()
{
    using input_base_t = double;
    using input_t = const input_base_t * const;
    using return_t = double;

    secdecutil::cuba::Vegas<return_t> vegas;
    vegas.epsrel = 1e-5;
    vegas.maxeval = 1e7;

    secdecutil::gsl::CQuad<return_t> cquad;
    cquad.epsrel = 1e-10;
    cquad.epsabs = 1e-13;

    secdecutil::MultiIntegrator<return_t,input_base_t> integrator(cquad,vegas,2);

    secdecutil::IntegrandContainer<return_t,input_t> one_dimensional(1, [] (input_t x, secdecutil::ResultInfo* result_info) { return x[0]; });
    secdecutil::IntegrandContainer<return_t,input_t> two_dimensional(2, [] (input_t x, secdecutil::ResultInfo* result_info) { return x[0]*x[1]; });

    secdecutil::UncorrelatedDeviation<return_t> result_1d = integrator.integrate(one_dimensional); // uses cquad
    secdecutil::UncorrelatedDeviation<return_t> result_2d = integrator.integrate(two_dimensional); // uses vegas

    std::cout << "result_1d: " << result_1d << std::endl;
    std::cout << "result_2d: " << result_2d << std::endl;
}

Compile/Run:

$ c++ -I${SECDEC_CONTRIB}/include -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lcuba -lgsl -lgslcblas -lm && ./example

Output:

result_1d: 0.5 +/- 9.58209e-17
result_2d: 0.25 +/- 5.28257e-06

4.6.4.4. Set the integral transform of the Qmc

Example:

#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/qmc.hpp>

using input_base_t = double;
using input_t = input_base_t const * const;
using return_t = double;
using container_t = secdecutil::IntegrandContainer<return_t,input_t>;
using result_t = secdecutil::UncorrelatedDeviation<return_t>;

const int seed = 12345, maxdim = 4;

int main()
{
    /*
     * minimal instantiation
     */
    secdecutil::integrators::Qmc
    <
        return_t, // the return type of the integrand
        maxdim, // the highest dimension this integrator will be used for
        ::integrators::transforms::Baker::type // the integral transform
    > integrator_baker;
    integrator_baker.randomgenerator.seed(seed);

    /*
     * disable adaptation
     */
    secdecutil::integrators::Qmc
    <
        return_t, // the return type of the integrand
        maxdim, // the highest dimension this integrator will be used for
        ::integrators::transforms::Korobov<4,1>::type, // the integral transform
        container_t, // the functor type to be passed to this integrator
        ::integrators::fitfunctions::None::type // the fit funtion
    > integrator_korobov4x1;
    integrator_korobov4x1.randomgenerator.seed(seed);

    /*
     * enable adaptation
     */
    secdecutil::integrators::Qmc
    <
        return_t, // the return type of the integrand
        maxdim, // the highest dimension this integrator will be used for
        ::integrators::transforms::Sidi<3>::type, // the integral transform
        container_t, // the functor type to be passed to this integrator
        ::integrators::fitfunctions::PolySingular::type // the fit funtion
    > integrator_sidi3_adaptive;
    integrator_sidi3_adaptive.randomgenerator.seed(seed);

    // define the integrand as a functor
    container_t integrand(
                             4, // dimension
                             [] (input_t x, secdecutil::ResultInfo* result_info) { return x[0]*x[1]*x[2]*x[3]; } // integrand function
                         );

    // compute the integral with different settings
    result_t result_baker = integrator_baker.integrate(integrand);
    result_t result_korobov4x1 = integrator_korobov4x1.integrate(integrand);
    result_t result_sidi3_adaptive = integrator_sidi3_adaptive.integrate(integrand);

    // print the results
    std::cout << "baker: " << result_baker << std::endl;
    std::cout << "Korobov (weights 4, 1): " << result_korobov4x1 << std::endl;
    std::cout << "Sidi (weight 3, adaptive): " << result_sidi3_adaptive << std::endl;
}

Compile/Run:

c++ -I${SECDEC_CONTRIB}/include -pthread -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lm && ./example

Output:

baker: 0.0625 +/- 7.93855e-08
Korobov (weights 4, 1): 0.0625108 +/- 2.97931e-05
Sidi (weight 3, adaptive): 0.0625 +/- 4.33953e-09

4.6.4.5. Run the Qmc on GPUs

Example:

#include <iostream>
#include <secdecutil/integrand_container.hpp>
#include <secdecutil/uncertainties.hpp>
#include <secdecutil/integrators/qmc.hpp>

using input_base_t = double;
using input_t = input_base_t const * const;
using return_t = double;
using container_t = secdecutil::IntegrandContainer<return_t,input_t>;
using result_t = secdecutil::UncorrelatedDeviation<return_t>;

/*
 * `container_t` cannot be used on the GPU --> define a different container type
 */
struct cuda_integrand_t
{
    const static unsigned number_of_integration_variables = 4;

    // integrand function
    #ifdef __CUDACC__
        __host__ __device__
    #endif
    return_t operator()(input_t x)
    {
        return x[0]*x[1]*x[2]*x[3];
    };

    void process_errors() const{ /* error handling */}
} cuda_integrand;

const int seed = 12345, maxdim = 4;

int main()
{
    /*
     * Qmc capable of sampling on the GPU
     */
    secdecutil::integrators::Qmc
    <
        return_t, // the return type of the integrand
        maxdim, // the highest dimension this integrator will be used for
        ::integrators::transforms::Sidi<3>::type, // the integral transform
        cuda_integrand_t, // the functor type to be passed to this integrator
        ::integrators::fitfunctions::PolySingular::type // the fit funtion (optional)
    > integrator_sidi3_adaptive_gpu;
    integrator_sidi3_adaptive_gpu.randomgenerator.seed(seed);

    // compute the integral with different settings
    result_t result_sidi3_adaptive_gpu = integrator_sidi3_adaptive_gpu.integrate(cuda_integrand);

    // print the results
    std::cout << "Sidi (weight 3, adaptive): " << result_sidi3_adaptive_gpu << std::endl;
}

Compile/Run:

nvcc -x cu -I${SECDEC_CONTRIB}/include -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lgsl -lgslcblas -lm && ./example # with GPU
c++ -I${SECDEC_CONTRIB}/include -pthread -L${SECDEC_CONTRIB}/lib -std=c++14 example.cpp -o example -lgsl -lgslcblas -lm && ./example # without GPU

Output:

Sidi (weight 3, adaptive): 0.0625 +/- 4.33953e-09