C++ std::exp方法代码示例

TEST(AgradFwdLog1pExp,Fvar) {
    using stan::agrad::fvar;
    using stan::math::log1p_exp;
    using std::exp;

    fvar<double> x(0.5,1.0);
    fvar<double> y(1.0,2.0);
    fvar<double> z(2.0,3.0);

    fvar<double> a = log1p_exp(x);
    EXPECT_FLOAT_EQ(log1p_exp(0.5), a.val_);
    EXPECT_FLOAT_EQ(exp(0.5) / (1 + exp(0.5)), a.d_);

    fvar<double> b = log1p_exp(y);
    EXPECT_FLOAT_EQ(log1p_exp(1.0), b.val_);
    EXPECT_FLOAT_EQ(2.0 * exp(1.0) / (1 + exp(1.0)), b.d_);

    fvar<double> a2 = log(1+exp(x));
    EXPECT_FLOAT_EQ(a.d_, a2.d_);

    fvar<double> b2 = log(1+exp(y));
    EXPECT_FLOAT_EQ(b.d_, b2.d_);
}

TEST(AgradFwdLogDiffExp,Double_FvarFvarVar_1stDeriv) {
    using stan::agrad::fvar;
    using stan::agrad::var;
    using stan::math::log_diff_exp;
    using std::exp;

    double x(9.0);
    fvar<fvar<var> > y;
    y.val_.val_ = 6.0;
    y.d_.val_ = 1.0;

    fvar<fvar<var> > a = log_diff_exp(x,y);

    EXPECT_FLOAT_EQ(log_diff_exp(9.0,6.0), a.val_.val_.val());
    EXPECT_FLOAT_EQ(0, a.val_.d_.val());
    EXPECT_FLOAT_EQ(-exp(6.0) / (exp(9.0) - exp(6.0)), a.d_.val_.val());
    EXPECT_FLOAT_EQ(0, a.d_.d_.val());

    AVEC p = createAVEC(y.val_.val_);
    VEC g;
    a.val_.val_.grad(p,g);
    EXPECT_FLOAT_EQ(-exp(6.0) / (exp(9.0) - exp(6.0)), g[0]);
}

TEST(AgradFvar, log1m_inv_logit){
  using stan::agrad::fvar;
  using stan::math::log1m_inv_logit;
  using std::exp;

  fvar<double> x(0.5);
  fvar<double> y(-1.0);
  fvar<double> z(0.0);
  x.d_ = 1.0;
  y.d_ = 2.0;
  z.d_ = 3.0;

  fvar<double> a = log1m_inv_logit(x);
  EXPECT_FLOAT_EQ(log1m_inv_logit(0.5), a.val_);
  EXPECT_FLOAT_EQ(-1.0 * exp(0.5) / (1 + exp(0.5)), a.d_);

  fvar<double> b = log1m_inv_logit(y);
  EXPECT_FLOAT_EQ(log1m_inv_logit(-1.0), b.val_);
  EXPECT_FLOAT_EQ(-2.0 * exp(-1.0) / (1 + exp(-1.0)), b.d_);

  fvar<double> c = log1m_inv_logit(z);
  EXPECT_FLOAT_EQ(log1m_inv_logit(0.0), c.val_);
  EXPECT_FLOAT_EQ(-3.0 * exp(0.0) / (1 + exp(0.0)), c.d_);
}

  void test_exp(const int fuzzy_bits)
  {
    // Use at least 8 resolution bits.
    // Use at least 7 range bits.

    BOOST_STATIC_ASSERT(-FixedPointType::resolution >= 8);
    BOOST_STATIC_ASSERT( FixedPointType::range      >= 7);

    const FixedPointType a1(+1L    );                                       const FloatPointType b1(+1L    );
    const FixedPointType a2(+2L    );                                       const FloatPointType b2(+2L    );
    const FixedPointType a3(+4.375L);                                       const FloatPointType b3(+4.375L);
    const FixedPointType a4(+1.125L);                                       const FloatPointType b4(+1.125L);
    const FixedPointType a5(-1.125L);                                       const FloatPointType b5(-1.125L);
    const FixedPointType a6(+0.875L);                                       const FloatPointType b6(+0.875L);
    const FixedPointType a7(FixedPointType( 1) /  3);                       const FloatPointType b7(FloatPointType( 1) /  3);
    const FixedPointType a8(FixedPointType(11) / 10);                       const FloatPointType b8(FloatPointType(11) / 10);
    const FixedPointType a9(boost::math::constants::phi<FixedPointType>()); const FloatPointType b9(boost::math::constants::phi<FloatPointType>());

    using std::exp;

    BOOST_CHECK_CLOSE_FRACTION(exp(a1), FixedPointType(exp(b1)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a2), FixedPointType(exp(b2)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a3), FixedPointType(exp(b3)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a4), FixedPointType(exp(b4)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a5), FixedPointType(exp(b5)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a6), FixedPointType(exp(b6)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a7), FixedPointType(exp(b7)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a8), FixedPointType(exp(b8)), tolerance_maker<FixedPointType>(fuzzy_bits));
    BOOST_CHECK_CLOSE_FRACTION(exp(a9), FixedPointType(exp(b9)), tolerance_maker<FixedPointType>(fuzzy_bits));
  }

inline fvar<T> expm1(const fvar<T>& x) {
    using std::exp;
    return fvar<T>(expm1(x.val_), x.d_ * exp(x.val_));
}

int test_cv_vib()
{
  using std::exp;

  const Scalar Mm_N  = 14.008e-3;   //in SI kg/mol
  const Scalar Mm_O  = 16e-3;       //in SI kg/mol
  const Scalar Mm_N2 = 2.L * Mm_N;  //in SI kg/mol
  const Scalar Mm_O2 = 2.L * Mm_O;  //in SI kg/mol
  const Scalar Mm_NO = Mm_O + Mm_N; //in SI kg/mol

  std::vector<std::string> species_str_list;
  const unsigned int n_species = 5;
  species_str_list.reserve(n_species);
  species_str_list.push_back( "N2" );
  species_str_list.push_back( "O2" );
  species_str_list.push_back( "N" );
  species_str_list.push_back( "O" );
  species_str_list.push_back( "NO" );

  Antioch::ChemicalMixture<Scalar> chem_mixture( species_str_list );

  // Can we instantiate it?
  Antioch::StatMechThermodynamics<Scalar> sm_thermo( chem_mixture );

  // Mass fractions
  std::vector<Scalar> mass_fractions( 5, 0.2 );
  mass_fractions[0] = 0.5;
  mass_fractions[1] = 0.2;
  mass_fractions[2] = 0.1;
  mass_fractions[3] = 0.1;
  mass_fractions[4] = 0.1;

  const Scalar R_N2 = Antioch::Constants::R_universal<Scalar>() / Mm_N2;
  const Scalar R_O2 = Antioch::Constants::R_universal<Scalar>() / Mm_O2;
  const Scalar R_NO = Antioch::Constants::R_universal<Scalar>() / Mm_NO;

  const Scalar th0_N2 = 3.39500e+03; // degeneracy = 1
  const Scalar th0_O2 = 2.23900e+03; // degeneracy = 1
  const Scalar th0_NO = 2.81700e+03; // degeneracy = 1

  // Tv
  const Scalar Tv = 1000.0;

  const Scalar tol = std::numeric_limits<Scalar>::epsilon() * 2;
  const Scalar ztol = std::numeric_limits<Scalar>::epsilon();

  int return_flag = 0;

  Scalar cv_vib_mix = 0.0;

  // N2
  {
    Scalar cv_vib_N2 = sm_thermo.cv_vib (0, Tv);

    const Scalar expv   = exp(th0_N2/Tv);
    const Scalar expvmi = expv - Scalar(1.0);
    Scalar cv_vib_N2_true = R_N2*th0_N2*th0_N2*expv/expvmi/expvmi/Tv/Tv;

    if( !test_relative(cv_vib_N2, cv_vib_N2_true, tol) )
    {
      std::cerr << std::scientific << std::setprecision(20);
      std::cerr << "Error: Mismatch in cv_vib for N2."
                << "\n Expected = " << cv_vib_N2_true
                << "\n Computed = " << cv_vib_N2
                << "\n Diff     = " << cv_vib_N2_true - cv_vib_N2
                << std::endl;
      return_flag += 1;
    }

    cv_vib_mix += mass_fractions[0]*cv_vib_N2_true;
  }

  // O2
  {
    Scalar cv_vib_O2 = sm_thermo.cv_vib (1, Tv);

    const Scalar expv   = exp(th0_O2/Tv);
    const Scalar expvmi = expv - Scalar(1.0);
    Scalar cv_vib_O2_true = R_O2*th0_O2*th0_O2*expv/expvmi/expvmi/Tv/Tv;

    if( !test_relative(cv_vib_O2, cv_vib_O2_true, tol) )
    {
      std::cerr << std::scientific << std::setprecision(20);
      std::cerr << "Error: Mismatch in cv_vib for O2."
                << "\n Expected = " << cv_vib_O2_true
                << "\n Computed = " << cv_vib_O2
                << "\n Diff     = " << cv_vib_O2_true - cv_vib_O2
                << std::endl;
      return_flag += 1;
    }

    cv_vib_mix += mass_fractions[1]*cv_vib_O2_true;
  }

  // O
  {
    Scalar cv_vib_O = sm_thermo.cv_vib (2, Tv);

    if( !test_zero(cv_vib_O, ztol) )
    {
//.........这里部分代码省略.........

    typename return_type<T_y, T_loc, T_scale, T_shape>::type
    skew_normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma,
                    const T_shape& alpha) {
      static const char* function("skew_normal_log");
      typedef typename stan::partials_return_type<T_y, T_loc,
                                                  T_scale, T_shape>::type
        T_partials_return;

      using std::log;
      using stan::is_constant_struct;
      using std::exp;

      if (!(stan::length(y)
            && stan::length(mu)
            && stan::length(sigma)
            && stan::length(alpha)))
        return 0.0;

      T_partials_return logp(0.0);

      check_not_nan(function, "Random variable", y);
      check_finite(function, "Location parameter", mu);
      check_finite(function, "Shape parameter", alpha);
      check_positive(function, "Scale parameter", sigma);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Location parameter", mu,
                             "Scale parameter", sigma,
                             "Shape paramter", alpha);

      if (!include_summand<propto, T_y, T_loc, T_scale, T_shape>::value)
        return 0.0;

      OperandsAndPartials<T_y, T_loc, T_scale, T_shape>
        operands_and_partials(y, mu, sigma, alpha);

      using std::log;

      VectorView<const T_y> y_vec(y);
      VectorView<const T_loc> mu_vec(mu);
      VectorView<const T_scale> sigma_vec(sigma);
      VectorView<const T_shape> alpha_vec(alpha);
      size_t N = max_size(y, mu, sigma, alpha);

      VectorBuilder<true, T_partials_return, T_scale> inv_sigma(length(sigma));
      VectorBuilder<include_summand<propto, T_scale>::value,
                    T_partials_return, T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto, T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

      for (size_t n = 0; n < N; n++) {
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return mu_dbl = value_of(mu_vec[n]);
        const T_partials_return sigma_dbl = value_of(sigma_vec[n]);
        const T_partials_return alpha_dbl = value_of(alpha_vec[n]);

        const T_partials_return y_minus_mu_over_sigma
          = (y_dbl - mu_dbl) * inv_sigma[n];
        const double pi_dbl = pi();

        if (include_summand<propto>::value)
          logp -=  0.5 * log(2.0 * pi_dbl);
        if (include_summand<propto, T_scale>::value)
          logp -= log(sigma_dbl);
        if (include_summand<propto, T_y, T_loc, T_scale>::value)
          logp -= y_minus_mu_over_sigma * y_minus_mu_over_sigma / 2.0;
        if (include_summand<propto, T_y, T_loc, T_scale, T_shape>::value)
          logp += log(erfc(-alpha_dbl * y_minus_mu_over_sigma
                           / std::sqrt(2.0)));

        T_partials_return deriv_logerf
          = 2.0 / std::sqrt(pi_dbl)
          * exp(-alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0)
                * alpha_dbl * y_minus_mu_over_sigma / std::sqrt(2.0))
          / (1 + erf(alpha_dbl * y_minus_mu_over_sigma
                     / std::sqrt(2.0)));
        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n]
            += -y_minus_mu_over_sigma / sigma_dbl
            + deriv_logerf * alpha_dbl / (sigma_dbl * std::sqrt(2.0));
        if (!is_constant_struct<T_loc>::value)
          operands_and_partials.d_x2[n]
            += y_minus_mu_over_sigma / sigma_dbl
            + deriv_logerf * -alpha_dbl / (sigma_dbl * std::sqrt(2.0));
        if (!is_constant_struct<T_scale>::value)
          operands_and_partials.d_x3[n]
            += -1.0 / sigma_dbl
            + y_minus_mu_over_sigma * y_minus_mu_over_sigma / sigma_dbl
            - deriv_logerf * y_minus_mu_over_sigma * alpha_dbl
            / (sigma_dbl * std::sqrt(2.0));
        if (!is_constant_struct<T_shape>::value)
          operands_and_partials.d_x4[n]
            += deriv_logerf * y_minus_mu_over_sigma / std::sqrt(2.0);
      }
      return operands_and_partials.value(logp);
    }

    typename return_type<T_y, T_scale_succ, T_scale_fail>::type
    beta_cdf_log(const T_y& y, const T_scale_succ& alpha,
                 const T_scale_fail& beta) {
      typedef typename stan::partials_return_type<T_y, T_scale_succ,
                                                  T_scale_fail>::type
        T_partials_return;

      // Size checks
      if ( !( stan::length(y) && stan::length(alpha)
              && stan::length(beta) ) )
        return 0.0;

      // Error checks
      static const char* function("stan::math::beta_cdf");

      using stan::math::check_positive_finite;
      using stan::math::check_not_nan;
      using stan::math::check_nonnegative;
      using stan::math::check_less_or_equal;
      using boost::math::tools::promote_args;
      using stan::math::check_consistent_sizes;
      using stan::math::value_of;

      T_partials_return cdf_log(0.0);

      check_positive_finite(function, "First shape parameter", alpha);
      check_positive_finite(function, "Second shape parameter", beta);
      check_not_nan(function, "Random variable", y);
      check_nonnegative(function, "Random variable", y);
      check_less_or_equal(function, "Random variable", y, 1);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "First shape parameter", alpha,
                             "Second shape parameter", beta);

      // Wrap arguments in vectors
      VectorView<const T_y> y_vec(y);
      VectorView<const T_scale_succ> alpha_vec(alpha);
      VectorView<const T_scale_fail> beta_vec(beta);
      size_t N = max_size(y, alpha, beta);

      OperandsAndPartials<T_y, T_scale_succ, T_scale_fail>
        operands_and_partials(y, alpha, beta);

      // Compute CDF and its gradients
      using stan::math::inc_beta;
      using stan::math::digamma;
      using stan::math::lbeta;
      using std::pow;
      using std::exp;
      using std::log;
      using std::exp;

      // Cache a few expensive function calls if alpha or beta is a parameter
      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_alpha_vec(max_size(alpha, beta));

      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_beta_vec(max_size(alpha, beta));

      VectorBuilder<contains_nonconstant_struct<T_scale_succ,
                                                T_scale_fail>::value,
                    T_partials_return, T_scale_succ, T_scale_fail>
        digamma_sum_vec(max_size(alpha, beta));

      if (contains_nonconstant_struct<T_scale_succ, T_scale_fail>::value) {
        for (size_t i = 0; i < N; i++) {
          const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
          const T_partials_return beta_dbl = value_of(beta_vec[i]);

          digamma_alpha_vec[i] = digamma(alpha_dbl);
          digamma_beta_vec[i] = digamma(beta_dbl);
          digamma_sum_vec[i] = digamma(alpha_dbl + beta_dbl);
        }
      }

      // Compute vectorized CDFLog and gradient
      for (size_t n = 0; n < N; n++) {
        // Pull out values
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
        const T_partials_return beta_dbl = value_of(beta_vec[n]);
        const T_partials_return betafunc_dbl = exp(lbeta(alpha_dbl, beta_dbl));
        // Compute
        const T_partials_return Pn = inc_beta(alpha_dbl, beta_dbl, y_dbl);

        cdf_log += log(Pn);

        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] += pow(1-y_dbl, beta_dbl-1)
            * pow(y_dbl, alpha_dbl-1) / betafunc_dbl / Pn;

        T_partials_return g1 = 0;
        T_partials_return g2 = 0;

        if (contains_nonconstant_struct<T_scale_succ, T_scale_fail>::value) {
//.........这里部分代码省略.........

DualNumber<T> exp(DualNumber<T> x)
{
    using std::exp;
    return DualNumber<T>(exp(x.first), x.second * exp(x.first));
}

    typename return_type<T_prob>::type
    binomial_cdf(const T_n& n, const T_N& N, const T_prob& theta) {
      static const char* function("stan::prob::binomial_cdf");
      typedef typename stan::partials_return_type<T_n,T_N,T_prob>::type
        T_partials_return;
          
      using stan::math::check_finite;
      using stan::math::check_bounded;
      using stan::math::check_nonnegative;
      using stan::math::value_of;
      using stan::math::check_consistent_sizes;
      using stan::prob::include_summand;
          
      // Ensure non-zero arguments lenghts
      if (!(stan::length(n) && stan::length(N) && stan::length(theta)))
        return 1.0;
          
      T_partials_return P(1.0);
          
      // Validate arguments
      check_nonnegative(function, "Population size parameter", N);
      check_finite(function, "Probability parameter", theta);
      check_bounded(function, "Probability parameter", theta, 0.0, 1.0);
      check_consistent_sizes(function, 
                             "Successes variable", n,
                             "Population size parameter", N, 
                             "Probability parameter", theta);
                             
          
      // Wrap arguments in vector views
      VectorView<const T_n> n_vec(n);
      VectorView<const T_N> N_vec(N);
      VectorView<const T_prob> theta_vec(theta);
      size_t size = max_size(n, N, theta);
          
      // Compute vectorized CDF and gradient
      using stan::math::value_of;
      using stan::math::inc_beta;
      using stan::math::lbeta;
      using std::exp;
      using std::pow;
          
      agrad::OperandsAndPartials<T_prob> operands_and_partials(theta);
          
      // Explicit return for extreme values
      // The gradients are technically ill-defined, but treated as zero
      for (size_t i = 0; i < stan::length(n); i++) {
        if (value_of(n_vec[i]) < 0) 
          return operands_and_partials.to_var(0.0,theta);
      }
        
      for (size_t i = 0; i < size; i++) {
              
        // Explicit results for extreme values
        // The gradients are technically ill-defined, but treated as zero
        if (value_of(n_vec[i]) >= value_of(N_vec[i])) {
          continue;
        }
          
        const T_partials_return n_dbl = value_of(n_vec[i]);
        const T_partials_return N_dbl = value_of(N_vec[i]);
        const T_partials_return theta_dbl = value_of(theta_vec[i]);
        const T_partials_return betafunc = exp(lbeta(N_dbl-n_dbl,n_dbl+1));
        const T_partials_return Pi = inc_beta(N_dbl - n_dbl, n_dbl + 1, 
                                              1 - theta_dbl);
          
        P *= Pi;

        if (!is_constant_struct<T_prob>::value)
          operands_and_partials.d_x1[i] -= pow(theta_dbl,n_dbl)
            * pow(1-theta_dbl,N_dbl-n_dbl-1) / betafunc / Pi;
      }
          
      if (!is_constant_struct<T_prob>::value) {
        for(size_t i = 0; i < stan::length(theta); ++i)
          operands_and_partials.d_x1[i] *= P;
      }
          
      return operands_and_partials.to_var(P,theta);
        
    }

    typename return_type<T_y, T_dof>::type
    inv_chi_square_cdf_log(const T_y& y, const T_dof& nu) {
      typedef typename stan::partials_return_type<T_y, T_dof>::type
        T_partials_return;

      // Size checks
      if ( !( stan::length(y) && stan::length(nu) ) ) return 0.0;

      // Error checks
      static const char* function("stan::math::inv_chi_square_cdf_log");

      using stan::math::check_positive_finite;
      using stan::math::check_not_nan;
      using stan::math::check_consistent_sizes;
      using stan::math::check_nonnegative;
      using boost::math::tools::promote_args;
      using stan::math::value_of;
      using std::exp;

      T_partials_return P(0.0);

      check_positive_finite(function, "Degrees of freedom parameter", nu);
      check_not_nan(function, "Random variable", y);
      check_nonnegative(function, "Random variable", y);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Degrees of freedom parameter", nu);

      // Wrap arguments in vectors
      VectorView<const T_y> y_vec(y);
      VectorView<const T_dof> nu_vec(nu);
      size_t N = max_size(y, nu);

      OperandsAndPartials<T_y, T_dof> operands_and_partials(y, nu);

      // Explicit return for extreme values
      // The gradients are technically ill-defined, but treated as zero

      for (size_t i = 0; i < stan::length(y); i++)
        if (value_of(y_vec[i]) == 0)
          return operands_and_partials.to_var(stan::math::negative_infinity(),
                                              y, nu);

      // Compute cdf_log and its gradients
      using stan::math::gamma_q;
      using stan::math::digamma;
      using boost::math::tgamma;
      using std::exp;
      using std::pow;
      using std::log;

      // Cache a few expensive function calls if nu is a parameter
      VectorBuilder<!is_constant_struct<T_dof>::value,
                    T_partials_return, T_dof> gamma_vec(stan::length(nu));
      VectorBuilder<!is_constant_struct<T_dof>::value,
                    T_partials_return, T_dof> digamma_vec(stan::length(nu));

      if (!is_constant_struct<T_dof>::value)  {
        for (size_t i = 0; i < stan::length(nu); i++) {
          const T_partials_return nu_dbl = value_of(nu_vec[i]);
          gamma_vec[i] = tgamma(0.5 * nu_dbl);
          digamma_vec[i] = digamma(0.5 * nu_dbl);
        }
      }

      // Compute vectorized cdf_log and gradient
      for (size_t n = 0; n < N; n++) {
        // Explicit results for extreme values
        // The gradients are technically ill-defined, but treated as zero
        if (value_of(y_vec[n]) == std::numeric_limits<double>::infinity()) {
          continue;
        }

        // Pull out values
        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return y_inv_dbl = 1.0 / y_dbl;
        const T_partials_return nu_dbl = value_of(nu_vec[n]);

        // Compute
        const T_partials_return Pn = gamma_q(0.5 * nu_dbl, 0.5 * y_inv_dbl);

        P += log(Pn);

        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] += 0.5 * y_inv_dbl * y_inv_dbl
            * exp(-0.5*y_inv_dbl) * pow(0.5*y_inv_dbl, 0.5*nu_dbl-1)
            / tgamma(0.5*nu_dbl) / Pn;
        if (!is_constant_struct<T_dof>::value)
          operands_and_partials.d_x2[n]
            += 0.5 * stan::math::grad_reg_inc_gamma(0.5 * nu_dbl,
                                                    0.5 * y_inv_dbl,
                                                    gamma_vec[n],
                                                    digamma_vec[n]) / Pn;
      }

      return operands_and_partials.to_var(P, y, nu);
    }

    typename return_type<T_y, T_dof, T_scale>::type
    scaled_inv_chi_square_ccdf_log(const T_y& y, const T_dof& nu,
                                   const T_scale& s) {
      typedef typename stan::partials_return_type<T_y, T_dof, T_scale>::type
        T_partials_return;

      if (!(stan::length(y) && stan::length(nu) && stan::length(s)))
        return 0.0;

      static const char* function("scaled_inv_chi_square_ccdf_log");

      using std::exp;

      T_partials_return P(0.0);

      check_not_nan(function, "Random variable", y);
      check_nonnegative(function, "Random variable", y);
      check_positive_finite(function, "Degrees of freedom parameter", nu);
      check_positive_finite(function, "Scale parameter", s);
      check_consistent_sizes(function,
                             "Random variable", y,
                             "Degrees of freedom parameter", nu,
                             "Scale parameter", s);

      VectorView<const T_y> y_vec(y);
      VectorView<const T_dof> nu_vec(nu);
      VectorView<const T_scale> s_vec(s);
      size_t N = max_size(y, nu, s);

      OperandsAndPartials<T_y, T_dof, T_scale>
        operands_and_partials(y, nu, s);

      // Explicit return for extreme values
      // The gradients are technically ill-defined, but treated as zero
      for (size_t i = 0; i < stan::length(y); i++) {
        if (value_of(y_vec[i]) == 0)
          return operands_and_partials.value(0.0);
      }

      using std::exp;
      using std::pow;
      using std::log;

      VectorBuilder<!is_constant_struct<T_dof>::value,
                    T_partials_return, T_dof> gamma_vec(stan::length(nu));
      VectorBuilder<!is_constant_struct<T_dof>::value,
                    T_partials_return, T_dof> digamma_vec(stan::length(nu));

      if (!is_constant_struct<T_dof>::value) {
        for (size_t i = 0; i < stan::length(nu); i++) {
          const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
          gamma_vec[i] = tgamma(half_nu_dbl);
          digamma_vec[i] = digamma(half_nu_dbl);
        }
      }

      for (size_t n = 0; n < N; n++) {
        // Explicit results for extreme values
        // The gradients are technically ill-defined, but treated as zero
        if (value_of(y_vec[n]) == std::numeric_limits<double>::infinity()) {
          return operands_and_partials.value(negative_infinity());
        }

        const T_partials_return y_dbl = value_of(y_vec[n]);
        const T_partials_return y_inv_dbl = 1.0 / y_dbl;
        const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
        const T_partials_return s_dbl = value_of(s_vec[n]);
        const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl
          * y_inv_dbl;
        const T_partials_return half_nu_s2_overx_dbl
          = 2.0 * half_nu_dbl * half_s2_overx_dbl;

        const T_partials_return Pn = gamma_p(half_nu_dbl,
                                             half_nu_s2_overx_dbl);
        const T_partials_return gamma_p_deriv = exp(-half_nu_s2_overx_dbl)
          * pow(half_nu_s2_overx_dbl, half_nu_dbl-1) / tgamma(half_nu_dbl);

        P += log(Pn);

        if (!is_constant_struct<T_y>::value)
          operands_and_partials.d_x1[n] -= half_nu_s2_overx_dbl * y_inv_dbl
            * gamma_p_deriv / Pn;
        if (!is_constant_struct<T_dof>::value)
          operands_and_partials.d_x2[n]
            -= (0.5 * grad_reg_inc_gamma(half_nu_dbl,
                                         half_nu_s2_overx_dbl,
                                         gamma_vec[n],
                                         digamma_vec[n])
                - half_s2_overx_dbl * gamma_p_deriv)
            / Pn;
        if (!is_constant_struct<T_scale>::value)
          operands_and_partials.d_x3[n] += 2.0 * half_nu_dbl * s_dbl * y_inv_dbl
            * gamma_p_deriv / Pn;
      }
      return operands_and_partials.value(P);
    }

 inline T0
 operator()(const T0& arg1) const {
   return exp(arg1);
 }

TEST(AgradFwdExp,Fvar) {
  using stan::math::fvar;
  using std::exp;

  fvar<double> x(0.5,1.0);
  
  fvar<double> a = exp(x);
  EXPECT_FLOAT_EQ(exp(0.5), a.val_);
  EXPECT_FLOAT_EQ(exp(0.5), a.d_);

  fvar<double> b = 2 * exp(x) + 4;
  EXPECT_FLOAT_EQ(2 * exp(0.5) + 4, b.val_);
  EXPECT_FLOAT_EQ(2 * exp(0.5), b.d_);

  fvar<double> c = -exp(x) + 5;
  EXPECT_FLOAT_EQ(-exp(0.5) + 5, c.val_);
  EXPECT_FLOAT_EQ(-exp(0.5), c.d_);

  fvar<double> d = -3 * exp(-x) + 5 * x;
  EXPECT_FLOAT_EQ(-3 * exp(-0.5) + 5 * 0.5, d.val_);
  EXPECT_FLOAT_EQ(3 * exp(-0.5) + 5, d.d_);

  fvar<double> y(-0.5,1.0);
  fvar<double> e = exp(y);
  EXPECT_FLOAT_EQ(exp(-0.5), e.val_);
  EXPECT_FLOAT_EQ(exp(-0.5), e.d_);

  fvar<double> z(0.0,1.0);
  fvar<double> f = exp(z);
  EXPECT_FLOAT_EQ(exp(0.0), f.val_);
  EXPECT_FLOAT_EQ(exp(0.0), f.d_);
}

double exact_soln(double x){
    return 1.0 - (1.0-exp(-10.0))*x - exp(-10.0*x);
}