Utility.h

#ifndef IPPL_RANDOM_UTILITY_H
#define IPPL_RANDOM_UTILITY_H
 
#include <Kokkos_MathematicalConstants.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
#include <Kokkos_Random.hpp>
 
#include "Types/ViewTypes.h"
 
namespace ippl {
    namespace random {
        namespace detail {
            template <typename T, class Distribution>
            struct NewtonRaphson {
                Distribution dist;
                double atol           = 1e-12;
                unsigned int max_iter = 20;
 
                KOKKOS_FUNCTION
                NewtonRaphson() = default;
 
                KOKKOS_FUNCTION
                ~NewtonRaphson() = default;
 
                KOKKOS_INLINE_FUNCTION NewtonRaphson(const Distribution& dist_)
                    : dist(dist_) {}
 
                KOKKOS_INLINE_FUNCTION void solve(unsigned int d, T& x, T& u) {
                    unsigned int iter = 0;
                    while (iter < max_iter && Kokkos::fabs(dist.getObjFunc(x, d, u)) > atol) {
                        // Find x, such that "cdf(x) - u = 0" for a given sample of u~uniform(0,1)
                        x = x - (dist.getObjFunc(x, d, u) / dist.getDerObjFunc(x, d));
                        iter += 1;
                    }
                }
            };
        }  // namespace detail
    }      // namespace random
}  // namespace ippl
 
#endif