CPP/C++ Math and Numerical Computing

  [S]ign  [E]xponent               [M]antissa (aka significand)
  1 bit   8 bits                   23 bits 
  +----+-------------------------+-------------------------------+
  | S  | E7 E6 E5 E4 E3 E2 E1 E0 | MMM.MMMM MMMM.MMMM MMMM.MMMM  |
  +----+-------------------------+-------------------------------+
    31 | 30                   23 | 22                           0  --- Bit position 

Mantissa: M22 M21 M20 M19 M18 M17 M16 .... M2 M1 M0
         --------                          ---------
           byte 4                            byte 0

              Sign         Exponent - Bias 
Value :=  (-1)     *  Radix        *        Mantissa

              S      E - B 
Value :=  (-1)  *  R    *  M

               S      (E - 127)
 Value :=  (-1)   *  2         * (1.M)
          -----    ------        ------
          Sign     Exponent     Mantissa 
         Factor     Factor       Factor 


              b31    b30,b29,..,b23 - 127           
Value :=  (-1)   *  2            *     (1.b22,b21,b20,...b0)

matissa = 1.0 + b22 / 2^1 + b21 / 2^2 + b20 / 2^3 + ... b0 / 2^22 

// IEEE754 - Float Point - 32 bits - single precision 
>> std::numeric_limits<float>::epsilon()
(float) 1.19209e-07f

// IEEE754 - Float Point - 64 bits - double precision 
>> std::numeric_limits<double>::epsilon()
(double) 2.2204460e-16

// IEE754 - Float point 128 bits - quadruple precision 
>> std::numeric_limits<long double>::epsilon()
(long double) 1.0842022e-19L

 // Machine precison. 
 >> constexpr auto eps = std::numeric_limits<double>::epsilon()
 (const double) 2.2204460e-16

 >> 10.0 == 10.0
 (bool) true

 // Two numbers are almost equal if their difference is less than or equal than the 
 // EPS - epsilon or the machine precison. 
 >> 10.0 == 10.0 + eps
 (bool) true
 >> 

 // =======>>> Double - 64 bits float point <<================

 >> 10.0 == 10.000000000000001
 (bool) false

 >> .000000000000001 < std::numeric_limits<double>::epsilon()
 (bool) false

 // Cannot distinguish 10.0 and 10.00000000000000001 since the difference 
 // between them is less than the EPSILON for double precision float points.
 >> 10.0 == 10.00000000000000001
 (bool) true

 >> .00000000000000001 < std::numeric_limits<double>::epsilon()
 (bool) true
 >> 
 >> .00000000000000001 < eps
(bool) true

 // =======>>> float - 32 bits floating point <<================ 

 >> 10.0f == 10.000001f
 (bool) false

 >> 0.000001f < std::numeric_limits<float>::epsilon()
 (bool) false

 // No longer can distinguish 10.0f and 10.0000001f since 
 // the difference between them is less than EPSILON for floating points. 
 >> 10.0f == 10.0000001f
 (bool) true

 >> .0000001f < std::numeric_limits<float>::epsilon()
 (bool) true

>> 0.1 + 0.1 + 0.1 == 0.3
(bool) false

>> 0.2 + 0.1 == 0.3
(bool) false

>> 1.0 / 3.0
(double) 0.33333333
>> 
>> 1.0 / 3.0 == 0.333333333
(bool) false
>

// Machine precison. - EPISILON 
>> constexpr auto eps = std::numeric_limits<double>::epsilon()
(const double) 2.2204460e-16

// NAN - Not A Number 
>> auto dnan = std::numeric_limits<double>::quiet_NaN()
(double) nan

// Positive Infinity 
>> constexpr auto pinf = std::numeric_limits<double>::infinity()
(const double) inf

// Negative Infinity 
>> constexpr auto ninf = -std::numeric_limits<double>::infinity()
(const double) -inf

>> q = sqrt(-10)
(double) -nan

// Division of zero by zero 
>> x = 0.0 / 0.0
(double) nan

>> q + 10.0
(double) nan

>> q - ninf
(double) nan

>> q + ninf
(double) nan

>> 0.0 == 0.0
(bool) true
>> 

>> 10.0 == 10.0
(bool) true

>> 9.81451654 == 9.81451654
(bool) true

>> dnan == dnan
(bool) false

>> xx = dnan;

>> double xx = 10.0;
>> if(xx) { std::puts(" [TRUE]"); } else { std::puts(" [FALSE] "); }
 [TRUE]

>> xx = -100.0;
>> if(xx) { std::puts(" [TRUE]"); } else { std::puts(" [FALSE] "); }
 [TRUE]

// Set to Zero 
>> xx = 0.0;
>> if(xx) { std::puts(" [TRUE]"); } else { std::puts(" [FALSE] "); }
 [FALSE] 

// Set to NAN 
>> xx = std::numeric_limits<double>::quiet_NaN();
>> if(xx) { std::puts(" [TRUE]"); } else { std::puts(" [FALSE] "); }
 [TRUE]

>> double x;

>> x = std::numeric_limits<double>::infinity()
(double) inf

>> x = -std::numeric_limits<double>::infinity()
(double) -inf

>> x = std::numeric_limits<double>::quiet_NaN()
(double) nan
>> 

>> x = std::numeric_limits<double>::signaling_NaN()
(double)

>>  x = std::numeric_limits<double>::max()
(double) 1.7976931e+308

>>  x = std::numeric_limits<double>::min()
(double) 2.2250739e-308

>> x = std::numeric_limits<double>::lowest()
(double) -1.7976931e+308

>> x = std::numeric_limits<double>::epsilon()
(double) 2.2204460e-16

>> float fx = std::numeric_limits<float>::epsilon()
(float) 1.19209e-07f

>> std::numeric_limits<long double>::epsilon()
(long double) 1.0842022e-19L
>> 

abs(x - y)  < tolerance

 abs(x - y)  < tolerance
 -----------
   abs(x)

Or:

   abs(x - y)  < tolerance * abs(x)

   abs(x - y)  < tolerance x max(abs(x), abs(y))

OR
   abs(x - y)  < tolerance x (abs(x) + abs(y)) / 2.0

def almostEqual(a, b, relTol, absTol):
    return abs(x - y) < max(relTol * max(abs(x), abs(y)), absTol)

/** Checks whether two floats X and Y are almost equal under given tolerance
  * Requires headers: <cmath> and <limits>
  * @tparam T - Any float point type (float, double, long double)
  */
template<typename T>
auto aeq_abs(T x, T y, T tol) -> bool
{
    constexpr auto eps = std::numeric_limits<T>::epsilon();
    // Tolerance can never be smaller than the Epsilon for this
    tol = std::max(tol, eps);
    // Missing: handle the case if X and/or Y are NaN
    return std::fabs(x - y) < tol;
}

template<typename T>
auto aeq_rel(T x, T y, T tol) -> bool
{
    constexpr auto eps = std::numeric_limits<T>::epsilon();
    // Tolerance can never be smaller than the Epsilon for this
    tol = std::max(tol, eps);
    // Missing: handle the case if X and/or Y are NaN
    return std::fabs(x - y) < tol * std::max(std::fabs(x), std::fabs(y)) ;
}

/** Checks wether two floating points are almost equal or close enough. 
 *  @tparam T        Any floating point type (float, double, long double)
 *  @param input     Input floating point value 
 *  @param expected  Expectd floating point value 
 *  @param relTol    Relative tolerance 
 *  @param absTol    Absolute tolerance 
 *  @return bool     Returns true if input and expected are close enough.
 *  Note: Requires headers: <cmath> and <limits>
 */
template<typename T>
bool almost_equal(T input, T expected
                  , T relTol, T absTol = 8 * std::numeric_limits<T>::eps())
{
    using std::fmax;
    using std::max;
    return fmax(expected - input) < max(relTol * max(input, expected), absTol);
}

/** Check whether two floating points are almost equal or close enough.
 *
 * @tparam T    Any floating point type (float, double, ...)
 * @param  x    Input floating point value
 * @param  y    Expected floating point value
 * @param  drel Precision digits of the relative tolerance.
 * @param  dabs Precision digits of the absolute tolerance.
 * @return bool Returns true if x and y are close enough.
 *
 * Requires headers: <cmath> and <limits>
 */
template<typename T>
auto almost_equal_digits(T x, T y, int drel = 8, int dabs = 0) -> bool
{
    using std::fabs;
    using std::max;
    static const T precisions [] = {
        1E-1, 1E-2, 1E-3, 1E-4, 1E-5, 1E-6, 1E-7, 1E-8
        , 1E-9, 1E-10, 1E-11, 1E-12, 1E-13, 1E-15, 1E-16 };

    constexpr size_t n = std::size(precisions);

    T tol_rel = drel > 0 && drel < n
                    ? precisions[drel - 1]
                    : std::numeric_limits<T>::epsilon();

    T tol_abs = dabs > 0 && dabs < n
                    ? precisions[dabs - 1]
                    : 8 & std::numeric_limits<T>::epsilon();

    return fabs(x - y) < max(tol_rel * max(fabs(x), fabs(y)), tol_abs);
}

$ g++ decompose_binary32flt.cpp -o out.bin -g -std=c++1z -Wall -Wextra 

#include <iostream>
#include <string>
#include <iomanip>
#include <limits>
#include <cmath>
#include <cstdint>
#include <map>

#define SHOW_FLOAT(expr) show_float(#expr, (expr))

// Decompose float point
void show_float(const char* expr, float x)
{
    constexpr int w1 = 34;
    constexpr int w2 = 4;

    // Float point classification constants
    static auto const fpmap = std::map<int, const char*>{
         {FP_INFINITE,  "FP_INFINITE"}
        ,{FP_NAN,       "FP_NAN"}
        ,{FP_NORMAL,    "PF_NORMAL"}
        ,{FP_SUBNORMAL, "FP_SUBNORMAL"}
        ,{FP_ZERO,      "FP_ZERO"}
    };

    // Integer represetation extracted
    // through type punning aka memory reinterpretation
    std::uint32_t* pn       = reinterpret_cast<std::uint32_t*>(&x);
    std::uint32_t  value    = *pn;
    std::uint32_t  mantissa = value & ((1 << 23) - 1);
    std::uint32_t  raw_exp  = (value >> 23) & 0xFF;
    std::uint32_t  sign     = (value >> 31) & 0x01;

    auto print_row = [&w1, &w2](const char* label, auto const& value)
    {
        std::cout << std::setw(w1) << std::right << label
                  << std::setw(w2) << std::left  << value
                  << "\n";
    };

    print_row("Input Expression: ", expr);
    print_row("Classification: ", fpmap.at(std::fpclassify(x)));
    print_row("Decimal Representation: ", x);
    print_row("Hexadecimal Representation: ", to_hexfloat(x));
    print_row("Integer Representation (hex): ", num2hexstr(value));
    print_row("Sign bit: ", sign);
    print_row("Exponent: ", static_cast<std::int32_t>(raw_exp) - 127);
    print_row("Raw Exponent: ", raw_exp);
    print_row("Raw Exponent (hex): ", num2hexstr(raw_exp, 0));
    print_row("Mantissa (hex): ", num2hexstr(mantissa, 0));    
    std::cout << "\n\n";
}

std::puts("\n==== Decomposition of Single Precision Float Points =====\n");

// Signed Zero
SHOW_FLOAT(+0.0f);
SHOW_FLOAT(-0.0f);

// NaN - Not a number
SHOW_FLOAT(+std::numeric_limits<float>::quiet_NaN());
SHOW_FLOAT(-std::numeric_limits<float>::quiet_NaN());

// Positive and negative Infinity
SHOW_FLOAT(+std::numeric_limits<float>::infinity());
SHOW_FLOAT(-std::numeric_limits<float>::infinity());

// Subnormal
SHOW_FLOAT(+std::numeric_limits<float>::min());
SHOW_FLOAT(+std::numeric_limits<float>::min() / 100.0f);
SHOW_FLOAT(-std::numeric_limits<float>::min() / 1000000.0f);

// Epsilon
SHOW_FLOAT(+std::numeric_limits<float>::epsilon());

// Normal Numbers
SHOW_FLOAT(1.0f);
SHOW_FLOAT(0.5f);
SHOW_FLOAT(1E-5f);
SHOW_FLOAT(1.051646E4f);
SHOW_FLOAT(-99000134.3401f);
SHOW_FLOAT(std::numeric_limits<float>::max());
SHOW_FLOAT(std::numeric_limits<float>::lowest());

$ ./out.bin 

==== Decomposition of Single Precision Float Points =====

     ...     ...     ...     ...     ...     ...     ...     ... 

                Input Expression: +std::numeric_limits<float>::infinity()
                  Classification: FP_INFINITE
          Decimal Representation: inf 
      Hexadecimal Representation: inf 
    Integer Representation (hex): 0x7f800000
                        Sign bit: 0   
                        Exponent: 128 
                    Raw Exponent: 255 
              Raw Exponent (hex): 0xff
                  Mantissa (hex): 0x0 


                Input Expression: -std::numeric_limits<float>::infinity()
                  Classification: FP_INFINITE
          Decimal Representation: -inf
      Hexadecimal Representation: -inf
    Integer Representation (hex): 0xff800000
                        Sign bit: 1   
                        Exponent: 128 
                    Raw Exponent: 255 
              Raw Exponent (hex): 0xff
                  Mantissa (hex): 0x0 


                Input Expression: +std::numeric_limits<float>::min()
                  Classification: PF_NORMAL
          Decimal Representation: 1.17549e-38
      Hexadecimal Representation: 0x1p-126
    Integer Representation (hex): 0x00800000
                        Sign bit: 0   
                        Exponent: -126
                    Raw Exponent: 1   
              Raw Exponent (hex): 0x1 
                  Mantissa (hex): 0x0 

  ... ....   ... ....   ... ....   ... ....   ... .... 

#include <iostream>
#include <random>

void generate_randoms1(size_t n){
     // Create default initialized engine object 
     std::default_random_engine engine{};
     // Generates uniformily distributed integers from 1 to 10, 
     // including both bounds of the interval.  
     std::uniform_int_distribution<int> dist(1, 10); 
     for(auto i = 0; i < n; i++)
       std::cout << " x[" << i << "] = " << dist(engine) << std::endl;
}

>> generate_randoms1(4)
 x[0] = 1
 x[1] = 2
 x[2] = 8
 x[3] = 5

>> generate_randoms1(4)
 x[0] = 1
 x[1] = 2
 x[2] = 8
 x[3] = 5

>> generate_randoms1(6)
 x[0] = 1
 x[1] = 2
 x[2] = 8
 x[3] = 5
 x[4] = 6
 x[5] = 3 

>> generate_randoms1(10)
 x[0] = 1
 x[1] = 2
 x[2] = 8
 x[3] = 5
 x[4] = 6
 x[5] = 3
 x[6] = 1
 x[7] = 7
 x[8] = 7
 x[9] = 10

#include <iostream>
#include <random>

void generate_randoms2(size_t n, unsigned int seed){
     // Create default initialized engine object 
     std::default_random_engine engine{seed};
     // Generates uniformily distributed integers from 1 to 10, 
     // including both bounds of the interval.  
     std::uniform_int_distribution<int> dist(1, 10); 
     for(auto i = 0; i < n; i++)
       std::cout << " x[" << i << "] = " << dist(engine) << std::endl;
}

// Generates a new seed whenever it is called 
std::random_device seed_gen{};

// A seed is a big integer used for initializing the generator 
>> seed_gen()
(unsigned int) 4133644392

>> seed_gen()
(unsigned int) 2900292039

// Save current seed 
>> auto seed = seed_gen()
(unsigned int) 2616205616

>> generate_randoms2(5, seed)
 x[0] = 4
 x[1] = 1
 x[2] = 5
 x[3] = 10
 x[4] = 6
>> 
>> generate_randoms2(6, seed)
 x[0] = 4
 x[1] = 1
 x[2] = 5
 x[3] = 10
 x[4] = 6
 x[5] = 10

>> generate_randoms2(10, seed)
 x[0] = 4
 x[1] = 1
 x[2] = 5
 x[3] = 10
 x[4] = 6
 x[5] = 10
 x[6] = 2
 x[7] = 8
 x[8] = 9
 x[9] = 6
>> 

>> auto seed2 = seed_gen()
(unsigned int) 3091314163


>> generate_randoms2(4, seed2)
 x[0] = 8
 x[1] = 5
 x[2] = 2
 x[3] = 7

>> generate_randoms2(8, seed2)
 x[0] = 8
 x[1] = 5
 x[2] = 2
 x[3] = 7
 x[4] = 6
 x[5] = 9
 x[6] = 4
 x[7] = 1

>> generate_randoms2(5, seed_gen())
 x[0] = 1
 x[1] = 9
 x[2] = 8
 x[3] = 10
 x[4] = 2

>> generate_randoms2(5, seed_gen())
 x[0] = 5
 x[1] = 6
 x[2] = 7
 x[3] = 8
 x[4] = 4

>> generate_randoms2(5, seed_gen())
 x[0] = 3
 x[1] = 4
 x[2] = 1
 x[3] = 8
 x[4] = 8

#include <iostream>
#include <random> 
#include <vector> 
#include <string> 
#include <functional>

int main(int argc, char** argv)
{       
    if(argc < 3)
    {
        std::cout << "Error: option provided. Try again!" << std::endl;
        return EXIT_FAILURE;        
    }   

    size_t quantity = std::stoi(argv[1]);
    std::string command = argv[2];

    // Subprogram 1 - Run without any provided seed 
    if(command == "-empty")
    {
        std::cout << " [INFO] Running without setting the generator's seed" << std::endl;       
        // Random number engine 
        std::default_random_engine rngen{};
        // mean = 10.0, stdev = 5.0 
        std::normal_distribution<double> ndist(10.0, 5.0);
        for(auto i = 0; i < quantity; i++)
            std::cout << "Next random is = " << ndist(rngen) << std::endl;
        return EXIT_SUCCESS;
    }
    // Subprogram 2 - User provides the seed 
    if(command == "-seed")
    {
        if(argc < 4){
                std::cout << "Error: missing seed. " << std::endl;
                return EXIT_FAILURE;
        }
        auto seed = std::atol(argv[3]);
        std::cout << " [INFO] Using seed = " << seed << std::endl;

        // Random number engine 
        std::default_random_engine rngen{seed};
        // mean = 10.0, stdev = 5.0 
        std::normal_distribution<double> ndist(10.0, 5.0);
        for(auto i = 0; i < quantity; i++)
           std::cout << "Next random is = " << ndist(rngen) << std::endl;               
    }
    // Subprogram 2  - Automatically generates a seed 
    if(command == "-auto")
    {
        // Object that generates seed (initial state of random generator)
        std::random_device seed_gen;

        // Get a seed (big random integer) used for initializing the random generator. 
        auto seed = seed_gen(); 
        std::cout << "seed = " <<  seed << std::endl;

        // Random number engine 
        std::default_random_engine rngen{seed};
        // mean = 10.0, stdev = 5.0 
        std::normal_distribution<double> ndist(10.0, 5.0);
        for(auto i = 0; i < quantity; i++)
           std::cout << "Next random is = " << ndist(rngen) << std::endl;               
    }       

    return 0;
}

# Using Clang 
$ clang++ randoms1.cpp -o randoms1.bin -std=c++1z -g -O0 -Wall 

# Using GCC
$ gcc++ randoms1.cpp -o randoms1.bin -std=c++1z -g -O0 -Wall 

# Fails! 

$ ./randoms1.bin 5 -empty
 [INFO] Running without setting the generator's seed
Next random is = 9.39017
Next random is = 4.56591
Next random is = 13.4214
Next random is = 4.62405
Next random is = 10.1663


$ ./randoms1.bin 5 -empty
 [INFO] Running without setting the generator's seed
Next random is = 9.39017
Next random is = 4.56591
Next random is = 13.4214
Next random is = 4.62405
Next random is = 10.1663
(base) 

$ ./randoms1.bin 8 -empty
 [INFO] Running without setting the generator's seed
Next random is = 9.39017
Next random is = 4.56591
Next random is = 13.4214
Next random is = 4.62405
Next random is = 10.1663
Next random is = 13.7242
Next random is = 10.168
Next random is = 7.36681

$ ./randoms1.bin 4 -auto
seed = 2138866676
Next random is = 6.58512
Next random is = 11.7181
Next random is = 9.55006
Next random is = 7.35961

$ ./randoms1.bin 8 -auto
seed = 4200712391
Next random is = 0.575461
Next random is = 11.82
Next random is = 6.07653
Next random is = 13.8118
Next random is = 11.0999
Next random is = 7.04851
Next random is = 7.56055
Next random is = 3.67443

$ ./randoms1.bin 3 -auto
seed = 1536173647
Next random is = 13.4482
Next random is = 12.2431
Next random is = 12.2183

$ ./randoms1.bin 4 -seed 2138866676
 [INFO] Using seed = 2138866676
Next random is = 6.58512
Next random is = 11.7181
Next random is = 9.55006
Next random is = 7.35961

$ ./randoms1.bin 8 -seed 2138866676
 [INFO] Using seed = 2138866676
Next random is = 6.58512
Next random is = 11.7181
Next random is = 9.55006
Next random is = 7.35961
Next random is = 3.621
Next random is = 8.41599
Next random is = 11.3628
Next random is = 8.27506

$ ./randoms1.bin 8 -seed 4200712391
 [INFO] Using seed = 4200712391
Next random is = 0.575461
Next random is = 11.82
Next random is = 6.07653
Next random is = 13.8118
Next random is = 11.0999
Next random is = 7.04851
Next random is = 7.56055
Next random is = 3.67443

#include <iostream>
#include <random>
#include <iomanip>
#include <string>

template<typename TInt = int>
class RandInt {
private:
    std::random_device                 m_nextSeed;  
    unsigned int                       m_seed;
    std::default_random_engine         m_engine;
    std::uniform_int_distribution<int> m_dist;
public:
     using TSeed = long;

     // Intialize with a known seed 
     RandInt(TInt min, TInt max, unsigned int seed)
      : m_seed(seed), m_engine(seed), m_dist(min, max)
     {      
     }

     // Initialize with a random seed 
     RandInt(TInt min, TInt max)
       : m_nextSeed{},
         m_seed(m_nextSeed()),
         m_engine(m_seed),
         m_dist(min, max)
     {      
     }
     TInt  Min()  const { return m_dist.min(); }
     TInt  Max()  const { return m_dist.max(); }
     TSeed Seed() const { return m_seed; }
     void  Seed(TSeed seed)
     {
          m_seed = seed;
          m_engine.seed(seed);
     }
     // Set seed to a new random (non-deterministic) value and return it 
     TSeed NextSeed()
     {
          m_seed = m_nextSeed();
          m_engine.seed(m_seed);
          return m_seed;
     }
     // Get NExt random 
     TInt operator()()
     {
          return m_dist(m_engine);
     }
};

std::cout << "\n ===== Random numbers with a random seed ===="
                  << std::endl;
{
     RandInt rnd(1, 10);
     std::cout << " => rnd.Seed() = " << rnd.Seed() << std::endl;
     std::cout << " => rnd.Min() = "  << rnd.Min() << std::endl;
     std::cout << " => rnd.Max() = "  << rnd.Max() << std::endl;

     for(int i = 0; i < 15; i++){
         auto field = std::string(" x[") + std::to_string(i) + "] = "; 
         std::cout << std::setw(10)
                   << field
                   << std::setw(5) << rnd() << std::endl;
     }      
}   

std::cout << "\n ===== Random numbers with a non-random seed ===="
          << std::endl;
{
     // Initialize Random generator object with known and fixed
     // seed to reproduce computation results. 
     long seed = 1195785783;
     RandInt rnd(1, 10, seed);
     std::cout << " => rnd.Seed() = " << rnd.Seed() << std::endl;
     std::cout << " => rnd.Min() = "  << rnd.Min() << std::endl;
     std::cout << " => rnd.Max() = "  << rnd.Max() << std::endl;

     /* Expected sequence: 7, 10, 9, 2, 2, 1, 3, 1, 6, 5, 1, 2, 3, 2 ... */ 
     for(int i = 0; i < 15; i++){
             auto field = std::string(" x[") + std::to_string(i) + "] = "; 
             std::cout << std::setw(10)
                       << field
                       << std::setw(5) << rnd() << std::endl;
     }      
}   

return 0;

$ g++ RandInt.cpp -o RandInt.bin -std=c++1z -g -O0 -Wall 
$ clang++ RandInt.cpp -o RandInt.bin -std=c++1z -g -O0 -Wall

$ clang++ RandInt.cpp -o RandInt.bin -std=c++1z -g -O0 -Wall && ./RandInt.bin

 ===== Random numbers with a random seed ====
 => rnd.Seed() = 2611535369
 => rnd.Min() = 1
 => rnd.Max() = 10
   x[0] =     9
   x[1] =     1
   x[2] =    10
   x[3] =     6
   x[4] =    10
   x[5] =     5
   x[6] =     2
   x[7] =     9
   x[8] =     9
   x[9] =    10
  x[10] =     5
  x[11] =     9
  x[12] =     8
  x[13] =     5
  x[14] =     9

 ===== Random numbers with a non-random seed ====
 => rnd.Seed() = 1195785783
 => rnd.Min() = 1
 => rnd.Max() = 10
   x[0] =     7
   x[1] =    10
   x[2] =     9
   x[3] =     2
   x[4] =     2
   x[5] =     1
   x[6] =     3
   x[7] =     1
   x[8] =     6
   x[9] =     5
  x[10] =     1
  x[11] =     2
  x[12] =     1
  x[13] =     3
  x[14] =     2

#include <iostream>
#include <random>
#include <iomanip>
#include <string>

// FUNCTION OVERLOAD 1 
//-----------------------
// Note: The return type of this function must not be set to std::function<TInt (void)>
// or that will have a performance penalty as this container uses type erasure and heap allocation.
// This function returns lambda or an anonymous callable object, in other words, an object that overloads
// the operator()(), generated by the compiler.
template<typename TInt> 
auto make_uniform_random_distribution(TInt a, TInt b, long seed) 
{       
    std::mt19937 engine(seed);
    std::uniform_int_distribution<TInt> dist(a, b);        
    return [=]() mutable -> TInt { return dist(engine); };
}

// FUNCTION OVERLOADING 2
//----------------------------
template<typename TInt> 
auto make_uniform_random_distribution(TInt a, TInt b) 
{
   std::random_device rd;
   return make_uniform_random_distribution(a, b, rd()) ;
}

#if 1 
 std::cout << "\n ===== Random numbers with a random seed ===="
           << std::endl;
 {
    auto rnd = make_uniform_random_distribution<int>(1, 10);
    for(int i = 0; i < 15; i++){
        auto field = std::string(" x[") + std::to_string(i) + "] = "; 
        std::cout << std::setw(10)
                  << field
                  << std::setw(5) << rnd() << std::endl;
        }       
 }  
#endif 

std::cout << "\n ===== Random numbers with a non-random seed ===="
          << std::endl;
{
   // Initialize Random generator object with known and fixed
   // seed to reproduce computation results. 
   unsigned int seed = 1195785783;
   auto rnd = make_uniform_random_distribution<int>(1, 10, seed);

   /* Expected sequence: 7, 10, 9, 2, 2, 1, 3, 1, 6, 5, 1, 2, 3, 2 ... */ 
   for(int i = 0; i < 15; i++)
   {
       auto field = std::string(" x[") + std::to_string(i) + "] = "; 
       std::cout << std::setw(10)
                 << field
                 << std::setw(5) << rnd() << std::endl;
    }       
}   

# Clang LLVM 
$ clang++ rand_int_fun.cpp -o rand_int_fun.bin -std=c++1z -g -O0 -Wall 

# GCC - GNU C/C++ Compiler 
$ g++ rand_int_fun.cpp -o rand_int_fun.bin -std=c++1z -g -O0 -Wall 

$ ./rand_int_fun.bin 

 ===== Random numbers with a random seed ====
   x[0] =     8
   x[1] =     4
   x[2] =     4
   x[3] =     7
   x[4] =     9
   x[5] =     2
   x[6] =     2
   x[7] =     6
   x[8] =     3
   x[9] =     8
  x[10] =     3
  x[11] =     1
  x[12] =     3
  x[13] =     1
  x[14] =    10

 ===== Random numbers with a non-random seed ====
   x[0] =     6
   x[1] =    10
   x[2] =     6
   x[3] =     2
   x[4] =     4
   x[5] =     1
   x[6] =     9
   x[7] =    10
   x[8] =     3
   x[9] =     6
  x[10] =    10
  x[11] =     6
  x[12] =     2
  x[13] =     1
  x[14] =     1

===== Random numbers with a non-random seed ====
   x[0] =     3
   x[1] =     6
   x[2] =     1
   x[3] =     5
   x[4] =     4
   x[5] =     7
   x[6] =     2
   x[7] =     6
   x[8] =     7
   x[9] =     4
  x[10] =     8
  x[11] =     7
  x[12] =     2
  x[13] =     1
  x[14] =    10

import sympy as sy 
from sympy.abc import x, y, z
from sympy import diff, expand, sin, exp, cos, tan, log 

sy.init_printing()

>>> f = sin(x) * exp(x) / (x**2  + 3 * x + 10)

>>> f
   x         
  ℯ ⋅sin(x)  
─────────────
 2           
x  + 3⋅x + 10

>>> df = diff(f, x)

>>> df
            x             x               x         
(-2⋅x - 3)⋅ℯ ⋅sin(x)     ℯ ⋅sin(x)       ℯ ⋅cos(x)  
──────────────────── + ───────────── + ─────────────
                 2      2               2           
  ⎛ 2           ⎞      x  + 3⋅x + 10   x  + 3⋅x + 10
  ⎝x  + 3⋅x + 10⎠  

>>> c = 0.25; (c, f.subs(x, c), df.subs(x, c))
(0.25, 0.0293801592482761, 0.134931846522447)

>>> c = 0.30; (c, f.subs(x, c), df.subs(x, c))
(0.3, 0.0362975936104176, 0.141747824460957)

>>> c = 0.60; (c, f.subs(x, c), df.subs(x, c))
(0.6, 0.0846090186079023, 0.179058168476784)

In [53]: print(sy.cxxcode(f))
std::exp(x)*std::sin(x)/(std::pow(x, 2) + 3*x + 10)

In [55]: print(sy.cxxcode(df))
(-2*x - 3)*std::exp(x)*std::sin(x)/std::pow(std::pow(x, 2) + 3*x + 10, 2) + std::exp(x)*std::sin(x)/(std::pow(x, 2) + 3*x + 10) + std::exp(x)*std::cos(x)/(std::pow(x, 2) + 3*x + 10)

cmake_minimum_required(VERSION 3.5)
project(Autodiff_experiments)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_VERBOSE_MAKEFILE ON)

#============================================================#

if(NOT EXISTS "${CMAKE_BINARY_DIR}/conan.cmake")
   message(STATUS "Downloading conan.cmake from https://github.com/conan-io/cmake-conan")
   file(DOWNLOAD "https://github.com/conan-io/cmake-conan/raw/v0.13/conan.cmake"
                 "${CMAKE_BINARY_DIR}/conan.cmake")
endif()

include(${CMAKE_BINARY_DIR}/conan.cmake)
set(CONAN_PROFILE default)

conan_cmake_run( REQUIRES
                 boost/1.71.0@conan/stable
                 BASIC_SETUP
                 BUILD missing )


#========== Target Configurations ================#

add_executable(main1 main1.cpp)

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>

#include <boost/math/differentiation/autodiff.hpp>

namespace bd = boost::math::differentiation;

template<typename T>
T func1(T x){
    return exp(x) * sin(x)/( pow(x, 2) + 3*x + 10);
}

// First derivate of function func1(x)
double func1_prime(double x)
{
    // Code generated by Python Sympy
    return (-2*x - 3)*std::exp(x)*std::sin(x) /
            std::pow(std::pow(x, 2) + 3*x + 10, 2)
            + std::exp(x)*std::sin(x)/(std::pow(x, 2) + 3*x + 10)
            + std::exp(x)*std::cos(x)/(std::pow(x, 2) + 3*x + 10);
}

int main()
{
    // Maximum order which the derivate will be computed
    constexpr unsigned Order = 5;
    auto xx = std::vector{0.25, 0.30, 0.60};
    std::cout << std::setprecision(10);

    std::cout << " ========== Experiment 1 ==========================" << std::endl;

    for(auto const& xc: xx)
    {
        auto x = bd::make_fvar<double, Order>(xc);
        auto y = func1(x);
        std::cout << "\n =>> For x = " << xc << std::endl;
        // 0-th derivate  => Same as f(x)
        std::cout << " [Auto] - f'[0](x) = " << y.derivative(0) << std::endl;
        // First derivate with autodiff
        std::cout << " [Auto] - f'[1](x) = " << y.derivative(1) << std::endl;
        // Firtst derivate with symbolic diff (Python Sympy)
        std::cout << " [Symb] - f'[1](x) = " << func1_prime(xc) << std::endl;
        // Second derivate
        std::cout << " [Auto] - f'[2](x) = " << y.derivative(2) << std::endl;
    }

    std::cout << "\n\n ==== Experiment 2 - Autodiff using lambda ====" << std::endl;

    // 'Auto variables' for lambdas works in a similar way to template type parameters
    auto fun1_lam = [](auto x){
        return exp(x) * sin(x)/( pow(x, 2) + 3*x + 10);
    };

    for(auto const& xc: xx)
      {
          auto x = bd::make_fvar<double, Order>(xc);
          auto y = fun1_lam(x);
          std::cout << "\n =>> For x = " << xc << std::endl;
          std::cout << " [Auto] - f'[0](x) = " << y.derivative(0)  << std::endl;
          std::cout << " [Auto] - f'[1](x) = " << y.derivative(1) << std::endl;
          std::cout << " [Symb] - f'[1](x) = " << func1_prime(xc)  << std::endl;
          std::cout << " [Auto] - f'[2](x) = " << y.derivative(2) << std::endl;
      }

    return 0;
}

========== Experiment 1 ==========================

=>> For x = 0.25
[Auto] - f'[0](x) = 0.02938015925
[Auto] - f'[1](x) = 0.1349318465
[Symb] - f'[1](x) = 0.1349318465
[Auto] - f'[2](x) = 0.1373348539

=>> For x = 0.3
[Auto] - f'[0](x) = 0.03629759361
[Auto] - f'[1](x) = 0.1417478245
[Symb] - f'[1](x) = 0.1417478245
[Auto] - f'[2](x) = 0.1352101205

=>> For x = 0.6
[Auto] - f'[0](x) = 0.08460901861
[Auto] - f'[1](x) = 0.1790581685
[Symb] - f'[1](x) = 0.1790581685
[Auto] - f'[2](x) = 0.1097378642


==== Experiment 2 - Autodiff using lambda ====

=>> For x = 0.25
[Auto] - f'[0](x) = 0.02938015925
[Auto] - f'[1](x) = 0.1349318465
[Symb] - f'[1](x) = 0.1349318465
[Auto] - f'[2](x) = 0.1373348539

=>> For x = 0.3
[Auto] - f'[0](x) = 0.03629759361
[Auto] - f'[1](x) = 0.1417478245
[Symb] - f'[1](x) = 0.1417478245
[Auto] - f'[2](x) = 0.1352101205

=>> For x = 0.6
[Auto] - f'[0](x) = 0.08460901861
[Auto] - f'[1](x) = 0.1790581685
[Symb] - f'[1](x) = 0.1790581685
[Auto] - f'[2](x) = 0.1097378642