00001 #ifndef MTH_AD_H
00002 #define MTH_AD_H
00003
00004 #include <cmath>
00005 #include <iostream>
00006 #include "canArray.h"
00007
00011 namespace mth {
00012
00014 template <class T=double, unsigned int N=0>
00015 class AD
00016 {
00017 public:
00019 double x_;
00021 T dx_[N];
00023 enum { degree = N };
00025 AD():x_(0.) {zero();}
00027 AD(double x):x_(x) {zero();}
00029 AD(AD<T, N> const& other) {copy(other);}
00030 template <class B>
00031 AD(AD<B, N> const& other) {copy(other);}
00033 unsigned int size() const {return N;}
00035 void diff(unsigned int i, unsigned int n=0)
00036 {
00037 (void)n;
00038 zero();
00039 dx_[i] = 1.;
00040 }
00042 double& val() {return x_;}
00044 const double& val() const {return x_;}
00046 T& dx(unsigned int i) {return dx_[i];}
00048 const T& dx(unsigned int i) const {return dx_[i];}
00050 inline void resize(unsigned int i) {(void)i;}
00052 operator double() const {return double(x_);}
00054 AD<T, N>& operator=(double other)
00055 {
00056 x_ = other;
00057 zero();
00058 return *this;
00059 }
00061 AD<T, N>& operator=(AD<T, N> const& other)
00062 {
00063 copy(other);
00064 return *this;
00065 }
00066
00067 template <class B>
00068 AD<T, N>& operator=(AD<B, N> const& other)
00069 {
00070 copy(other);
00071 return *this;
00072 }
00074 AD<T, N>& operator+=(double other)
00075 {
00076 x_ += other;
00077 return *this;
00078 }
00080 template< class B>
00081 AD<T, N>& operator+=(AD<B, N> const& other)
00082 {
00083 x_ += other.x_;
00084 for (unsigned int i=0; i < N; ++i)
00085 dx_[i] += other.dx_[i];
00086 return *this;
00087 }
00089 AD<T, N>& operator-=(double other)
00090 {
00091 x_ -= other;
00092 return *this;
00093 }
00095 template< class B>
00096 AD<T, N>& operator-=(AD<B, N> const& other)
00097 {
00098 x_ -= other.x_;
00099 for (unsigned int i=0; i < N; ++i)
00100 dx_[i] -= other.dx_[i];
00101 return *this;
00102 }
00104 AD<T, N>& operator*=(double other)
00105 {
00106 x_ *= other;
00107 for (unsigned int i=0; i < N; ++i)
00108 dx_[i] *= other;
00109 return *this;
00110 }
00112 template<class B>
00113 AD<T, N>& operator*=(AD<B, N> const& other)
00114 {
00115 x_ *= other.x_;
00116 for (unsigned int i=0; i < N; ++i)
00117 dx_[i] = dx_[i]*other + x_*other.dx_[i];
00118 }
00120 AD<T, N>& operator/=(double other)
00121 {
00122 x_ /= other;
00123 for (unsigned int i=0; i < N; ++i)
00124 dx_[i] /= other;
00125 return *this;
00126 }
00128 template<class B>
00129 AD<T, N>& operator/=(AD<B, N> const& other)
00130 {
00131 x_ /= other.x_;
00132 for (unsigned int i=0; i < N; ++i)
00133 dx_[i] = (dx_[i]*other - other*other.dx_[i]) / (other*other);
00134 }
00135 private:
00136 void zero()
00137 {
00138 for (unsigned int i=0; i < N; ++i)
00139 dx_[i] = 0.;
00140 }
00141
00142 template <class B>
00143 void copy(AD<B, N> const& other)
00144 {
00145 this->x_ = other.x_;
00146 for (unsigned int i=0; i < N; ++i)
00147 this->dx_[i] = (T)other.dx_[i];
00148 }
00149 };
00150
00152 template <class T>
00153 class AD<T,0>
00154 {
00155 public:
00157 double x_;
00159 can::Array<T> dx_;
00160 const static double _zero_;
00162 AD():x_(0), dx_() {zero();}
00164 AD(double val):x_(val), dx_() {zero();}
00166 template<class B>
00167 AD(AD<B, 0> const& other) {zero(); copy(other);}
00168 unsigned size() { return dx_.size();}
00169 unsigned size() const { return dx_.size();}
00170 void diff(unsigned int i, unsigned int n)
00171 {
00172 dx_.resize_copy(n);
00173 zero();
00174 dx_[i] = 1;
00175 }
00177 double& val() {return x_;}
00179 const double& val() const {return x_;}
00181 T& dx(unsigned int i)
00182 {
00183 if(i >= size())
00184 {
00185 unsigned int temp_size = size();
00186 dx_.resize_copy(i + 1);
00187 for(unsigned int x = temp_size; x <= i; x++)
00188 dx_[x] = 0.;
00189 }
00190 return dx_[i];
00191 }
00193 const T dx(unsigned int i) const {
00194 if(i >= size())
00195 return T(_zero_);
00196 return dx_[i];
00197 }
00199 operator double() const {return double(x_);}
00201 AD<T, 0>& operator=(double other)
00202 {
00203 x_ = other;
00204 zero();
00205 return *this;
00206 }
00208 template<class B>
00209 AD<T, 0>& operator=(AD<B, 0> const& other)
00210 {
00211 resize(other.size());
00212 copy(other);
00213 return *this;
00214 }
00216 AD<T, 0>& operator+=(double other)
00217 {
00218 x_ += other;
00219 return *this;
00220 }
00222 AD<T, 0>& operator+=(AD<T, 0> const& other)
00223 {
00224 resize(other.size());
00225 x_ += other.x_;
00226 for (unsigned int i=0; i < size(); ++i)
00227 dx_[i] += other.dx_[i];
00228 return *this;
00229 }
00231 AD<T, 0>& operator-=(double other)
00232 {
00233 x_ -= other;
00234 return *this;
00235 }
00237 AD<T, 0>& operator-=(AD<T, 0> const& other)
00238 {
00239 if (other.size() > size())
00240 resize(other.size());
00241 x_ -= other.x_;
00242 for (unsigned int i=0; i < size(); ++i)
00243 dx_[i] -= other.dx_[i];
00244 return *this;
00245 }
00247 AD<T, 0>& operator*=(double other)
00248 {
00249 x_ *= other;
00250 for (unsigned int i=0; i < size(); ++i)
00251 dx_[i] *= other;
00252 return *this;
00253 }
00255 AD<T, 0>& operator*=(AD<T, 0> const& other)
00256 {
00257 resize(other.size());
00258 x_ *= other.x_;
00259 for (unsigned int i=0; i < size(); ++i)
00260 dx_[i] = dx_[i]*other.x_ + x_*other.dx_[i];
00261 return *this;
00262 }
00264 AD<T, 0>& operator/=(double other)
00265 {
00266 x_ /= other;
00267 for (unsigned int i=0; i < size(); ++i)
00268 dx_[i] /= other;
00269 return *this;
00270 }
00272 AD<T, 0>& operator/=(AD<T, 0> const& other)
00273 {
00274 resize(other.size());
00275 x_ /= other.x_;
00276 for (unsigned int i=0; i < size(); ++i)
00277 dx_[i] = (dx_[i]*other.x_ - x_*other.dx_[i]) / (other.x_*other.x_);
00278 return *this;
00279 }
00280 void resize(unsigned int n)
00281 {
00282 dx_.resize(n);
00283 }
00284 private:
00285 void zero()
00286 {
00287 for (unsigned int i=0; i < size(); ++i)
00288 dx_[i] = 0.;
00289 }
00290 template<class B>
00291 void copy(AD<B, 0> const& other)
00292 {
00293 if(size() != other.size())
00294 dx_.resize(other.size());
00295 x_ = other.x_;
00296 for (unsigned int i=0; i < size(); ++i)
00297 dx_[i] = (T)other.dx_[i];
00298 }
00299 };
00300
00301 template<class T>
00302 const double mth::AD<T, 0>::_zero_ = 0.;
00303
00304
00305
00306
00307
00309 template <class T, unsigned int N>
00310 AD<T, N> operator-(AD<T, N> const& A)
00311 {
00312 AD<T, N> tmp;
00313 tmp.resize(A.size());
00314 tmp = AD<T, N>(-1.) * A;
00315 return tmp;
00316 }
00317
00318
00319
00320
00321
00323 template <class T, unsigned int N>
00324 AD<T, N> operator+(double L, AD<T, N> const& R)
00325 {
00326 AD<T, N> tmp;
00327 tmp.resize(R.size());
00328 for (unsigned int i=0; i < tmp.size(); ++i)
00329 tmp.dx(i) = T(R.dx(i));
00330 tmp.val() = L + R.val();
00331 return tmp;
00332 }
00333
00335 template <class T, unsigned int N>
00336 AD<T, N> operator+(AD<T, N> const& L, double R)
00337 {
00338 AD<T, N> tmp;
00339 tmp.resize(L.size());
00340 for (unsigned int i=0; i < tmp.size(); ++i)
00341 tmp.dx(i) = T(L.dx(i));
00342 tmp.val() = L.val() + R;
00343 return tmp;
00344 }
00346 template <class T, class B, unsigned int N>
00347 AD<T, N> operator+(AD<T, N> const& L, AD<B, N> const& R)
00348 {
00349 AD<T, N> tmp;
00350 unsigned int max = L.size() > R.size() ? L.size() : R.size();
00351 tmp.resize(max);
00352 for (unsigned int i=0; i < tmp.size(); ++i)
00353 tmp.dx(i) = T(L.dx(i) + R.dx(i));
00354 tmp.val() = L.val() + R.val();
00355 return tmp;
00356 }
00357
00359 template <class T, unsigned int N>
00360 AD<T, N> operator-(double L, AD<T, N> const& R)
00361 {
00362 AD<T, N> tmp;
00363 tmp.resize(R.size());
00364 for (unsigned int i=0; i < tmp.size(); ++i)
00365 tmp.dx(i) = T(-R.dx(i));
00366 tmp.val() = L - R.val();
00367 return tmp;
00368 }
00369
00371 template <class T, unsigned int N>
00372 AD<T, N> operator-(AD<T, N> const& L, double R)
00373 {
00374 AD<T, N> tmp;
00375 tmp.resize(L.size());
00376 for (unsigned int i=0; i < tmp.size(); ++i)
00377 tmp.dx(i) = T(L.dx(i));
00378 tmp.val() = L.val() - R;
00379 return tmp;
00380 }
00381
00383 template <class T, class B, unsigned int N>
00384 AD<T, N> operator-(AD<T, N> const& L, AD<B, N> const& R)
00385 {
00386 AD<T, N> tmp;
00387 unsigned int max = L.size() > R.size() ? L.size() : R.size();
00388 tmp.resize(max);
00389 for (unsigned int i=0; i < tmp.size(); ++i)
00390 tmp.dx(i) = T(L.dx(i) - R.dx(i));
00391 tmp.val() = L.val() - R.val();
00392 return tmp;
00393 }
00394
00396 template <class T, unsigned int N>
00397 AD<T, N> operator*(double L, AD<T, N> const& R)
00398 {
00399 AD<T, N> tmp;
00400 tmp.resize(R.size());
00401 for (unsigned int i=0; i < tmp.size(); ++i)
00402 tmp.dx(i) = T(L * R.dx(i));
00403 tmp.val() = L * R.val();
00404 return tmp;
00405 }
00406
00408 template <class T, unsigned int N>
00409 AD<T, N> operator*(AD<T, N> const& L, double R)
00410 {
00411 AD<T, N> tmp;
00412 tmp.resize(L.size());
00413 for (unsigned int i=0; i < tmp.size(); ++i)
00414 tmp.dx(i) = T(L.dx(i) * R);
00415 tmp.val() = L.val() * R;
00416 return tmp;
00417 }
00418
00420 template <class T, class B, unsigned int N>
00421 AD<T, N> operator*(AD<T, N> const& L, AD<B, N> const& R)
00422 {
00423 AD<T, N> tmp;
00424 unsigned int max = L.size() > R.size() ? L.size() : R.size();
00425 tmp.resize(max);
00426 for (unsigned int i=0; i < tmp.size(); ++i)
00427 tmp.dx(i) = T(L.dx(i) * R + L * R.dx(i));
00428 tmp.val() = L.val() * R.val();
00429 return tmp;
00430 }
00431
00433 template <class T, unsigned int N>
00434 AD<T, N> operator/(double L, AD<T, N> const& R)
00435 {
00436 AD<T, N> tmp;
00437 tmp.resize(R.size());
00438 T R_tmp = R;
00439 for (unsigned int i = 0; i < R.size(); i++)
00440 tmp.dx(i) = T(-L * R.dx(i) * (1. / R_tmp) * (1. / R_tmp));
00441 tmp.val() = L / R.val();
00442 return tmp;
00443 }
00444
00446 template <class T, unsigned int N>
00447 AD<T, N> operator/(AD<T, N> const& L, double R)
00448 {
00449 AD<T, N> tmp;
00450 tmp.resize(L.size());
00451 for (unsigned int i=0; i < L.size(); ++i)
00452 tmp.dx(i) = T(L.dx(i) / R);
00453 tmp.val() = L.val() / R;
00454 return tmp;
00455 }
00456
00458 template <class T, class B, unsigned int N>
00459 AD<T, N> operator/(AD<B, N> const& L, AD<T, N> const& R)
00460 {
00461 AD<T, N> tmp;
00462 unsigned int max = L.size() > R.size() ? L.size() : R.size();
00463 tmp.resize(max);
00464 for (unsigned int i=0; i < tmp.size(); ++i)
00465 tmp.dx(i) = T(((L.dx(i) * R) - (L * R.dx(i))) * (1. / R) * (1. / R));
00466 tmp.val() = L.val() / R.val();
00467 return tmp;
00468 }
00469
00470
00471
00472
00473
00475 double exp(double x)
00476 {
00477 return std::exp(x);
00478 }
00479
00481 template <class T, unsigned int N>
00482 AD<T, N> exp(AD<T, N> const& A)
00483 {
00484 AD<T, N> tmp;
00485 tmp.resize(A.size());
00486 T A_tmp = A;
00487 tmp.val() = std::exp(A.val());
00488 for (unsigned int i=0; i < A.size(); ++i)
00489 tmp.dx(i) = T(A.dx(i) * exp(A_tmp));
00490 return tmp;
00491 }
00492
00494 double log(double A)
00495 {
00496 return std::log(A);
00497 }
00498
00500 template <class T, unsigned int N>
00501 AD<T, N> log(AD<T, N> const& A)
00502 {
00503 AD<T, N> tmp;
00504 tmp.resize(A.size());
00505 T A_tmp = A;
00506 tmp.val() = std::log(A.val());
00507 for (unsigned int i=0; i < A.size(); ++i)
00508 tmp.dx(i) = T(A.dx(i) / A_tmp);
00509 return tmp;
00510 }
00511
00513 double pow(double A, double e)
00514 {
00515 return std::pow(A, e);
00516 }
00517
00519 template <class T, unsigned int N>
00520 AD<T, N> pow(AD<T, N> const& A, const int e)
00521 {
00522 AD<T, N> tmp;
00523 tmp.resize(A.size());
00524 T A_tmp = A;
00525 tmp.val() = std::pow(A.val(), e);
00526 for (unsigned int i=0; i < tmp.size(); ++i)
00527 tmp.dx(i) = T(AD<T, N>(e)*A.dx(i)*pow(A_tmp, (double)e-1.));
00528 return tmp;
00529 }
00530
00532 template <class T, unsigned int N>
00533 AD<T, N> pow(AD<T, N> const& A, const double e)
00534 {
00535 AD<T, N> tmp;
00536 tmp.resize(A.size());
00537 T A_tmp = A;
00538 tmp.val() = pow(A.val(), e);
00539 for (unsigned int i=0; i < tmp.size(); ++i)
00540 tmp.dx(i) = T(e*A.dx(i)*pow(A_tmp, e - 1.));
00541 return tmp;
00542 }
00543
00545 template <class T, unsigned int N>
00546 AD<T, N> pow(const int base, AD<T, N> const& A)
00547 {
00548 AD<T, N> tmp;
00549 tmp.resize(A.size());
00550 T A_tmp = A;
00551 tmp.val() = std::pow((double)base, A.val());
00552 for (unsigned int i=0; i < tmp.size(); ++i)
00553 tmp.dx(i) = T(std::log((double)base) * pow((double)base, A_tmp) * A.dx(i));
00554 return tmp;
00555 }
00556
00558 template <class T, unsigned int N>
00559 AD<T, N> pow(const double base, AD<T, N> const& A)
00560 {
00561 AD<T, N> tmp;
00562 tmp.resize(A.size());
00563 T A_tmp = A;
00564 tmp.val() = std::pow((double)base, A.val());
00565 for (unsigned int i=0; i < tmp.size(); ++i)
00566 tmp.dx(i) = T(std::log(base) * pow(base, A_tmp) * A.dx(i));
00567 return tmp;
00568 }
00569
00570
00572 template <class T, unsigned int N>
00573 AD<T, N> pow(AD<T, N> const& A, AD<T, N> const& e)
00574 {
00575 AD<T, N> tmp;
00576 unsigned int max = A.size() > e.size() ? A.size() : e.size();
00577 tmp.resize(max);
00578 T A_tmp = A;
00579 T e_tmp = e;
00580 tmp.val() = std::pow(A.val(), e.val());
00581 for (unsigned int i=0; i < tmp.size(); ++i)
00582 tmp.dx(i) = T(e.dx(i) * log(A_tmp) * pow(A_tmp, e_tmp) +
00583 e_tmp * A.dx(i) * pow(A_tmp, e_tmp - 1.));
00584 return tmp;
00585 }
00586
00588 double sqrt(double A)
00589 {
00590 return std::sqrt(A);
00591 }
00592
00594 template <class T, unsigned int N>
00595 AD<T, N> sqrt(AD<T, N> const& A)
00596 {
00597 AD<T, N> tmp;
00598 tmp.resize(A.size());
00599 tmp.val() = std::sqrt(A.val());
00600 T A_tmp = A;
00601 for (unsigned int i=0; i < tmp.size(); ++i)
00602 tmp.dx(i) = T(.5 * A.dx(i) / sqrt(A_tmp));
00603 return tmp;
00604 }
00605
00607 double sin(double A)
00608 {
00609 return std::sin(A);
00610 }
00611
00613 double cos(double A)
00614 {
00615 return std::cos(A);
00616 }
00617
00619 template <class T, unsigned int N>
00620 AD<T, N> sin(AD<T, N> const& A)
00621 {
00622 AD<T, N> tmp;
00623 tmp.resize(A.size());
00624 tmp.val() = std::sin(A.val());
00625 T A_tmp = A;
00626 for(unsigned int i = 0; i < tmp.size(); i++)
00627 tmp.dx(i) = T(cos(A_tmp) * A.dx(i));
00628 return tmp;
00629 }
00630
00632 template <class T, unsigned int N>
00633 AD<T, N> cos(AD<T, N> const& A)
00634 {
00635 AD<T, N> tmp;
00636 tmp.val() = std::cos(A.val());
00637 T A_tmp = A;
00638 tmp.resize(A.size());
00639 for(unsigned int i = 0; i < tmp.size(); i++)
00640 tmp.dx(i) = T(-sin(A_tmp) * A.dx(i));
00641 return tmp;
00642 }
00643
00645 template <class T, unsigned int N>
00646 AD<T, N> tan(AD<T, N> const& A)
00647 {
00648 AD<T, N> tmp;
00649 tmp.val() = std::tan(A.val());
00650 T A_tmp = A;
00651 tmp.resize(A.size());
00652 for(unsigned int i = 0; i < tmp.size(); i++)
00653 tmp.dx(i) = T(A.dx(i) * (1. / (cos(A_tmp) * cos(A_tmp))));
00654 return tmp;
00655 }
00656
00658 template <class T, unsigned int N>
00659 AD<T, N> abs(AD<T, N> const& A)
00660 {
00661 int sign = A.val() > 0 ? 1 : 0;
00662 if (sign) return A;
00663 else return (-A);
00664 }
00665
00666
00667
00668
00669
00671 template <class T, unsigned int N>
00672 bool operator<(double L, AD<T, N> const& R)
00673 {
00674 return L < R.val();
00675 }
00676
00678 template <class T, unsigned int N>
00679 bool operator<(AD<T, N> const& R, double L)
00680 {
00681 return R.val() < L;
00682 }
00683
00685 template <class T, unsigned int N>
00686 bool operator<(AD<T, N> const& R, AD<T, N> const& L)
00687 {
00688 return R.val() < L.val();
00689 }
00690
00692 template <class T, unsigned int N>
00693 bool operator<=(double L, AD<T, N> const& R)
00694 {
00695 return L <= R.val();
00696 }
00697
00699 template <class T, unsigned int N>
00700 bool operator<=(AD<T, N> const& R, double L)
00701 {
00702 return R.val() <= L;
00703 }
00704
00706 template <class T, unsigned int N>
00707 bool operator<=(AD<T, N> const& R, AD<T, N> const& L)
00708 {
00709 return R.val() <= L.val();
00710 }
00711
00713 template <class T, unsigned int N>
00714 bool operator>(double L, AD<T, N> const& R)
00715 {
00716 return L > R.val();
00717 }
00718
00720 template <class T, unsigned int N>
00721 bool operator>(AD<T, N> const& R, double L)
00722 {
00723 return R.val() > L;
00724 }
00725
00727 template <class T, unsigned int N>
00728 bool operator>(AD<T, N> const& R, AD<T, N> const& L)
00729 {
00730 return R.val() > L.val();
00731 }
00732
00734 template <class T, unsigned int N>
00735 bool operator>=(double L, AD<T, N> const& R)
00736 {
00737 return L >= R.val();
00738 }
00739
00741 template <class T, unsigned int N>
00742 bool operator>=(AD<T, N> const& R, double L)
00743 {
00744 return R.val() >= L;
00745 }
00746
00748 template <class T, unsigned int N>
00749 bool operator>=(AD<T, N> const& R, AD<T, N> const& L)
00750 {
00751 return R.val() >= L.val();
00752 }
00753 }
00754 #endif