# include <cppad/core/cppad_assert.hpp>
# include <cppad/utility/check_numeric_type.hpp>
# include <cppad/utility/vector.hpp>

// following needed by gcc under fedora 17 so that exp(double) is defined
# include <cppad/base_require.hpp>

namespace CppAD {
    template <class Float, class FloatVector>
    void sparse_hes_fun(
        size_t                       n    ,
        const FloatVector&           x    ,
        const CppAD::vector<size_t>& row  ,
        const CppAD::vector<size_t>& col  ,
        size_t                       p    ,
        FloatVector&                fp    )
    {
        // check numeric type specifications
        CheckNumericType<Float>();

        // check value of p
        CPPAD_ASSERT_KNOWN(
            p == 0 || p == 2,
            "sparse_hes_fun: p != 0 and p != 2"
        );

        size_t K = row.size();
        size_t i, j, k;
        if( p == 0 )
            fp[0] = Float(0);
        else
        {   for(k = 0; k < K; k++)
                fp[k] = Float(0);
        }

        // determine which diagonal entries are present in row[k], col[k]
        CppAD::vector<size_t> diagonal(n);
        for(i = 0; i < n; i++)
            diagonal[i] = K;   // no diagonal entry for this row
        for(k = 0; k < K; k++)
        {   if( row[k] == col[k] )
            {   CPPAD_ASSERT_UNKNOWN( diagonal[row[k]] == K );
                // index of the diagonal entry
                diagonal[ row[k] ] = k;
            }
        }

        // determine which entries must be multiplied by a factor of two
        CppAD::vector<Float> factor(K);
        for(k = 0; k < K; k++)
        {   factor[k] = Float(1);
            for(size_t k1 = 0; k1 < K; k1++)
            {   bool reflected = true;
                reflected &= k != k1;
                reflected &= row[k] != col[k];
                reflected &= row[k] == col[k1];
                reflected &= col[k] == row[k1];
                if( reflected )
                    factor[k] = Float(2);
            }
        }

        Float t;
        for(k = 0; k < K; k++)
        {   i    = row[k];
            j    = col[k];
            t    = exp( x[i] * x[j] );
            switch(p)
            {
                case 0:
                fp[0] += t;
                break;

                case 2:
                if( i == j )
                {   // second partial of t w.r.t. x[i], x[i]
                    fp[k] += ( Float(2) + Float(4) * x[i] * x[i] ) * t;
                }
                else // (i != j)
                {   //
                    // second partial of t w.r.t x[i], x[j]
                    fp[k] += factor[k] * ( Float(1) + x[i] * x[j] ) * t;
                    if( diagonal[i] != K )
                    {   // second partial of t w.r.t x[i], x[i]
                        size_t ki = diagonal[i];
                        fp[ki] += x[j] * x[j] * t;
                    }
                    if( diagonal[j] != K )
                    {   // second partial of t w.r.t x[j], x[j]
                        size_t kj = diagonal[j];
                        fp[kj] += x[i] * x[i] * t;
                    }
                }
                break;
            }
        }

    }
}