master: pow.

pow.sh: demonstration of unexpected results.
pow_op.hpp: fix case where computing higher order reverse.
pow.cpp: test higher order reverse.

bug/pow.sh

@@ -0,0 +1,85 @@
+#! /bin/bash -e
+# vim: set expandtab:
+# -----------------------------------------------------------------------------
+# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-18 Bradley M. Bell
+#
+# CppAD is distributed under the terms of the
+#              Eclipse Public License Version 2.0.
+#
+# This Source Code may also be made available under the following
+# Secondary License when the conditions for such availability set forth
+# in the Eclipse Public License, Version 2.0 are satisfied:
+#       GNU General Public License, Version 2.0 or later.
+# -----------------------------------------------------------------------------
+name=`echo $0 | sed -e 's|^bug/||' -e 's|\.sh$||'`
+if [ "$0" != "bug/$name.sh" ]
+then
+    echo 'usage: bug/pow.sh'
+    exit 1
+fi
+# -----------------------------------------------------------------------------
+if [ -e build/bug ]
+then
+    rm -r build/bug
+fi
+mkdir -p build/bug
+cd build/bug
+# cmake ../..
+# -----------------------------------------------------------------------------
+cat << EOF
+Issue 43:
+f(x) = (x^0.5) * (x^0.5) = x, and differentiate at x = 0.
+Expect some NaN, Inf, or possibly one, but not zero.
+EOF
+cat << EOF > $name.cpp
+# include <cstdio>
+# include "cppad/cppad.hpp"
+
+
+int main(int argc, char** argv)
+{   bool ok = true;
+
+    using std::cout;
+    using CppAD::AD;
+    using CppAD::vector;
+    //
+    vector< double> x(1), y(1), w(1), dw(1);
+    vector< AD<double> > ax(1), ay(1);
+    //
+    ax[0] = 0.0;
+    //
+    CppAD::Independent(ax);
+    ay[0] = pow(ax[0], 0.5) * pow(ax[0], 0.5);
+    CppAD::ADFun<double> f(ax, ay);
+    //
+    x[0]  = 0.0;
+    y     = f.Forward(0, x);
+    w[0]  = 1.0;
+    dw    = f.Reverse(1, w);
+    //
+    cout << "dw = " << dw << "\n";
+    //
+    ok &= y[0] == 0.0;
+    ok &= dw[0] == 1.0 || ! std::isfinite( dw[0] );
+    //
+    if( ! ok )
+        return 1;
+    return 0;
+}
+EOF
+cxx_flags='-Wall -pedantic-errors -std=c++11 -Wshadow -Wconversion -g -O0'
+eigen_dir="$HOME/prefix/eigen/include"
+echo "g++ -I../../include -isystem $eigen_dir $cxx_flags $name.cpp -o $name"
+g++ -I../../include -isystem $eigen_dir $cxx_flags $name.cpp -o $name
+#
+echo "build/bug/$name"
+if ! ./$name
+then
+    echo
+    echo "build/bug/$name: Error"
+    exit 1
+fi
+echo
+# -----------------------------------------------------------------------------
+echo "bug/$name.sh: OK"
+exit 0

include/cppad/local/pow_op.hpp

@@ -667,13 +667,17 @@ void reverse_powvp_op(
         d, i_z, size_t(arg[0]), cap_order, taylor, nc_partial, partial
     );
 
-    // Special case where partial[ i_z + 2 ] != 0 and x <= 0.
-    // In this case partial[i_z] * (1 / x) is zero, but should not be zero.
-    // (This is a subtitle issue related to the azmul operator.)
+    // For sake of this discussion, consider case where nc_partial = 1.
+    // There is a special case when partial[i_z + 2] != 0 and x <= 0.
+    // In this case partial[i_z] is zero and hence partial[i_z] * (1 / x)
+    // is zero, because reverse_log_op is using the azmul operator for the
+    // multiply. We only want absolute zero multiply when partia[i_z + 2] == 0.
     Base zero(0);
-    if( partial[i_z + 2] != zero  )
+    if( partial[ (i_z + 2) * nc_partial ] != zero  )
     {   if( ! GreaterThanZero( taylor[ size_t(arg[0]) * cap_order ] ) )
-            partial[arg[0]] = numeric_limits<Base>::quiet_NaN();
+            for(size_t k = 0; k < nc_partial; ++k)
+                partial[k + size_t(arg[0]) * nc_partial] =
+                    numeric_limits<Base>::quiet_NaN();
     }
 }
 

test_more/general/pow.cpp

@@ -418,24 +418,28 @@ bool PowTestSeven(void)
     using CppAD::AD;
     using CppAD::vector;
     //
-    vector< double> x(1), y(1), w(1), dw(1);
+    vector< double> x(1), y(1), dx(1), dy(1), w(1), dw(2);
     vector< AD<double> > ax(1), ay(1);
     //
     ax[0] = 0.0;
     //
     CppAD::Independent(ax);
     ay[0] = pow(ax[0], 0.5);
     CppAD::ADFun<double> f(ax, ay);
+    f.check_for_nan(false);
     //
     x[0]  = 0.0;
     y     = f.Forward(0, x);
     //
-    f.check_for_nan(false);
+    dx[0] = 1.0;
+    dy    = f.Forward(1, dx);
+    //
     w[0]  = 1.0;
-    dw    = f.Reverse(1, w);
+    dw    = f.Reverse(2, w);
     //
     ok &= y[0] == 0.0;
     ok &= ! std::isfinite( dw[0] );
+    ok &= ! std::isfinite( dw[1] );
     //
     return ok;
 }