pow_nan.cpp

View page source

pow: Nan in Result of Pow Function: Example and Test

Purpose

The pow(x, y) function will work when \(x < 0\) and \(y\) is a parameter. It will often generate nan or infinity when \(x < 0\) and one tries to compute a derivatives (even if \(y\) is a positive integer). This is because the derivative of the log is \(1 / x\) and the power function uses the representation

\[\R{pow}(x, y) = \exp [ y \cdot \log(x) ]\]

Problem

There is a problem with this representation when \(y\) is a parameter and \(x = 0\). For example, when \(x = 0\) and \(y = 1\), it returns zero for the derivative, but the actual derivative w.r.t \(x\) is one.

# include <cppad/cppad.hpp>
# include <cmath>

bool pow_nan(void)
{   bool ok = true;

    using std::cout;
    using CppAD::AD;
    using CppAD::vector;
    //
    vector<double>       x(2), y(2), dx(2), dy(2), ddx(2), ddy(2);
    vector< AD<double> > ax(2), ay(2);
    //
    // variable vector
    ax[0] = x[0]  = -1.0;
    ax[1] = x[1] = 2.0;
    //
    CppAD::Independent(ax);
    //
    ay[0] = pow( ax[0], ax[1] );
    ay[1] = pow( ax[0], 2.0 );
    //
    CppAD::ADFun<double> f(ax, ay);
    //
    // check_for_nan is set false so it does not generate an assert
    // when compiling with debugging
    f.check_for_nan(false);
    //
    // Zero order forward does not generate nan
    y  = f.Forward(0, x);
    ok &= y[0] == 1.0;
    ok &= y[1] == 1.0;
    //
    // First order forward generates a nan
    dx[0] = 1.0;
    dx[1] = 1.0;
    dy = f.Forward(1, dx);
    ok &= std::isnan( dy[0] );
    ok &= dy[1] == -2.0;
    //
    // Second order Taylor coefficient is 1/2 times second derivative
    ddx[0] = 0.0;
    ddx[1] = 0.0;
    ddy = f.Forward(2, ddx);
    ok &= std::isnan( ddy[0] );
    ok &= ddy[1] == 1.0;
    //
    return ok;
}