Commit d34085a9 authored by Michael Tesch's avatar Michael Tesch
Browse files

fix for #5. restore signgam.

parent f68e7db8

duals/dual

+35 −43
Original line number Diff line number Diff line
@@ -790,11 +790,10 @@ template<class T> dual<T> exp(const dual<T> & x) { 
template<class T> dual<T> log(const dual<T> & x) {

  using std::log;

  T v = log(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / x.rpart());

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / x.rpart());

}



template<class T> dual<T> log10(const dual<T> & x) {
@@ -894,19 +893,19 @@ template<class T> bool (signbit)(const duals::dual<T> & d) { using std::signbit; 
template<class T> dual<T> sqrt(const dual<T> & x) {

  using std::sqrt;

  T v = sqrt(x.rpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / (T(2) * v) );

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / (T(2) * v) );

}



template<class T> dual<T> cbrt(const dual<T> & x) {

  using std::cbrt;

  T v = cbrt(x.rpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / (T(3) * v * v) );

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / (T(3) * v * v) );

}



template<class T> dual<T> sin(const dual<T> & x) {
@@ -933,40 +932,38 @@ template<class T> dual<T> asin(const dual<T> & x) { 
  using std::asin;

  using std::sqrt;

  T v = asin(x.rpart());

  if (x.dpart())

    return dual<T>(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

}



template<class T> dual<T> acos(const dual<T> & x) {

  using std::acos;

  using std::sqrt;

  T v = acos(x.rpart());

  if (x.dpart())

    return dual<T>(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, -x.dpart() / sqrt(1 - x.rpart()*x.rpart()));

}



template<class T> dual<T> atan(const dual<T> & x) {

  using std::atan;

  T v = atan(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / (1 + x.rpart()*x.rpart()));

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / (1 + x.rpart()*x.rpart()));

}



template<class T> dual<T> atan2(const dual<T> & x, const dual<T> & y) {

  using std::atan2;

  T v = atan2(x.rpart(), y.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() / (1 + x.rpart()*x.rpart()));

  else

  if (x.dpart() == T(0))

    return v;

  else

    return dual<T>(v, x.dpart() / (1 + x.rpart()*x.rpart()));

}



// TODO
@@ -986,12 +983,8 @@ template<class T> dual<T> erf(const dual<T> & x) { 
  using std::sqrt;

  using std::pow;

  using std::exp;

  T v = erf(x.rpart());

  if (x.dpart())

    return dual<T>(v,

  return dual<T>(erf(x.rpart()),

                 x.dpart() * T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));

  else

    return v;

}



/// Error function complement (1 - erf()).
@@ -1000,12 +993,8 @@ template<class T> dual<T> erfc(const dual<T> & x) { 
  using std::sqrt;

  using std::pow;

  using std::exp;

  T v = erfc(x.rpart());

  if (x.dpart())

    return dual<T>(v,

  return dual<T>(erfc(x.rpart()),

                 x.dpart() * -T(2)/sqrt(T(MY_PI)) * exp(-pow(x.rpart(),T(2))));

  else

    return v;

}



/// Gamma function.  Approximation of the dual part.
@@ -1013,13 +1002,13 @@ template<class T> dual<T> erfc(const dual<T> & x) { 
template<class T> dual<T> tgamma(const dual<T> & x) {

  using std::tgamma;

  T v = tgamma(x.rpart());

  if (x.dpart()) {

  if (x.dpart() == T(0))

    return v;

  else {

    T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

    return dual<T>(v,

                   x.dpart() * (tgamma(x.rpart()+h) - tgamma(x.rpart()-h))/(2*h));

  }

  else

    return v;

}



/// Log of absolute value of gamma function.  Approximation of the dual part.
@@ -1027,11 +1016,14 @@ template<class T> dual<T> lgamma(const dual<T> & x) { 
  using std::lgamma;

  T h(T(1) / (1ull << (std::numeric_limits<T>::digits / 3)));

  T v = lgamma(x.rpart());

  if (x.dpart())

    return dual<T>(v,

                   x.dpart() * (lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h));

  else

  if (x.dpart() == T(0))

    return v;

  else {

    int signgam_saved = signgam;

    T w = (lgamma(x.rpart()+h) - lgamma(x.rpart()-h))/(2*h);

    signgam = signgam_saved;

    return dual<T>(v, x.dpart() * w);

  }

}



/// Putto operator

tests/sandbox.cpp

+18 −171
Original line number Diff line number Diff line
@@ -84,16 +84,20 @@ int main(int argc, char * argv[]) 
#else



template <class T> T   f(T x) { return pow(x,pow(x,x)); }

template <class T> T df(T x) { return pow(x,-1. + x + pow(x,x)) * (1. + x*log(x) + x*pow(log(x),2.)); }



template <class T> T  df(T x) { return pow(x,-1 + x + pow(x,x)) * (1 + x*log(x) + x*pow(log(x),2)); }

template <class T> T ddf(T x) { return (pow(x,pow(x,x)) * pow(pow(x,x - 1) + pow(x,x)*log(x)*(log(x) + 1), 2) +

                                        pow(x,pow(x,x)) * (pow(x,x - 1) * log(x) +

                                                           pow(x,x - 1) * (log(x) + 1) +

                                                           pow(x,x - 1) * ((x - 1)/x + log(x)) +

                                                           pow(x,x) * log(x) * pow(log(x) + 1, 2) )); }

int main(int argc, char * argv[])

{

  dualf h;

  dualf x(1);

  dualf xx(1);

  hyperdualf y;

  hyperdualf w(1);

  hyperdualf z(x,h);

  hyperdualf a(x);

  hyperdualf z(xx,h);

  hyperdualf a(xx);

  emtx<double> ed;

  emtx<float> ef;

  emtx<complexd> ecd;
@@ -105,171 +109,14 @@ int main(int argc, char * argv[]) 


  std::cout << "  f(2.)            = " << f(2.)    << "\n";

  std::cout << " df(2.)            = " << df(2.)   << "\n";

  std::cout << "ddf(2.)            = " << ddf(2.)  << "\n";

  std::cout << "  f(2+1_e)         = " << f(2+1_e) << "\n";

  std::cout << "  f(2+1_e).dpart() = " << f(2+1_e).dpart() << "\n";



  x = h;

  a = y;

  std::cout << x << "\n" << h << "\n";

  std::cout << w << "\n" << "unitx:" << Eigen::Vector2f::UnitX() << "\n";

  std::cout << z << "\n" << "unity:" << Eigen::Vector2f::UnitY() << "\n";





  std::cout << "M:" << type_name<decltype (duals::rpart(ecf_))>() << "\n";

  std::cout << "M:" << type_name<decltype (duals::rpart(edf_))>() << "\n";

  std::cout << "M:" << type_name<decltype (duals::rpart(ecdf_))>() << "\n";



  std::cout << "M:" << type_name<

    typename Eigen::CwiseUnaryOp<duals::CwiseRpartOp<dualf>,const Eigen::Matrix<duals::dual<float>,2,2>>::Scalar

    >() << "\n";



  std::cout << "N:" << type_name<

    typename Eigen::internal::generic_xpr_base<Eigen::CwiseUnaryOp<duals::CwiseRpartOp<dualf>,

                                                     const Eigen::Matrix<duals::dual<float>,2,2>> >::type

    >() << "\n";



  std::cout << "O:" << type_name<

    typename Eigen::internal::traits<

      Eigen::CwiseUnaryOp<duals::CwiseRpartOp<dualf>, const Eigen::Matrix<duals::dual<float>,2,2>>

      >::Scalar

    >() << "\n";



  std::cout << "o:" << type_name<

    typename Eigen::internal::traits<

      Eigen::CwiseUnaryOp<duals::CwiseRpartOp<dualf>, Eigen::Matrix<duals::dual<float>,2,2>>

      >::Scalar

    >() << "\n";



  std::cout << "a:" << duals::CwiseRpartOp<dualf>()(duals::dual<float>(3,4)) << "\n";



  std::cout << "P:" << type_name<

    typename std::result_of<

      duals::CwiseRpartOp<dualf>(duals::dual<float> &)

      >::type

    >() << "\n";



  std::cout << "Q:" << type_name<

    typename std::result_of<

      duals::CwiseRpartOp<dualf>(const duals::dual<float> &)

      >::type

    >() << "\n";



  std::cout << "S:" << type_name<

    typename std::result_of<

      duals::CwiseRpartOp<dualf>(const typename Eigen::Matrix<duals::dual<float>,2,2>::Scalar &)

      >::type

    >() << "\n";



  return 0;

  duals::hyperduald x(2+1_e,1+0_e);

  std::cout << "  f((2+1_e) + (1+0_e)_e).dpart().dpart() = " << f(x).dpart().dpart() << "\n";

  std::cout << "  c((2+1_e) + (1+0_e)_e).dpart().dpart() = " << cbrt(x).dpart().dpart() << "\n";



  std::cout << type_name<duals::promote<double,complexf>::type>() << "\n";

  std::cout << type_name<duals::promote<dualf,complexf>::type>() << "\n";

  std::cout << type_name<duals::promote<dualf,cdualf>::type>() << "\n";

  std::cout << type_name<duals::promote<dualf,chyperdualf>::type>() << "\n";

  std::cout << type_name<duals::promote<complexf,dualf>::type>() << "\n";

  std::cout << "f-bool:" << type_name<duals::promote<float,std::true_type>::type>() << "\n";

  complexf cf;

  std::cout << type_name<decltype(cf)>() << "\n";

  std::cout << "ecf:" << type_name<decltype(ecf_)>() << "\n";

  std::cout << "ecf:" << type_name<

    typename Eigen::internal::scalar_product_op<complexf,int>>() << "\n";

  std::cout << "dpart(ecdf_)" << type_name<decltype((dpart(ecdf_)).derived())>() << "\n";

  std::cout << "ecf*1:" << type_name<decltype((ecf_ * 1).derived())>() << "\n";

  std::cout << "dpart(ecdf_) po:"

            << type_name<typename std::decay<decltype((rpart(ecdf_)).matrix())>::type::PlainObject>() << "\n";

  std::cout << "ecf*1 po:"

            << type_name<typename std::decay<decltype((ecf_ * 1).matrix())>::type::PlainObject>() << "\n";

  std::cout << "edf:" << type_name<decltype(edf_)>() << "\n";

  std::cout << "edf*1:" << type_name<decltype((edf_ * 1).derived())>() << "\n";

  std::cout << "edf*1 po:"

            << type_name<typename std::decay<decltype((edf_ * 1).matrix())>::type::PlainObject>() << "\n";



  std::cout << type_name<decltype(1_ef)>() << "\n";

  std::cout << type_name<decltype(cf + 1_ef)>() << "\n";

  std::cout << "a+h2 " << type_name<decltype(cf + hyperduald(2))>() << "\n";

  std::cout << "d:"

            << dual_traits<decltype(cf)>::depth << ", "

            << dual_traits<decltype(1_ef)>::depth << "\n";

  std::cout << "---\n";

  auto xx = cf + 1_ef;

  std::cout << "---\n";

  std::cout << xx << "---\n";

  //std::cout << "ecf*edf_:" << type_name<PO_EXPR_TYPE(ecf_ * edf_)>() << "\n";



  std::cout << edf<3>::Constant(3,3, 5 - 8_ef) << "\n";



  std::cout << ecf<3>::Random() + 1 * ecf<3>::Random() << "\n";

  std::cout << ecf<3>::Random() + complexf(1,2) * ecf<3>::Random() << "\n";

  std::cout << edf<3>::Random() + 1_ef * edf<3>::Random() << "\n";

  std::cout << "has ct:" << duals::detail::has_member_type<Randa, int>::value << "\n";

  //std::cout << "has ct:" << duals::detail::has_common_type<decltype(1_ef),

  //                                                         decltype(edf<3>::Random())>::value << "\n";

  std::cout << ecdf<3>::Random(3,3) << "\n";

  auto xyz = ecdf<6,1>::LinSpaced(6,cdualf(3+4_ef,3+4_ef), cdualf(1-5_ef,1-5_ef) );

  std::cout << xyz << "\n-\n";

  std::cout << conj(xyz.array()) << "\n";

  std::cout << "iscomplex? " << Eigen::NumTraits<cdualf>::IsComplex << "\n";

#if 0

  using Eigen::internal::Packet2cf;

  using Eigen::internal::pload;



  std::complex<duals::dual<float>> cd1 = cdualf(1+2_ef,3+4_ef);

  std::complex<duals::dual<float>> cd2 = cdualf(5+6_ef,7+8_ef);

  //std::complex<duals::dual<float>> cd1 = ecdf<1>::Random(1,1)(0) + 1_ef * ecdf<1>::Random(1,1)(0);

  //std::complex<duals::dual<float>> cd2 = ecdf<1>::Random(1,1)(0) + 1_ef * ecdf<1>::Random(1,1)(0);

  std::complex<duals::dual<float>> cd3;

  Packet1cdf p1 = pload<Packet1cdf>(&cd1);

  Packet1cdf p2 = pload<Packet1cdf>(&cd2);

  Packet1cdf p3;

  p3 = pconj(p1);

  pstore(&cd3, p3);

  std::cout << "conj=" << cd3 << "\n";

  p3 = pnegate(p1);

  pstore(&cd3, p3);

  std::cout << "nega=" << cd3 << "\n";

  p3 = pdiv(p1,p2);

  pstore(&cd3, p3);

  std::cout << "div =" << cd3 << "\n";

  std::cout << cd1 << "*\n" << cd2 << " = \n" << cd1 / cd2 << "\n";

  //std::cout << cd1 << "*" << cd2 << " = " << cd3 << "\n";

#endif

#ifdef EIGEN_VECTORIZE_AVX

  using Eigen::internal::Packet4df;

  using namespace Eigen::internal;

  std::vector<cdualf,Eigen::aligned_allocator<cdualf>> da(4), db(4), dc(4);

  da[0] = 1 + 2_ef;

  da[0] = duals::randos::random2<cdualf>();

  da[1] = 3 + 4_ef;

  da[2] = 5 + 6_ef;

  da[3] = 7 + 8_ef;

  db[0] = 1.1 + 2.2_ef;

  db[0] = duals::randos::random2<cdualf>();

  db[1] = 3.3 + 4.4_ef;

  db[2] = 5.5 + 6.6_ef;

  db[3] = 7.7 + 8.8_ef;



  //std::complex<duals::dual<float>> cd1 = ecdf<1>::Random(1,1)(0) + 1_ef * ecdf<1>::Random(1,1)(0);

  //std::complex<duals::dual<float>> cd2 = ecdf<1>::Random(1,1)(0) + 1_ef * ecdf<1>::Random(1,1)(0);



  Packet2cdf p1 = pload<Packet2cdf>(&da[0]);

  Packet2cdf p2 = pload<Packet2cdf>(&db[0]);

  Packet2cdf p3;

  p3 = pmul(p1,p2);

  pstore(&dc[0], p3);

  for (int i = 0; i < 4; i++) {

    std::cout << "conj =" << da[i] << "*" << db[i] << " = " << dc[i] << "\n";

    //std::cout << da[i] << ". " << " = " << conj(da[i]) << " err=" << 0 << "\n";

    //std::cout << da[i] << "." << db[i] << " = " << da[i]/db[i] << " err=" << (dc[i] - da[i]/db[i]) << "\n";

  }

#endif

#if 1

  std::stringstream Xx("(0.24151021242141724+0.31506165862083435_ef)");

  std::stringstream Yy("(3.9031729102134705e-06+0.06223597377538681_ef)");

  dualf xxxx, yyyy;

  Xx >> xxxx;

  Yy >> yyyy;

  std::cout << std::setprecision(10) << "x=" << xxxx << " y=" << yyyy << "\nx*y=" << xxxx / yyyy << "\n";

#endif

}



#endif