Commit 827442af authored by Michael Tesch's avatar Michael Tesch
Browse files

Merge branch 'ft-fix-questionable' into 'master' 

convenient but questionable maths

See merge request !34
parents 317bb1ff 546b51fe

.gitlab-ci.yml

+1 −0
Original line number Diff line number Diff line
@@ -91,6 +91,7 @@ upload: 
    - if: $CI_COMMIT_TAG

  script:

    - |

      tar czvf cppduals-h-${CI_BUILD_TAG#v}.tgz duals/

      curl --header "JOB-TOKEN: ${CI_JOB_TOKEN}" \

        --upload-file cppduals-h-${CI_BUILD_TAG#v}.tgz \

                      ${PACKAGE_REGISTRY_URL}/cppduals-h-${CI_BUILD_TAG#v}.tgz

CMakeLists.txt

+1 −1
Original line number Diff line number Diff line
@@ -3,7 +3,7 @@ 
#

cmake_minimum_required (VERSION 3.12)

project (cppduals

  VERSION 0.5.3

  VERSION 0.5.4

  LANGUAGES C CXX

  )

include (GNUInstallDirs)

duals/dual

+92 −48
Original line number Diff line number Diff line
@@ -728,6 +728,45 @@ CPPDUALS_LHS_COMPARISON(!=) 


#endif // PARSED_BY_DOXYGEN



} // namespace duals



// some useful functions of solely the real part

#ifdef _LIBCPP_BEGIN_NAMESPACE_STD

_LIBCPP_BEGIN_NAMESPACE_STD

#else

namespace std {

#endif



#ifndef PARSED_BY_DOXYGEN



#define make_math(T)                                                    \

  inline T (frexp)(const duals::dual<T> & arg, int* exp ) { return (frexp)(arg.rpart(), exp); } \

  inline duals::dual<T> (ldexp)(const duals::dual<T> & arg, int exp ) { return arg * std::pow((T)2,exp); } \

  inline T (trunc)(const duals::dual<T> & d) { return (trunc)(d.rpart()); } \

  inline T (floor)(const duals::dual<T> & d) { return (floor)(d.rpart()); } \

  inline T (ceil)(const duals::dual<T> & d)  { return (ceil)(d.rpart()); } \

  inline T (round)(const duals::dual<T> & d) { return (round)(d.rpart()); } \

  inline int (fpclassify)(const duals::dual<T> & d) { return (fpclassify)(d.rpart()); } \

  inline bool (isfinite)(const duals::dual<T> & d) { return (isfinite)(d.rpart()); } \

  inline bool (isnormal)(const duals::dual<T> & d) { return (isnormal)(d.rpart()); } \

  inline bool (isinf)(const duals::dual<T> & d) { return (isinf)(d.rpart()); } \

  inline bool (isnan)(const duals::dual<T> & d) { return (isnan)(d.rpart()); } \



make_math(float)

make_math(double)

make_math(long double)



#undef make_math



#endif // PARSED_BY_DOXYGEN



#ifdef _LIBCPP_BEGIN_NAMESPACE_STD

_LIBCPP_END_NAMESPACE_STD

#else

} // namespace std

#endif

namespace duals {



namespace randos {



// Random real value between a and b.
@@ -833,24 +872,15 @@ template<class T> dual<T> log2(const dual<T> & x) { 
  return log(x) / log(static_cast<T>(2));

}



template<class T, class U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)>

common_t

pow(const dual<T> & f, const dual<U> & g) {

  using std::pow;

  using std::log;

  T v = pow(f.rpart(), g.rpart());

  return common_t(v,

                  pow(f.rpart(), g.rpart() - T(1)) *

                  (g.rpart() * f.dpart()

                   + f.rpart() * log(f.rpart()) * g.dpart()));

}



template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>

common_t

pow(const dual<T> & x, const U & y) {

  using std::pow;

  typedef typename common_t::value_type V;



  return common_t(pow(x.rpart(), y),

                  x.dpart() * y * pow(x.rpart(), y - U(1)));

                  x.dpart() == V(0) ? V(0)

                  : x.dpart() * y * pow(x.rpart(), y - U(1)));

}



template<class T, class U, CPPDUALS_ENABLE_LEQ_DEPTH_AND_COMMON_T(T,U)>
@@ -859,6 +889,55 @@ pow(const U & x, const dual<T> & y) { 
  return pow(common_t(x), y);

}



template<class T, class U, CPPDUALS_ENABLE_SAME_DEPTH_AND_COMMON_T(T,U)>

common_t

pow(const dual<T> & f, const dual<U> & g) {

  using std::pow;

  using std::log;

#if 1

  using std::floor;

  typedef typename common_t::value_type V;

  common_t result;



 if (f.rpart() == T(0) && g.rpart() >= U(1)) {

    if (g.rpart() > U(1)) {

      result = common_t(0);

    } else {

      result = f;

    }



  } else {

    if (f.rpart() < T(0) && g.rpart() == floor(g.rpart())) {

      V const tmp = g.rpart() * pow(f.rpart(), g.rpart() - U(1.0));

      result = common_t(pow(f.rpart(), g.rpart()),

                        f.dpart() == T(0) ? T(0) : f.dpart() * tmp);

      if (g.dpart() != U(0.0)) {

        // Return a NaN when g.dpart() != 0.

        result.dpart(std::numeric_limits<T>::quiet_NaN());

      }

    } else {

      // Handle the remaining cases. For cases 4,5,6,9 we allow the log()

      // function to generate -HUGE_VAL or NaN, since those cases result in a

      // nonfinite derivative.

      V const tmp1 = pow(f.rpart(), g.rpart());

      V const tmp2 = g.rpart() * pow(f.rpart(), g.rpart() - T(1.0));

      V const tmp3 = tmp1 * log(f.rpart());

      result = common_t(tmp1,

                        f.dpart() == T(0) && g.dpart() == U(0) ? V(0)

                        : tmp2 * f.dpart() + tmp3 * g.dpart());

    }

  }



  return result;

#else

  T v = pow(f.rpart(), g.rpart());

  return common_t(v,

                  pow(f.rpart(), g.rpart() - T(1)) *

                  (g.rpart() * f.dpart()

                   + f.rpart() * log(f.rpart()) * g.dpart()));

#endif

}



namespace utils {

template <typename T> int sgn(T val) {

  return (T(0) < val) - (val < T(0));
@@ -1372,39 +1451,4 @@ struct fmt::formatter<std::complex<T>,Char> : public fmt::formatter<T,Char> 
};

#endif



#ifdef _LIBCPP_BEGIN_NAMESPACE_STD

_LIBCPP_BEGIN_NAMESPACE_STD

#else

namespace std {

#endif



#ifndef PARSED_BY_DOXYGEN



#define make_math(T)                                                    \

  inline T (frexp)(const duals::dual<T> & arg, int* exp ) { return (frexp)(arg.rpart(), exp); } \

  inline duals::dual<T> (ldexp)(const duals::dual<T> & arg, int exp ) { return arg * std::pow((T)2,exp); } \

  inline T (trunc)(const duals::dual<T> & d) { return (trunc)(d.rpart()); } \

  inline T (floor)(const duals::dual<T> & d) { return (floor)(d.rpart()); } \

  inline T (ceil)(const duals::dual<T> & d)  { return (ceil)(d.rpart()); } \

  inline T (round)(const duals::dual<T> & d) { return (round)(d.rpart()); } \

  inline int (fpclassify)(const duals::dual<T> & d) { return (fpclassify)(d.rpart()); } \

  inline bool (isfinite)(const duals::dual<T> & d) { return (isfinite)(d.rpart()); } \

  inline bool (isnormal)(const duals::dual<T> & d) { return (isnormal)(d.rpart()); } \

  inline bool (isinf)(const duals::dual<T> & d) { return (isinf)(d.rpart()); } \

  inline bool (isnan)(const duals::dual<T> & d) { return (isnan)(d.rpart()); } \



make_math(float)

make_math(double)

make_math(long double)



#undef make_math



#endif // PARSED_BY_DOXYGEN



#ifdef _LIBCPP_BEGIN_NAMESPACE_STD

_LIBCPP_END_NAMESPACE_STD

#else

} // namespace std

#endif



#endif // CPPDUALS_DUAL

tests/test_dual.cpp

+12 −0
Original line number Diff line number Diff line
@@ -186,6 +186,8 @@ TEST(template_, common_type) { 
  _EXPECT_TRUE(std::is_same<decltype(cd*cd), cdualf>);

  _EXPECT_TRUE(std::is_same<decltype(cd*d), cdualf>);

  _EXPECT_TRUE(std::is_same<decltype(d*cd), cdualf>);

  _EXPECT_TRUE(std::is_same<decltype(cd*a), cdualf>);

  _EXPECT_TRUE(std::is_same<decltype(a*cd), cdualf>);



  _EXPECT_TRUE(std::is_same<decltype(cd*1), cdualf>);

  _EXPECT_TRUE(std::is_same<decltype(1*cd), cdualf>);
@@ -328,6 +330,7 @@ TEST(members, rpart) { 
  EXPECT_EQ(z.rpart(), 4);

  EXPECT_EQ(z.dpart(), 3);

}



TEST(members, dpart) {

  EXPECT_EQ(dpart(3), 0);

  dualf z(2,3);
@@ -555,6 +558,15 @@ TEST(comparison, ge) { 
  EXPECT_FALSE(1 >= a);

}



#if 0

TEST(simple_ops, add) {

  // https://gitlab.com/tesch1/cppduals/-/issues/11

  duals::dual<std::complex<double>> x;

  x = x+std::complex<double>(1., 2.);

  duals::dual<std::complex<double>> y = sqrt(x+std::complex<double>(1., 2.));

}

#endif



TEST(compound_assign, same_type) {

  // OP=

  dualf x = 2 + 4_e;

tests/test_funcs.cpp

+13 −0
Original line number Diff line number Diff line
@@ -129,6 +129,19 @@ TEST(func, abs) { 
TEST(func, arg) {

  // TODO

}

TEST(func, pow) {

  dualf x = pow(0 + 0_e, 0.);

  EXPECT_EQ(rpart(x), 1);

  EXPECT_EQ(dpart(x), 0);



  dualf y = pow(0 + 0_e, 0. + 0_e);

  EXPECT_EQ(rpart(y), 1);

  EXPECT_EQ(dpart(y), 0);



  dualf z = pow(0, 0. + 0_e);

  EXPECT_EQ(rpart(z), 1);

  EXPECT_EQ(dpart(z), 0);

}

TEST(func, norm) {

  // TODO

}