# CPP/C++ Math and Numerical Computing

## 1 Math and Numerical Computing

### 1.1 Numerical Libraries and Headers

All numeric libraries:

• <cmath> - C++ version of the C-header <math.h>. Contains basic trasncedental functions, sin, cos, tan and so on.
• C++11 Floating point classification functions:
• Widely used trasncedental functions: (C++11)
• std::erf - Error functions, can be used for computing normal function.
• std::erfc - Complementary error function.
• std::tgamma - Gamma function.
• Special math functions: (C++17)
• <complex> - Complex number library.
• <random> - High quality C++11 random generator library. Provides lots of distributions and random number engines.
• <limits> - Provides numeric limits of all C++ numeric types. For instance, minimum value, maximum value, number of digits, precision and episilon values of a given numeric type.
• <ratio> - Compile-time rational arithmetic library.
• <numeric> - Numerical "algorithms" or numerical functions for processing containers. std::accumulate, std::iner_product, std::partial_sum and std::iota.

### 1.2 Numerical Constants

Notes: Those constants are defined at the header <cmath>

Table 1: C++ Numerical constants in header <cmath>
Constant Value Description
General Math Constants
`M_E` 2.7182818 Euler's number or exp(1)
`M_LN2` 0.6931472 Natural logarithm of 2 or log(2)
`M_LN10` 2.3025851 Natural logarithm of 10 or log(10)

`M_LOG10E` 0.4342945 Log of E (Euler's number) at base 10
`M_LOG2E` 1.4426950 Log of E (Euler's number) at base 2

`M_SQRT2` 1.4142136 Square root of 2 or sqrt(2)
`M_SQRT1_2` 0.7071068 Square root of 1/2 or sqrt(1/2) or 1/sqrt(2)

`M_PI` 3.1415927 PI number
`M_PI_2` 1.5707963 PI/2
`M_PI_4` 0.7853982 PI/3
`M_1_PI` 0.3183099 1/PI
`M_2_PI` 0.6366197 2/PI
`M_2_SQRTPI` 1.1283792 2/sqrt(PI)

Example in CERN's Root REPL - Math Constants

```#include <cmath>

>> M_PI
(double) 3.1415927

>> sin(M_PI)
(double) 1.2246468e-16

>> cos(M_PI)
(double) -1.0000000

>> cos(M_PI_2)
(double) 6.1232340e-17

>> sin(M_PI_2)
(double) 1.0000000
>>
```

See:

### 1.3 Floating Point

#### 1.3.1 Overview

Computer Representation of Real Numbers

Main Representations:

Fixed Point

Most used in some financial calculations for handling monetary values or in embedded devices such as microcontrollers and some DSP - Digital Signal Processors or in legacy systems without a FPU Floating Point Unit where the cost of software emulation of floating points is not acceptable.

Types of Fixed Point:

• Binary Fixed Points
• Decimal Fixed Points
• Arbitrary Precision Fixed Points
• Arbitrary Precision Decimal Fixed Point (aka Big-Decimal)
• Better for handling monetary values such as prices.

Note:

• A systems without a FPU can emulate using software or it can use an external FPU chip.

Float Points

Most used in mainstream processors or systems with systems with FPU - (Floating Point Unit). Nowadays most of float point implementations are conformant to the IEEE754 standard.

• Binary Float Points (IEEE754)
• binary32 - (C++ floating) - single precision
• binary64 - (C++ double) - double precision
• binary128 - (C++ long double) - quadruple precision
• Decimal Floating Point (IEEE754)
• decimal32
• decimal64
• decimal128
• Arbitrary Precision Floating Points (Non IEEE754)

Technical Standards for Float Points

• IEEE 754:1985 - Standard for Binary Float-Point Arithmetic
• IEEE 754:1987 - Standard for Radix-Indenpendet Float-Point Arithmetic
• IEEE 754:2008 - Standard for Floating-Point Arithmetic
• ISO standard ISO/IEC/IEEE 60559:2011

Problems of Floating Points

• Floating Points - FP are not real numbers, they are an approximation or a computer representation for real numbers.
• Floating Points can only represent a finite and limited amount of real numbers.
• Not all real numbers can be exactly represented due to the quantiazation or discretization error.
• Operations are not commutative
• Loss of precision
• Cancellation aka catastrophic cancellation which is the total loss of precision when subtracting two large and very close real numbers.

Category of Floating Points Values

Any floating point values can fall in one of the following categories:

1. Signed Zeror (+0 or -0)
2. Normal numbers
3. Subnormals
4. Inifnities (positive and negative infinity)
5. NaN - Not a number.

Float Point Types in C++ X IEEE754

Name C++ Type Bits S E M B Description
binary16 N/A 16 1 5 10 15 IEE754 half precision floating point
binary32 float 32 1 8 23 127 IEE754 single precision floating point
binary64 double 64 1 11 52 1023 IEE754 double precision floating point
binary128 long double 128 1 15 112 16383 IEE754 quadruple precision floating point
- long double 80 1 15 64 16383 80 bit extended precision.

Symbols:

• S => Sign bits => 0 for positive numbers and 1 for negative numbers.
• E => Exponent bits
• B => Exponent Bias
• M => Mantiassa Bits (aka Fraction or significand)

Notes:

• long double => In some platforms, the type long double can be an 80-bits extended precision, not binary128. However, most computers and CPUs support binary128.
• binary16 (half precision floating point) => There is no C++ fundamental type with binary16 float point format. It mostly used in computer graphics and GPUs. The format was intruduced in NVidia GPUs and defined by the IEEE754-2008.

Exponent range:

Name   C++ Type ERmin ERmax Emin Emax
binar16 half precision -
binary32 single precision float 1 254 -126 127
binary64 double precision double 1 2026 -1022 1023
binary128 quadruple precision long double 1 32766 -16382 16383

Symbols:

• ERmin => Minimum value of raw expoent (without bias).
• ERMax => Maximum value of raw exponent (without bias)
• Emax = ERmax - bias => Maximum value of exponent.

Storage Layout

Float Point Type Bits Sign Exponent Mantissa
float 32 1  8 [23-30] 32 [0-22]
double 64 1  11 [52-62] 51 [0-51]
long double 128 1 
• 1  - Means that sign bit of the float point is the 31th bit.
• 8 [23-30] - Means that the exponents bits are are all bits between 23th and 30th bit (including both sides).

Binary Format for 32 bits floating points binary32, aka float

```  [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
```

Representation of Floating Points

• General Format:
```              Sign         Exponent - Bias
Value :=  (-1)     *  Radix        *        Mantissa

S      E - B
Value :=  (-1)  *  R    *  M
```
• Single Precision (floating)
• Where:
• S - sign bit - 0 for positive number, 1 for negative
• E - Exponent bits
• M - Mantissa bits
• Bias: 127
• b[i] i-th bit, b22 => bit 22.
```               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
```
• Double Precision (double)
```              S      (E - 1023)
Value :=  (-1)   X  2         X (1.M)
```

Special Values for Single Precision floating point number in Hexadecimal (floating type)

Case Case Bytes Values
Positive Zero +0.0 0x00000000
Negative Zero -0.0 0x80000000
Positive Infinity +INF 0x7F800000
Negative Infinity -INF 0xFF800000
Positive NaN +NaN 0x7FC00000
Negative NaN -NaN 0xFFC00000
Maximum positive 3.40E38 0x7F7FFFFF
Minimum negative -3.40E38 0xFF7FFFFF

Floating Point IEEE754 Parameters and Constants

Table 2: C++ Floating Point Parameters and Special Values in header <limits>
Constant Description

std::numeric_limits<T>::max() Maximum absolute value that a floating can represent.
std::numeric_limits<T>::min() Minimum absolute value that a floating can represent.
std::numeric_limits<T>::lowest() Minimum negative value that a floating can represent.

std::numeric_limits<T>::epsilon() Machine epsilon - maximum difference between two representable values.
std::numeric_limits<T>::quiet_NaN() NaN - Not a Number floating point value
std::numeric_limits<T>::signaling_NaN() NaN - Not a Number floating point value
std::numeric_limits<T>::infinity() inf - Maximum positive infinity
-std::numeric_limits<T>::infinity() -inf - Miminum negative infinity

Notes:

• eps (epsilon) is defined as the positive number x, such that x + 1.0 is not equal to 1.0. Two numbers X and Y, cannot be distinguished if ||X - Y|| < eps.
• The eps is the highest precision or the smallest tolerance tha any numerical algorithm or numerical method implementation can achieve.
• NaN is a special float error value indicating that the value cannot be represented. Any computation with NaN yields a NaN that is propagated to all other computations. This value is returned from operations such as dividing zero by zero, taking a square root of a negative number … and so on.

Examples in CERN's Root REPL - Float point parameters.

EPS (Epsilon) value for several floating point types.

```// 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
```

Testing float point equality machine precision EPS - epsilon.

``` // 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
>
```

Testing NaN constant - Not a Number

Alias values:

```// 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
```

NaN propagation: Any operation with a NaN yields a NaN.

```>> 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
```

Equality comparison between any floating point number and itself evaluates to true while equality between NaNs evaluates to false.

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

>> 10.0 == 10.0
(bool) true

>> 9.81451654 == 9.81451654
(bool) true

>> dnan == dnan
(bool) false
```

If-else statement:

```>> 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]

```

Testing floating point limits:

```>> 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
>>
```

#### 1.3.2 Equality

Equality based on absolute error:

• Two floating points are almost equal under a given tolerance when the relative error is smaller than some tolerance.
```abs(x - y)  < tolerance
```

Equality based on relative error (Better):

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

Or:

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

Or:

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

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

Equality based on relative error and asbsolute error:

• absTol => Absolute tolerance
• relTol => Relative tolerance
```def almostEqual(a, b, relTol, absTol):
return abs(x - y) < max(relTol * max(abs(x), abs(y)), absTol)
```

C++ Implementations

Equality based on absolute error:

• Note: It can be misleading for very small numbers.
```/** 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;
}
```

Equality based on relative error:

```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)) ;
}
```

Equality based on absolute and relative tolerance:

```/** 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);
}
```

Equality based on absolute and tolerance given as number of precision digits. 1 precision digit is 0.1 tolerance, 5 precision digits is a tolerance of 1E-5 or 0.00001.

```/** 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);
}
```

#### 1.3.3 Functions

1. Classify float point numbers

• int std::fclasify(float arg)
• int std::fclasify(double arg)
• int std::fclasify(long double arg)

See:

Probe function:

```#define CLASSIFY_FLOAT(expr) classify_float(#expr, (expr))

template<typename FLOAT>
void classify_float(const char* expr, FLOAT number)
{
int x = std::fpclassify(number);
if(x == FP_INFINITE)
{
std::cout << " Number x = " << expr << " is infinity" << "\n";
return;
}
if(x == FP_NAN)
{
std::cout << " Number x = " << expr  << " ; NAN not a number" << "\n";
return;
}
if(x == FP_NORMAL)
{
std::cout << " Number x = " << expr << " is normal" << std::endl;
return;
}
if(x == FP_SUBNORMAL)
{
std::cout << " Number x = " << expr << " is sub-normal" << std::endl;
return;
}
if(x == FP_ZERO)
{
std::cout << " Number x = " << expr << " is zero" << std::endl;
return;
}
std::cout << " [ERROR] Number x = " << expr << " cannot be classified. " << std::endl;
}
```

Testing code (main function):

```CLASSIFY_FLOAT(3451.0523);
CLASSIFY_FLOAT(1e100);
CLASSIFY_FLOAT(1e500);
CLASSIFY_FLOAT(1e-100);
CLASSIFY_FLOAT(1e-500);
CLASSIFY_FLOAT(0.0);
CLASSIFY_FLOAT(+0.0);
CLASSIFY_FLOAT(-0.0);
CLASSIFY_FLOAT(0.0 / 0.0);
CLASSIFY_FLOAT(1.0 / 0.0);
CLASSIFY_FLOAT(-1.0 / 0.0);
CLASSIFY_FLOAT( std::numeric_limits<double>::min());
CLASSIFY_FLOAT( std::numeric_limits<double>::min() / 10.0);
CLASSIFY_FLOAT( std::numeric_limits<double>::min() / 1000.0);
```

Output:

```Number x = 3451.0523 is normal
Number x = 1e100 is normal
Number x = 1e500 is infinity
Number x = 1e-100 is normal
Number x = 1e-500 is zero
Number x = 0.0 is zero
Number x = +0.0 is zero
Number x = -0.0 is zero
Number x = 0.0 / 0.0 ; NAN not a number
Number x = 1.0 / 0.0 is infinity
Number x = -1.0 / 0.0 is infinity
Number x = std::numeric_limits<double>::min() is normal
Number x = std::numeric_limits<double>::min() / 10.0 is sub-normal
Number x = std::numeric_limits<double>::min() / 1000.0 is sub-normal
```
2. Check if float is Finite, NaN or Infinity

Check if float point is NaN - Not-a-Number.

• std::isnan
```>> float x = 4.5;
>>
>> std::isnan(x)
(bool) false
>>
>> x = 4.0 / 0.0
(float) inff
>>
>> std::isnan(x)
(bool) false
>>
>> x = std::sqrt(-1.0)
(float) -nanf
>>
>> std::isnan(x)
(bool) true
>>
>> x = 0.0 / 0.0
(float) nanf
>>
>> std::isnan(x)
(bool) true
```

Check whether a floating point is infinity

• std::isinf
```>> std::isinf(1e100)
(bool) false

>> std::isinf(1e500)
(bool) true

>> std::isinf(1 / 0.0)
(bool) true

>> std::isinf(-1 / 0.0)
(bool) true

>> std::isinf(0.0 / 0.0)
(bool) false
```

Check whether floating point is finite

• std::isfinite
```>> double q = 10.0
(double) 10.000000
>> std::isfinite(q)
(bool) true

>> q = 0.0 / 0.0
(double) nan
>> std::isfinite(q)
(bool) false

>> q = 1.0 / 0.0
(double) inf
>> std::isfinite(q)
(bool) false

>> q = -1.0 / 0.0
(double) -inf
>> std::isfinite(q)
(bool) false
```
3. Maximum, minimum and absolute value

Maximum

```float       fmax( float x, float y );               // (1)  (since C++11)
double      fmax( double x, double y );             // (2)  (since C++11)
long double fmax( long double x, long double y );   // (3)  (since C++11)
```

Minimum

```float       fmin( float x, float y );              // (1)   (since C++11)
double      fmin( double x, double y );            // (2)   (since C++11)
long double fmin( long double x, long double y );  // (3)   (since C++11)
```

Absolute Value

• See: fabs (cppreference)

```float       fabs ( float arg );
float       fabsf( float arg );
double      fabs ( double arg );
long double fabs ( long double arg );
long double fabsl( long double arg );
```

4. Rounding and truncation

#### 1.3.4 Decompose 32 bits floating points

The following application decomposes 32 bits IEE754 floating point numbers showing the integer data representation in hexadecimal, sign bit, exponent and mantissa.

Source:

Compilation:

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

The results can be checked with the online tool:

Program Source Listing

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

Preprocessor Macros:

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

Function: show_float - shows decomposed floating point number.

```// 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("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";
}
```

Function main:

```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());
```

Program Output:

```\$ ./out.bin

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

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

Input Expression: +std::numeric_limits<float>::infinity()
Classification: FP_INFINITE
Decimal 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
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
Integer Representation (hex): 0x00800000
Sign bit: 0
Exponent: -126
Raw Exponent: 1
Raw Exponent (hex): 0x1
Mantissa (hex): 0x0

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

```

#### 1.3.5 Libraries

• Boost Multiprecision
• site: https://www.boost.org/doc/libs/1_71_0/libs/multiprecision/doc/html/boost_multiprecision/intro.html
• "The Multiprecision Library provides integer, rational, floating-point, and complex types in C++ that have more range and precision than C++'s ordinary built-in types. The big number types in Multiprecision can be used with a wide selection of basic mathematical operations, elementary transcendental functions as well as the functions in Boost.Math. The Multiprecision types can also interoperate with the built-in types in C++ using clearly defined conversion rules. This allows Boost.Multiprecision to be used for all kinds of mathematical calculations involving integer, rational and floating-point types requiring extended range and precision. Multiprecision consists of a generic interface to the mathematics of large numbers as well as a selection of big number back ends, with support for integer, rational, floating-point, and complex types. Boost.Multiprecision provides a selection of back ends provided off-the-rack in including interfaces to GMP, MPFR, MPIR, MPC, TomMath as well as its own collection of Boost-licensed, header-only back ends for integers, rationals and floats. In addition, user-defined back ends can be created and used with the interface of Multiprecision, provided the class implementation adheres to the necessary concepts."
• half - IEEE 754-based half-precision floating-point library
• site: http://half.sourceforge.net/index.html
• "This is a C++ header-only library to provide an IEEE 754 conformant 16-bit half-precision floating-point type along with corresponding arithmetic operators, type conversions and common mathematical functions. It aims for both efficiency and ease of use, trying to accurately mimic the behaviour of the built-in floating-point types at the best performance possible. It is hosted on SourceForge.net."

#### 1.3.6 Literature and Further Reading

Fundamental:

• FAQ What Every Computer Scientist Should Know About Floating-Point Arithmetic
• The Perils of Float Point
• A Tutorial on Data Representation Integers, Floating-point Numbers, and Characters
• An Introduction to Float-Point Arithmetic and Computation - CERN, Jeff Arnold

Float Point Comparison, Assertions and Testing

• PEP 485 – A Function for testing approximate equality [BEST] - Cristopher Baker - Noaa.gov
• https://www.python.org/dev/peps/pep-0485/
• "This PEP proposes the addition of an isclose() function to the standard library math module that determines whether one value is approximately equal or "close" to another value."
• Boost Libraries Floating Point Comparison [BEST]
• Comparing Float Point Numbers, 2021 Edition
• Comparing Floting Point Numbers - Bruce Dawson
• Comparing Floating-Point Numbers Is Tricky
• Boost - Float Point Comparison - Boost C++ Libraries
• Float Point Comparison in Rust

Floating Point Exceptions - SIGFPE:

Note - floating point exceptions are not C++ exceptions, they are signal.

Platform Specific Behavior

Floating point Binary16

Miscellaneous

• Real Close to the Machine: Floating Point in D
• Making floating point math highly efficient for AI hardware - Facebook
• Subnormal Floating-Point Degradation in Parallel Applications
• Float Point Numbers & Currency Rouding Errors
• Symbolic Execution for Checking Accuracy of Float-Point Programs (Paper)
• Five Tips for Floating Point Programs - John D. Coock
• Hydra Chronicles, Part III - Catastrophic Imprecision
• Number Representation and Errors - Jun Zhang
• Why Fixed Point Won't Cure Your Float Point Blues - Richard Harris
• Fixed-Point Math in C

### 1.4 C++11 Random Number Library

#### 1.4.1 Overview

C++11 provides a safer, more convenient and useful random number facilities than the old C-API rand() that fits many different domains such as numerical simulations, games, quizzes, password generation, security and so on.

• Predecessor: Boost.Random
• Advatanges over the old C-API rand():
• Safer and Stronger
• More user control
• Multiple random number engines
• Many probability distributions

Note:

• The C++11 random numbers' API is not object oriented, there is no class hierarchy. The library is generic and uses lots of C++'s generic programming (template meta programming) features. This design choice allows the library to be header-only and to be reused for many different types without virtual-function calls overhead.

Uses of Random Numbers

• Many Science, Physics and Engineering fields:
• Monte-Carlo Simulations
• Brownian Motion Simulation
• Simulation of noises and disturbances
• SDE
• Cryptography
• Generation of keys
• Games and RPG Role Play Games
• Shuffling
• Lottery

C++11 Random Number System:

The random number API has three parts, a seed generator, a pseudo-random number engine and a probability distribution.

• seed - Initial state of a pseudo-random number generator, a big random integer.
• seed generator: A non-deterministic random number generator that yields new seed random whenever it is invoked. Its purpose is to initialize the random engine.
• pseudo-random generator engine: Generates a reproducible sequence of big integers that are fed to a probability distribution. A pseudo-random generator always yields the same numerical sequence for the same seed or initial state, therefore, in order to make the generator really random, it is necessary to initialize it with a different seed during the program initialization. While the deterministic feature of the random engines may look like a disadvantange, they are useful for testing and reproducibility of results, simulations and experiments. For instance, by using the same seed and engine, it is possible to reproduce and test the results of monte-carlo numerical simulation.
• Note: The engines are deterministic, but can be made non-deterministic by changing the seed during the initialization.
• probability distribution: A probability distribution object takes the generator output, a big number, and returns a random number specified by the distribution.

Some Engines:

Library Class Description
Seed Generator
std::random_device Generates a seed, aka non-deterministic random number.

Random Generator Engines
std::default_random_engine Default engine, the algorithm used is implementation defined.
std::mt19937 32-bit Mersenne Twister by Matsumoto and Nishimura, 1998
std::mt19937_64 64-bit Mersenne Twister by Matsumoto and Nishimura, 2000
std::minstd_rand0 Linear Congruential
std::minstd_rand Linear Congruential
std::ranlux48
std::knuth_b

Probability Distributions
std::uniform_int_distribution Generates uniformly distributed integers in the interval [a, b]
std::uniform_real_distribution<T> Generates uniformly distributed float point numbers.

std::normal_distribution Yields normally distributed numbers with some mean and standard deviation.
std::lognormal_distribution Variant of the normal distribution.
-
bernouli_distributions<T> -
binomial_distribution<T> -
negative_binonmial_distribution<T> -
geometric_distribution<T> -

Reproducibility in Science and Research

#### 1.4.2 Basic Examples

Example in CERN's ROOT REPL:

• Example 1: The function generate_randoms1 fails to generate a random sequence as whenever the function is called, it generates the same sequence of numbers because the seed was not set.
```#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;
}
```

Testing:

```>> generate_randoms1(4)
x = 1
x = 2
x = 8
x = 5

>> generate_randoms1(4)
x = 1
x = 2
x = 8
x = 5

>> generate_randoms1(6)
x = 1
x = 2
x = 8
x = 5
x = 6
x = 3

>> generate_randoms1(10)
x = 1
x = 2
x = 8
x = 5
x = 6
x = 3
x = 1
x = 7
x = 7
x = 10
```
• Example 2: The function generate_randoms2, generates quasi-random numbers whenever the seed is changed, if it is kept unchanged, the distribution and engine will always generate the same results.
```#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;
}
```

Testing:

```// 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 = 4
x = 1
x = 5
x = 10
x = 6
>>
>> generate_randoms2(6, seed)
x = 4
x = 1
x = 5
x = 10
x = 6
x = 10

>> generate_randoms2(10, seed)
x = 4
x = 1
x = 5
x = 10
x = 6
x = 10
x = 2
x = 8
x = 9
x = 6
>>
```

Changing the seed:

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

>> generate_randoms2(4, seed2)
x = 8
x = 5
x = 2
x = 7

>> generate_randoms2(8, seed2)
x = 8
x = 5
x = 2
x = 7
x = 6
x = 9
x = 4
x = 1
```

Making the sequence "really" random:

```>> generate_randoms2(5, seed_gen())
x = 1
x = 9
x = 8
x = 10
x = 2

>> generate_randoms2(5, seed_gen())
x = 5
x = 6
x = 7
x = 8
x = 4

>> generate_randoms2(5, seed_gen())
x = 3
x = 4
x = 1
x = 8
x = 8
```

#### 1.4.3 Example - Normally Distributed Randoms

The following sample code generates N normally distributed random numbers with mean 10.0 and standard deviation 5 in three different ways: first with seed not set; then with seed set by used and third with seed generated automatically.

• File: randoms1.cpp
```#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);
std::string command = argv;

// 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);
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;
}

```

Compiling:

```# 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
```
• Experiment 1: Running command empty => generates N normally distributed numbers with mean 4 and standard deviation 10 without setting the engine rngen's seed, aka initial state.
```# 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
```
• Experiment 2: Running command auto => generate N random numbers with same features, but with a seed generated automatically by the object seed_gen with type random_device.

Simulation run A:

```\$ ./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
```

Simulation run B:

```\$ ./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
```

Simulation run C:

```\$ ./randoms1.bin 3 -auto
seed = 1536173647
Next random is = 13.4482
Next random is = 12.2431
Next random is = 12.2183
```
• Experiment 3: The results of the simulation runs of experiment 2 can be reproduced by reusing the seed with command -seed.

Reproduce simulation run A:

```\$ ./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
```

Reproduce simulation run B:

```\$ ./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
```

#### 1.4.4 Example - Class encapsulating uniform random distribution

Source:

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

Class RndInt => Class for generating random numbers between an integer interval including both bounds in a simplified way similar to the old C's rand() API:

```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);
}
};
```

Main function:

```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;
```

Compilation:

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

Sample program outpu:

```\$ 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 =     9
x =     1
x =    10
x =     6
x =    10
x =     5
x =     2
x =     9
x =     9
x =    10
x =     5
x =     9
x =     8
x =     5
x =     9

===== Random numbers with a non-random seed ====
=> rnd.Seed() = 1195785783
=> rnd.Min() = 1
=> rnd.Max() = 10
x =     7
x =    10
x =     9
x =     2
x =     2
x =     1
x =     3
x =     1
x =     6
x =     5
x =     1
x =     2
x =     1
x =     3
x =     2

```

#### 1.4.5 Example - Functional wrapper for random uniform distribution

The function make_uniform_random_distribution returns a lambda object (anonymous callable object) that generates uniformly distributed randoms with a given seed.

Sources:

```#include <iostream>
#include <random>
#include <iomanip>
#include <string>
```
• Function overloading 1 of make_uniform_random_distribution - Returns a fuction (lambda object) that generates random integers with a seed supplied by the calling code.
```// 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 of make_uniform_random_distribution - Returns a "fuction" (lambda object) that generates random integers with a random seed.
```// 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()) ;
}
```
• Main Function:
```#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;
}
}
```

Sample Output:

Compiling:

```# 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
```

Running: (64 bits processor)

```\$ ./rand_int_fun.bin

===== Random numbers with a random seed ====
x =     8
x =     4
x =     4
x =     7
x =     9
x =     2
x =     2
x =     6
x =     3
x =     8
x =     3
x =     1
x =     3
x =     1
x =    10

===== Random numbers with a non-random seed ====
x =     6
x =    10
x =     6
x =     2
x =     4
x =     1
x =     9
x =    10
x =     3
x =     6
x =    10
x =     6
x =     2
x =     1
x =     1

```

Second for-loop output in 32 bits processor (x86) or with 32 bits compiler:

```===== Random numbers with a non-random seed ====
x =     3
x =     6
x =     1
x =     5
x =     4
x =     7
x =     2
x =     6
x =     7
x =     4
x =     8
x =     7
x =     2
x =     1
x =    10
```

### 1.5 Numerical Libraries Evaluation

#### 1.5.1 Boost Autodiff - Automatic Differentiation

Automatic differentiation is a third approach for computing exact derivatives, gradient, Jacobian matrix and Hessian matrix without any truncation error, roundoff error or loss of precision.

Automatic differentiation, aka algorithm differentiation is not:

• Symbolic differentiation => No symbolic manipulation is performed.
• Numeric / finite difference differentiation

• Faster than finite difference approximation
• No numerically unstable such as finite difference approximations.
• Avoids computing derivatives manually
• Avoids loss of precision, truncation errors and possible catastrophic cancellation.

Use-cases:

• Many numerical algorithms such as Newton-Raphson method for root finding; numerical optimization algorithms and simulations require computing derivatives.
• Solving nonlinear equations
• Optimization
• Newton Raphson method
• BFGS method (Broyden–Fletcher–Goldfarb–Shanno algorithm)
• Computer Simulations
• Computing sensitivities

Generating a test function and its derivative

Generate a function and its symbolic derivate in Python Sympy:

• Define a function
```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
```
• Compute symbolic derivate (first derivative):
```>>> 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⎠
```
• Compute value of the function and its derivative in a set of pre-defined points:
```>>> 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)
```
• Generate C++ code for the test function
```In : print(sy.cxxcode(f))
std::exp(x)*std::sin(x)/(std::pow(x, 2) + 3*x + 10)
```
• Generate C++ code for the derivate of the test function
```In : 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)
```

Sample Project with Boost Autodiff Library

File: CMakeLists.txt (Requires Conan)

```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")
"\${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 ================#

```

File: 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'(x) = " << y.derivative(0) << std::endl;
// First derivate with autodiff
std::cout << " [Auto] - f'(x) = " << y.derivative(1) << std::endl;
// Firtst derivate with symbolic diff (Python Sympy)
std::cout << " [Symb] - f'(x) = " << func1_prime(xc) << std::endl;
// Second derivate
std::cout << " [Auto] - f'(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'(x) = " << y.derivative(0)  << std::endl;
std::cout << " [Auto] - f'(x) = " << y.derivative(1) << std::endl;
std::cout << " [Symb] - f'(x) = " << func1_prime(xc)  << std::endl;
std::cout << " [Auto] - f'(x) = " << y.derivative(2) << std::endl;
}

return 0;
}
```

Program output:

```========== Experiment 1 ==========================

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

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

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

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

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

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

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

#### 1.5.2 Eigen - Linear Algebra

1. Overview

Eigen is a high-performance and header-only linear algebra library suitable for a wide variety of applications, including, numerical linear algebra, numerical statistics and numerical methods.

• Many types of linear algebra operations:
• scalar X vector operations
• scalar X matrix operations
• vector X matrix operations
• matrix X matrix operations
• solution of linear systems
• norm of vectors and norm of matrices
• eigenvalues and eigenvectors
• Matrix decomposition: QR, SVD, LU, Cholesky and so on.

Applications:

• Science
• Physics
• Engineering
• Implementation of numerical methods which requires linear algebra.
• Multivariate statistics
• PDE - Partial Differential Equations
• FEM - Finite Elements Methods
• Computer graphics

Perofmance: Eigen uses the following techniques for achieving high performance.

• Expression templates and lazy evaluation for optmizing loops. For instance, with just operator overloadings, the sum of three vectors would require three loops. Due the usage of expression templates, the sum of three vectors can be computed in a single loop.
• CRTP Curious Recurring Template for virtual member-function call overhead, specially in loops.
• SIMD - Single Instruction Multiple Data (Data Parallelism) - for speeding up linear algebra routines. SIMD instructions allows processing several values with a single instruction.
• Classes for sparse matrices which avoids storing all non-zero elements.
• Specific matrices and vector classes for computer graphics of 2, 3 and 4 dimensions.
• Fixed-size matrices without dynamic meory allocation.

Some high-profile projects using Eigen:

• Google's TensorFlow => Library for Machine Learning
• Shogun => a large scale machine learning toolbox.
• Gnu Data Language, a GPL interpretor of IDL syntax codes.

Official Repository

Conan Recipe

Documentation / Doxygen

Miscellaneous

2. Sample CMake Project

Project Files

File: CMakeLists.txt / Version 1 - without Conan

• This sample CmakeLists.txt requires any dependencies and no previous installation of the Eigen library. It has a macro for automating the library.
```cmake_minimum_required(VERSION 3.5)
project(Eigen_test)

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

include(FetchContent)
IF(\${is_via_git})
FetchContent_Declare(
eigenlib
GIT_REPOSITORY  https://gitlab.com/libeigen/eigen
GIT_TAG         \${version}
)
ELSE()
FetchContent_Declare(
eigenlib
URL https://gitlab.com/libeigen/eigen/-/archive/\${version}/eigen-\${version}.zip
)
ENDIF()
FetchContent_GetProperties(eigenlib)
if(NOT eigenlib_POPULATED)
FetchContent_Populate(eigenlib)
endif()
include_directories(\${CMAKE_BINARY_DIR}/_deps/eigenlib-src )
endmacro()

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

File: CMakeLists.txt / Version 2 - using Conan package manager.

• The advantage of this CMakeLists.txt file is that the library is donwload only once saving disk space, while the first version of the CMakeLists.txt downloads the entire library, about 160 MB.
```cmake_minimum_required(VERSION 3.5)
project(Eigen_test)

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

# ============= Conan Bootstrap =============================#

if(NOT EXISTS "\${CMAKE_BINARY_DIR}/conan.cmake")
"\${CMAKE_BINARY_DIR}/conan.cmake")
endif()

include(\${CMAKE_BINARY_DIR}/conan.cmake)

conan_cmake_run(REQUIRES
eigen/3.3.7@conan/stable

BASIC_SETUP
BUILD missing
)

# find_package (Eigen3 3.3 REQUIRED NO_MODULE)

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

File: main1.cpp

```#include <iostream>
#include <iomanip>
#include <Eigen/Dense>

using Eigen::MatrixXd;

int main()
{
std::puts("\n >>>[EXPERIMENT 1] ========= Create Matrix ================");

// Matrix of variable size.
MatrixXd m(3,3);

std::cout << " Matrix features "
<< "=> rows = " << m.rows()
<< "; cols = " << m.cols()
<< std::endl;

std::cout << "\n  m = \n" << m << std::endl;

std::puts("\n >>>[EXPERIMENT 3] ====== Initialize Matrix  ================");
{

m(0,0) = 3;
m(1,0) = 2.5;
m(0,1) = -1;
m(0,2) = 9.23;

m(1,1) = m(1,0) + m(0,1);
m(1,2) = 10.523;

m(2, 0) = 7.24;
m(2, 2) = 4.56;
m(2, 1) = 2.5;
std::cout << " m = \n" << m << std::endl;

std::cout << "\n Transpose matrix =>>\n m^T = \n" << m.transpose() << std::endl;

std::cout << "\n Inverse matrix =>>\n m^-1 = \n" << m.inverse() << std::endl;
}

std::puts("\n >>>[EXPERIMENT 4] === Create Initialized Matrix =========");
{
MatrixXd A(3, 4);
A <<  10.255, -21.5, 5.15, 9.5
,9.251,   18.4, 7.25, 25.1
,56.13,   24.5, 8.25, 17.55;

std::cout << " Matrix A  = \n" << A << std::endl;
}

std::puts("\n >>>[EXPERIMENT 5] === Multiply matrix by scalar =========");
{
std::cout << "\n 10 * m = \n" << 10.0 * m << std::endl;
std::cout << "\n m  / 5 = \n" << m / 5 << std::endl;
}

std::puts("\n >>>[EXPERIMENT 6] === Matrix-Matrix operations ==========");
{
// Create a 3 x 3 random matrix multiplied by 10
auto rndm1 =  10 * MatrixXd::Random(3, 3);

std::cout << "\n rndm1 = \n" << rndm1 << std::endl;
std::cout << "\n m + rndm1 = \n" << m + rndm1 << std::endl;

std::cout << "\n m * rndm1 = \n" << m * rndm1 << std::endl;
std::cout << "\n rndm1 * m = \n" << m * rndm1 << std::endl;
}

std::puts("\n >>>[EXPERIMENT 7] === Eigenvalues and Eigenvectors ==========");
{
MatrixXd A(3, 3);

A <<  2.50000E-02, 3.25000E-02,   1.75000E-03
,3.25000E-02, 2.17000E-01,  -2.15000E-03
,1.75000E-03, -2.15000E-03,  4.30000E-04;

std::cout << std::fixed << std::setprecision(5);

std::cout << " A = \n" << A << std::endl;

assert(eigsolver.info() == Eigen::Success
&& "[ERROR] Unable to find eigenvalues");

std::cout << "\n Eigenvalues = \n"
<< eigsolver.eigenvalues() << std::endl;

std::cout << "\n Eigenvectiors = \n"
<< eigsolver.eigenvectors() << std::endl;
}

}
```

Program Output

``` >>>[EXPERIMENT 1] ========= Create Matrix ================
Matrix features => rows = 3; cols = 3

m =
0 0 0
0 0 0
0 0 0

>>>[EXPERIMENT 3] ====== Initialize Matrix  ================
m =
3     -1   9.23
2.5    1.5 10.523
7.24    2.5   4.56

Transpose matrix =>>
m^T =
3    2.5   7.24
-1    1.5    2.5
9.23 10.523   4.56

Inverse matrix =>>
m^-1 =
0.117459 -0.166738  0.147026
-0.390894  0.320655 0.0512492
0.0278148 0.0889348 -0.042235

>>>[EXPERIMENT 4] === Create Initialized Matrix =========
Matrix A  =
10.255  -21.5   5.15    9.5
9.251   18.4   7.25   25.1
56.13   24.5   8.25  17.55

>>>[EXPERIMENT 5] === Multiply matrix by scalar =========

10 * m =
30    -10   92.3
25     15 105.23
72.4     25   45.6

m  / 5 =
0.6   -0.2  1.846
0.5    0.3 2.1046
1.448    0.5  0.912

>>>[EXPERIMENT 6] === Matrix-Matrix operations ==========

rndm1 =
6.80375   5.9688 -3.29554
-2.11234  8.23295  5.36459
5.66198 -6.04897 -4.44451

m + rndm1 =
4.0794 -3.70431  17.5539
2.04794  1.76802  13.2372
9.81742  11.5446  8.90594

m * rndm1 =
-112.934  47.9796 -86.9588
-116.51  40.2783  -98.052
-90.6609 -27.6275 -88.4288

rndm1 * m =
-85.4597    14.787  -10.5057
-63.8875   34.5262 -0.960369
-57.3932    29.101   -20.442

>>>[EXPERIMENT 7] === Eigenvalues and Eigenvectors ==========
A =
0.02500  0.03250  0.00175
0.03250  0.21700 -0.00215
0.00175 -0.00215  0.00043

Eigenvalues =
0.00019
0.01987
0.22237

Eigenvectiors =
0.10336 -0.98130  0.16240
-0.02535  0.16062  0.98669
-0.99432 -0.10610 -0.00828

```

Project in QtCreator IDE

Just open the CMakeLists.txt file or run the following comamnd from command line:

```\$ qtcreator /path/to/project/CMakeLists.txt &
```

Project in CodeBlocks IDE

Generate CodeBlocks project:

```\$ cmake -B_build -H. -G "CodeBlocks - Unix Makefiles" -DCMAKE_BUILD_TYPE=Debug
-- The C compiler identification is GNU 8.3.1
-- The CXX compiler identification is GNU 8.3.1
-- Check for working C compiler: /usr/lib64/ccache/cc
-- Check for working C compiler: /usr/lib64/ccache/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /usr/lib64/ccache/c++
-- Check for working CXX compiler: /usr/lib64/ccache/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
....  ....  ....  ....  ....  ....  ....  ....  ....
....  ....  ....  ....  ....  ....  ....  ....  ....
```

List generated project files:

```\$ ls _build/
CMakeFiles/  _deps/  CMakeCache.txt  cmake_install.cmake  Eigen_test.cbp  Makefile

\$ file _build/Eigen_test.cbp
_build/Eigen_test.cbp: XML 1.0 document, ASCII text

\$ ls _build/_deps/eigenlib-src/
bench/  doc/       test/           COPYING.LGPL       CTestCustom.cmake.in
blas/   Eigen/     unsupported/    COPYING.MINPACK    eigen3.pc.in
cmake/  failtest/  CMakeLists.txt  COPYING.MPL2       INSTALL
demos/  scripts/   COPYING.GPL     CTestConfig.cmake  signature_of_eigen3_matrix_library
```

Show CodeBlocks project file:

```\$ head -n 10 _build/Eigen_test.cbp

<?xml version="1.0" encoding="UTF-8"?>
<CodeBlocks_project_file>
<FileVersion major="1" minor="6"/>
<Project>
<Option title="Eigen_test"/>
<Option makefile_is_custom="1"/>
<Option compiler="gcc"/>
<Option virtualFolders="CMake Files\;"/>
<Build>
<Target title="all">
```

Open the following file in CodeBlocks IDE

• '/path/to/project-root-dir/_build/Eigen_test.cbp'

Compiling from command line

Configuration step:

```\$ cmake -B_build -H. -DCMAKE_BUILD_TYPE=Debug
-- The C compiler identification is GNU 8.3.1
-- The CXX compiler identification is GNU 8.3.1
-- Check for working C compiler: /usr/lib64/ccache/cc
-- Check for working C compiler: /usr/lib64/ccache/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /usr/lib64/ccache/c++
-- Check for working CXX compiler: /usr/lib64/ccache/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /home/user/projects/eigen/_build

```

Building step:

```\$ cmake --build _build --target
.... .... .... .... ... ... ... ...
/usr/lib64/ccache/c++  -g   CMakeFiles/main1.dir/main1.cpp.o  -o main1
gmake: Leaving directory '/home/user/projects/eigen/_build'
[100%] Built target main1
gmake: Leaving directory '/home/user/projects/eigen/_build'
/usr/bin/cmake -E cmake_progress_start /home/user/projects/eigen/_build/CMakeFiles
```

Running application:

```\$ _build/main1

>>>[EXPERIMENT 1] ========= Create Matrix ================
Matrix features => rows = 3; cols = 3

... .... ... ... .... ... ... .... ... ... .... ... ... .... ...
```

#### 1.5.3 Armadillo - Linear Algebra

1. Overview

Armadillo is high performance numerical linear algebra library built on top of battle tested BLAS, Lapack and OpenBLAS libraries. The library provides an easy-to-use interface and functionality similar to ©Matlab and Octave (Matlab open source clone).

Benefits

• Expression templates for compile-time optimization.
• BLAS and LAPACK backend which are the de-factor standard for numerical linear algebra and scientific computing.
• Support for OpenMP
• Classes for:
• Vectors, column vectors
• Matrices
• Sparces matrices
• Tensors (multi-dimensional arrays)
• Functionalities similar to Matlab and Octave

Drawbacks

• Requires BLAS, OpenBLAS, SuperLU and arpack binary depedencies. Some, of those depedencies are written in Fortran which makes harder to build them from sources on Windows. Altough, it is not hard to install Armadillo via Linux distributions' package managers.

Official Web Site

Documentation

• General Documentation:
• Matlab/Octave translation table:

• Armadillo: An Open Source C++ Linear AlgebraLibrary for Fast Prototyping and Computationally Intensive Experiments
• Matrix Computing in C++ADVANCED PROGRAMMING IN SCIENCE AND TECHNOLOGY 2012
• RcppArmadillo: Accelerating Rwith High-Performance C++ Linear Algebra1
• Seamless R and C++Integration with Rcpp:Part 2 – RcppArmadillo Examples - Dirk Eddelbuettel
• FindArmadillo — CMake 3.17.0-rc1 Documentation
2. Sample CMake Project

Armadillo can be installed by building from sources. However, it is not so easy due to the BLAS and LAPACK binary dependencies. Another easier way to install the library by using Linux distributions' package managers.

Install on Fedora Linux:

```\$ sudo dnf install armadillo-devel.x86_64
```

Install on Ubuntu or Debian:

```\$ apt-get install libarmadillo-dev
```

Project Files

File: CMakeListst.txt

```make_minimum_required(VERSION 3.5)

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

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

```#include <iostream>
#include <iomanip>

#include <cmath>
#include <cassert>

using arma::mat;
using arma::randu;

int main()
{
std::cout << " [INFO] Armadillo version = " << arma::arma_version::as_string()
<< "\n\n";

std::cout << "======== EXPERIMENT 1 === Create Random Matrices ==========\n\n";
{

// Random 3 x 4 matrix
mat matrixA = randu<mat>(3, 4);
std::cout << " =>    Number of rows of matrixA = " << matrixA.n_rows << "\n";
std::cout << " => Number of columns of matrixA = " << matrixA.n_cols << "\n";

//std::cout << " Dimensions of matrix A: rows = " << matrixA.rows()  << std::endl;

// Random 4 x 5 matrix
auto matrixB = randu<mat>(4, 5);
std::cout << "matrixB = \n" << matrixA << std::endl;

auto matrixC = matrixA * matrixB;
matrixC.print(" C = A * B = ");
}

std::cout << "======== EXPERIMENT 2 === Special Matrices ====================\n\n";
{
auto identity3 = arma::mat(3, 3, arma::fill::eye);

std::cout << " Identity(3x3) = \n"  << identity3;

auto zeromatrix = arma::mat(4, 5);
zeromatrix.fill(0.0);
std::cout << " ZeroMatrix = \n" << zeromatrix << "\n";

auto zeromatrix2 = arma::zeros<arma::mat>(3, 3);
std::cout << " ZeroMatrix2 = \n" << zeromatrix2 << "\n";

}

std::cout << "======== EXPERIMENT 3 === Initialize matrix manually ==========\n\n";
{
// Note: matrices are zero-based not one-based like Fortran or Matlab
mat A(3, 3);
A(0, 0) = 10.0;
A(0, 1) = 5.6;
A(0, 2) = -8.24;
A(1, 0) = 7.251;
A(1, 1) = 18.62;
A(1, 2) = 20.51;
A(2, 0) = -8.0;
A(2, 1) = 1.351E-3;
A(2, 2) = 100.0;

mat B;
B << 2.541 <<  9.251 << 7.51  << arma::endr
<< 9.871 << -12.56 << 81.25 << arma::endr
<< 10.52 <<  72.15 << -67.23 << arma::endr;

std::cout << " Determinant of A = " << arma::det(A) << "\n";
std::cout << " Determinant of B = " << arma::det(B) << "\n";
std::cout << " Sum of elements of A => SUM(A) = " << arma::accu(A) << "\n";
std::cout << " Sum of elements of B => SUM(B) = " << arma::accu(B) << "\n";
std::cout << " Trace of elements of A => TRACE(A) = " << arma::trace(A) << "\n";

A.print("\n matrix A = ");

B.print("\n matrix B = ");

// Tranpose of matrix A
std::cout << "\n transpose(A) = A^t = \n" << A.t() << "\n";

// Inverse of matrix A
auto Ainv = arma::inv(A);
Ainv.print("\n  A^-1 = inv(A) = ");

std::cout << "\n  A^-1 = inv(A) = \n" << A.i() << "\n";

// Linear comnbination of A and B
auto C = 3.0 * A + 4.0 * B;
C.print("\n C = 3A + 4B = ");
}

std::cout << "\n======== EXPERIMENT 4 === Solving Linear Systems ==========\n\n";
{
auto b = arma::vec{51.0,  61.0, 21.0};
b.print(" b = ");

arma::mat A{  {3.0,  5.0, 6.0}
, {8.0, 10.0, 3.0}
, {4.0, 1.0, 2.0}};

A.print(" A = ");

// Solve linear equation: A.x = b
// or perform the operation x = A \ b (Matlab notation)
// Expects results [2.0 3.0 5.0]'
arma::vec x = arma::solve(A, b);
x.print(" x = ");

std::cout << "normL2(x) = euclidianNorm(x) " << arma::norm(x) << "\n";

auto x_expected = arma::vec{2.0, 3.0, 5.0};
assert(arma::norm(x_expected - x) < 1e-3);

std::cout << " [INFO] Tests passed Ok." << "\n";
}
}
```

Program output:

``` [INFO] Armadillo version = 9.400.2 (Surrogate Miscreant)

======== EXPERIMENT 1 === Create Random Matrices ==========

=>    Number of rows of matrixA = 3
=> Number of columns of matrixA = 4
matrixB =
0.7868   0.9467   0.2513   0.3447
0.2505   0.0193   0.0227   0.2742
0.7107   0.4049   0.5206   0.5610

C = A * B =
1.0516   1.0632   0.7769   0.8944   0.8861
0.2924   0.2185   0.1320   0.3470   0.2753
1.0723   1.0523   0.8662   1.1705   0.9633
======== EXPERIMENT 2 === Special Matrices ====================

Identity(3x3) =
1.0000        0        0
0   1.0000        0
0        0   1.0000
ZeroMatrix =
0        0        0        0        0
0        0        0        0        0
0        0        0        0        0
0        0        0        0        0

ZeroMatrix2 =
0        0        0
0        0        0
0        0        0

======== EXPERIMENT 3 === Initialize matrix manually ==========

Determinant of A = 12412.8
Determinant of B = 7637.21
Sum of elements of A => SUM(A) = 145.742
Sum of elements of B => SUM(B) = 113.303
Trace of elements of A => TRACE(A) = 128.62

matrix A =
1.0000e+01   5.6000e+00  -8.2400e+00
7.2510e+00   1.8620e+01   2.0510e+01
-8.0000e+00   1.3510e-03   1.0000e+02

matrix B =
2.5410    9.2510    7.5100
9.8710  -12.5600   81.2500
10.5200   72.1500  -67.2300

transpose(A) = A^t =
1.0000e+01   7.2510e+00  -8.0000e+00
5.6000e+00   1.8620e+01   1.3510e-03
-8.2400e+00   2.0510e+01   1.0000e+02

 A^-1 = inv(A) =
0.1500  -0.0451   0.0216
-0.0716   0.0753  -0.0213
0.0120  -0.0036   0.0117

 A^-1 = inv(A) =
0.1500  -0.0451   0.0216
-0.0716   0.0753  -0.0213
0.0120  -0.0036   0.0117

C = 3A + 4B =
4.0164e+01   5.3804e+01   5.3200e+00
6.1237e+01   5.6200e+00   3.8653e+02
1.8080e+01   2.8860e+02   3.1080e+01

======== EXPERIMENT 4 === Solving Linear Systems ==========

b =
51.0000
61.0000
21.0000
A =
3.0000    5.0000    6.0000
8.0000   10.0000    3.0000
4.0000    1.0000    2.0000
x =
2.0000
3.0000
5.0000
normL2(x) = euclidianNorm(x) 6.16441
[INFO] Tests passed Ok.
```

Created: 2021-06-04 Fri 15:07

Validate