NAG Library Manual, Mark 27.3
```/* nag::ad::e01bf Adjoint Example Program.
*/

#include <dco.hpp>
#include <iostream>

// Function which calls NAG AD routines.
template <typename T> void func(const T &xint, T &fint);

// Driver with the adjoint calls.
//  Evaluates a piecewise cubic Hermite interpolant and computes function value
//  at point xint.
// Also computes the derivative d(fint)/d(xint).
void driver(const double &xint, double &fint, double &d);

int main()
{

// Point to interpolate at
double xint = 8.4;
// Interpolated function value
double fint;

// Derivative
double d;

// Call driver
driver(xint, fint, d);

// Print outputs
std::cout << " Derivatives calculated: First order adjoints\n";
std::cout << " Computational mode    : algorithmic\n";

std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(6);
std::cout << "\n fint = " << fint << std::endl;

// Print derivative
std::cout << "\n Derivative of fint w.r.t. interpolation point xint :\n";
std::cout << " d(fint) / d(xint) = " << d << std::endl;

return 0;
}

// Driver with the adjoint calls.
//  Evaluates a piecewise cubic Hermite interpolant and computes function value
//  at point xint.
// Also computes the derivative d(fint)/d(xint).
void driver(const double &xintv, double &fintv, double &d)
{
using mode = dco::ga1s<double>;
using T    = mode::type;

mode::global_tape = mode::tape_t::create();

// Variable to differentiate w.r.t.
T xint = xintv;

// Register variable to differentiate w.r.t.
mode::global_tape->register_variable(xint);

// Interpolated function value
T fint;

// Call the NAG AD Lib functions
func(xint, fint);

// Extract the computed solutions
fintv = dco::passive_value(fint);
mode::global_tape->register_output_variable(fint);
dco::derivative(fint) = 1.0;

// Interpret the tape

d = dco::derivative(xint);

// Remove tape
mode::tape_t::remove(mode::global_tape);
}

// Function which calls NAG AD Library routines
template <typename T> void func(const T &xint, T &fint)
{

// Create AD configuration data object
Integer           ifail = 0;

// Initialize abscissa and ordinate data
std::vector<T> x{7.99, 8.09, 8.19, 8.70, 9.20, 10.00, 12.00, 15.00, 20.00};
std::vector<T> f{0.00000E+0, 0.27643E-4, 0.43750E-1, 0.16918E+0, 0.46943E+0,
0.94374E+0, 0.99864E+0, 0.99992E+0, 0.99999E+0};

// Local output of e01be; approximation of derivatives at points x
std::vector<T> d(x.size());
ifail = 0;