Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Sacado_Tay_TaylorImp.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Sacado Package
4 //
5 // Copyright 2006 NTESS and the Sacado contributors.
6 // SPDX-License-Identifier: LGPL-2.1-or-later
7 // *****************************************************************************
8 // @HEADER
9 
10 // disable clang warnings
11 #if defined (__clang__) && !defined (__INTEL_COMPILER)
12 #pragma clang system_header
13 #endif
14 
15 #include "Sacado_ConfigDefs.h"
17 #include <ostream> // for std::ostream
18 
19 namespace Sacado {
20 namespace Tay {
21 
22 template <typename T>
25  coeff_(NULL), deg_(-1), len_(0)
26 {
27 }
28 
29 template <typename T>
31 TaylorData(const T& x) :
32  coeff_(), deg_(0), len_(1)
33 {
35  coeff_[0] = x;
36 }
37 
38 template <typename T>
40 TaylorData(int d, const T& x) :
41  coeff_(), deg_(d), len_(d+1)
42 {
44  coeff_[0] = x;
45 }
46 
47 template <typename T>
49 TaylorData(int d) :
50  coeff_(), deg_(d), len_(d+1)
51 {
53 }
54 
55 template <typename T>
57 TaylorData(int d, int l) :
58  coeff_(), deg_(d), len_(l)
59 {
61 }
62 
63 template <typename T>
65 TaylorData(const typename Taylor<T>::TaylorData& x) :
66  coeff_(), deg_(x.deg_), len_(x.deg_+1)
67 {
69 }
70 
71 template <typename T>
74 {
75  if (len_ > 0)
77 }
78 
79 template <typename T>
80 typename Taylor<T>::TaylorData&
83 {
84  if (len_ < x.deg_+1) {
86  len_ = x.deg_+1;
87  deg_ = x.deg_;
88  coeff_ = Sacado::ds_array<T>::get_and_fill(x.coeff_, len_);
89  }
90  else {
91  deg_ = x.deg_;
92  Sacado::ds_array<T>::copy(x.coeff_, coeff_, deg_+1);
93  }
94 
95  return *this;
96 }
97 
98 template <typename T>
101  th(new TaylorData)
102 {
103 }
104 
105 template <typename T>
107 Taylor(const T& x) :
108  th(new TaylorData(x))
109 {
110 }
111 
112 template <typename T>
115  th(new TaylorData(value_type(x)))
116 {
117 }
118 
119 template <typename T>
121 Taylor(int d, const T& x) :
122  th(new TaylorData(d, x))
123 {
124 }
125 
126 template <typename T>
128 Taylor(int d) :
129  th(new TaylorData(d))
130 {
131 }
132 
133 template <typename T>
135 Taylor(const Taylor<T>& x) :
136  th(x.th)
137 {
138 }
139 
140 template <typename T>
143 {
144 }
145 
146 template <typename T>
147 void
149 resize(int d, bool keep_coeffs)
150 {
151  if (d+1 > length()) {
153  if (keep_coeffs)
154  Sacado::ds_array<T>::copy(th->coeff_, h->coeff_, th->deg_+1);
155  th = h;
156  }
157  else {
158  th.makeOwnCopy();
159  if (!keep_coeffs)
160  Sacado::ds_array<T>::zero(coeff(), degree()+1);
161  th->deg_ = d;
162  }
163 }
164 
165 template <typename T>
166 void
168 reserve(int d)
169 {
170  if (d+1 > length()) {
171  Sacado::Handle<TaylorData> h(new TaylorData(th->deg_,d+1));
172  Sacado::ds_array<T>::copy(th->coeff_, h->coeff_, th->deg_+1);
173  th = h;
174  }
175 }
176 
177 template <typename T>
178 Taylor<T>&
180 operator=(const T& v)
181 {
182  th.makeOwnCopy();
183 
184  if (th->len_ == 0) {
185  th->len_ = 1;
186  th->deg_ = 0;
187  th->coeff_ = Sacado::ds_array<T>::get_and_fill(th->len_);
188  }
189 
190  th->coeff_[0] = v;
191  Sacado::ds_array<T>::zero(th->coeff_+1, th->deg_);
192 
193  return *this;
194 }
195 
196 template <typename T>
197 Taylor<T>&
200 {
201  return operator=(value_type(v));
202 }
203 
204 template <typename T>
205 Taylor<T>&
208 {
209  th = x.th;
210  return *this;
211 }
212 
213 template <typename T>
214 Taylor<T>
216 operator+() const
217 {
218  return *this;
219 }
220 
221 template <typename T>
222 Taylor<T>
224 operator-() const
225 {
226  Taylor<T> x;
227  x.th->deg_ = th->deg_;
228  x.th->len_ = th->deg_+1;
229  x.th->coeff_ = Sacado::ds_array<T>::get_and_fill(x.th->len_);
230  for (int i=0; i<=th->deg_; i++)
231  x.th->coeff_[i] = -th->coeff_[i];
232 
233  return x;
234 }
235 
236 template <typename T>
237  Taylor<T>&
239 operator+=(const T& v)
240 {
241  th.makeOwnCopy();
242 
243  th->coeff_[0] += v;
244 
245  return *this;
246 }
247 
248 template <typename T>
249 Taylor<T>&
251 operator-=(const T& v)
252 {
253  th.makeOwnCopy();
254 
255  th->coeff_[0] -= v;
256 
257  return *this;
258 }
259 
260 template <typename T>
261 Taylor<T>&
263 operator*=(const T& v)
264 {
265  th.makeOwnCopy();
266 
267  for (int i=0; i<=th->deg_; i++)
268  th->coeff_[i] *= v;
269 
270  return *this;
271 }
272 
273 template <typename T>
274 Taylor<T>&
276 operator/=(const T& v)
277 {
278  th.makeOwnCopy();
279 
280  for (int i=0; i<=th->deg_; i++)
281  th->coeff_[i] /= v;
282 
283  return *this;
284 }
285 
286 template <typename T>
287 Taylor<T>&
290 {
291  th.makeOwnCopy();
292 
293  int d = degree();
294  int xd = x.degree();
295  int dmin = xd < d ? xd : d;
296 
297  int l = xd + 1;
298  bool need_resize = l > length();
299  T* c;
300  if (need_resize) {
302  for (int i=0; i<=d; i++)
303  c[i] = fastAccessCoeff(i);
304  }
305  else
306  c = th->coeff_;
307  T* xc = x.th->coeff_;
308 
309  for (int i=0; i<=dmin; i++)
310  c[i] += xc[i];
311  if (th->deg_ < xd) {
312  for (int i=d+1; i<=xd; i++)
313  c[i] = xc[i];
314  th->deg_ = xd;
315  }
316 
317  if (need_resize) {
318  Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
319  th->len_ = l;
320  th->coeff_ = c;
321  }
322 
323  return *this;
324 }
325 
326 template <typename T>
327 Taylor<T>&
330 {
331  th.makeOwnCopy();
332 
333  int d = degree();
334  int xd = x.degree();
335  int dmin = xd < d ? xd : d;
336 
337  int l = xd + 1;
338  bool need_resize = l > length();
339  T* c;
340  if (need_resize) {
342  for (int i=0; i<=d; i++)
343  c[i] = fastAccessCoeff(i);
344  }
345  else
346  c = th->coeff_;
347  T* xc = x.th->coeff_;
348 
349  for (int i=0; i<=dmin; i++)
350  c[i] -= xc[i];
351  if (d < xd) {
352  for (int i=d+1; i<=xd; i++)
353  c[i] = -xc[i];
354  th->deg_ = xd;
355  }
356 
357  if (need_resize) {
358  Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
359  th->len_ = l;
360  th->coeff_ = c;
361  }
362 
363  return *this;
364 }
365 
366 template <typename T>
367 Taylor<T>&
370 {
371  int d = degree();
372  int xd = x.degree();
373 
374 #ifdef SACADO_DEBUG
375  if ((xd != d) && (xd != 0) && (d != 0))
376  throw "Taylor Error: Attempt to assign with incompatible degrees";
377 #endif
378 
379  T* c = th->coeff_;
380  T* xc = x.th->coeff_;
381 
382  if (d > 0 && xd > 0) {
384  T* cc = h->coeff_;
385  T tmp;
386  for (int i=0; i<=d; i++) {
387  tmp = T(0.0);
388  for (int k=0; k<=i; ++k)
389  tmp += c[k]*xc[i-k];
390  cc[i] = tmp;
391  }
392  th = h;
393  }
394  else if (d > 0) {
395  th.makeOwnCopy();
396  for (int i=0; i<=d; i++)
397  c[i] = c[i]*xc[0];
398  }
399  else if (xd >= 0) {
400  if (length() < xd+1) {
402  T* cc = h->coeff_;
403  for (int i=0; i<=xd; i++)
404  cc[i] = c[0]*xc[i];
405  th = h;
406  }
407  else {
408  th.makeOwnCopy();
409  for (int i=d; i>=0; i--)
410  c[i] = c[0]*xc[i];
411  }
412  }
413 
414  return *this;
415 }
416 
417 template <typename T>
418 Taylor<T>&
421 {
422  int d = degree();
423  int xd = x.degree();
424 
425 #ifdef SACADO_DEBUG
426  if ((xd != d) && (xd != 0) && (d != 0))
427  throw "Taylor Error: Attempt to assign with incompatible degrees";
428 #endif
429 
430  T* c = th->coeff_;
431  T* xc = x.th->coeff_;
432 
433  if (d > 0 && xd > 0) {
435  T* cc = h->coeff_;
436  T tmp;
437  for(int i=0; i<=d; i++) {
438  tmp = c[i];
439  for (int k=1; k<=i; ++k)
440  tmp -= xc[k]*cc[i-k];
441  cc[i] = tmp / xc[0];
442  }
443  th = h;
444  }
445  else if (d > 0) {
446  th.makeOwnCopy();
447  for (int i=0; i<=d; i++)
448  c[i] = c[i]/xc[0];
449  }
450  else if (xd >= 0) {
452  T* cc = h->coeff_;
453  T tmp;
454  cc[0] = c[0] / xc[0];
455  for (int i=1; i<=xd; i++) {
456  tmp = T(0.0);
457  for (int k=1; k<=i; ++k)
458  tmp -= xc[k]*cc[i-k];
459  cc[i] = tmp / xc[0];
460  }
461  th = h;
462  }
463 
464  return *this;
465 }
466 
467 template <typename T>
468 void
470 resizeCoeffs(int len)
471 {
472  if (th->coeff_)
473  Sacado::ds_array<T>::destroy_and_release(th->coeff_, th->len_);
474  th->len_ = len;
475  th->coeff_ = Sacado::ds_array<T>::get_and_fill(th->len_);
476 }
477 
478 template <typename T>
479 Taylor<T>
480 operator+(const Base< Taylor<T> >& aa,
481  const Base< Taylor<T> >& bb)
482 {
483  const Taylor<T>& a = aa.derived();
484  const Taylor<T>& b = bb.derived();
485 
486  int da = a.degree();
487  int db = b.degree();
488  int dc = da > db ? da : db;
489 
490 #ifdef SACADO_DEBUG
491  if ((da != db) && (da != 0) && (db != 0))
492  throw "operator+(): Arguments have incompatible degrees!";
493 #endif
494 
495  Taylor<T> c(dc);
496  const T* ca = a.coeff();
497  const T* cb = b.coeff();
498  T* cc = c.coeff();
499 
500  if (da > 0 && db > 0) {
501  for (int i=0; i<=dc; i++)
502  cc[i] = ca[i] + cb[i];
503  }
504  else if (da > 0) {
505  cc[0] = ca[0] + cb[0];
506  for (int i=1; i<=dc; i++)
507  cc[i] = ca[i];
508  }
509  else if (db >= 0) {
510  cc[0] = ca[0] + cb[0];
511  for (int i=1; i<=dc; i++)
512  cc[i] = cb[i];
513  }
514 
515  return c;
516 }
517 
518 template <typename T>
519 Taylor<T>
521  const Base< Taylor<T> >& bb)
522 {
523  const Taylor<T>& b = bb.derived();
524 
525  int dc = b.degree();
526 
527  Taylor<T> c(dc);
528  const T* cb = b.coeff();
529  T* cc = c.coeff();
530 
531  cc[0] = a + cb[0];
532  for (int i=1; i<=dc; i++)
533  cc[i] = cb[i];
534 
535  return c;
536 }
537 
538 template <typename T>
539 Taylor<T>
540 operator+(const Base< Taylor<T> >& aa,
541  const typename Taylor<T>::value_type& b)
542 {
543  const Taylor<T>& a = aa.derived();
544 
545  int dc = a.degree();
546 
547  Taylor<T> c(dc);
548  const T* ca = a.coeff();
549  T* cc = c.coeff();
550 
551  cc[0] = ca[0] + b;
552  for (int i=1; i<=dc; i++)
553  cc[i] = ca[i];
554 
555  return c;
556 }
557 
558 template <typename T>
559 Taylor<T>
560 operator-(const Base< Taylor<T> >& aa,
561  const Base< Taylor<T> >& bb)
562 {
563  const Taylor<T>& a = aa.derived();
564  const Taylor<T>& b = bb.derived();
565 
566  int da = a.degree();
567  int db = b.degree();
568  int dc = da > db ? da : db;
569 
570 #ifdef SACADO_DEBUG
571  if ((da != db) && (da != 0) && (db != 0))
572  throw "operator+(): Arguments have incompatible degrees!";
573 #endif
574 
575  Taylor<T> c(dc);
576  const T* ca = a.coeff();
577  const T* cb = b.coeff();
578  T* cc = c.coeff();
579 
580  if (da > 0 && db > 0) {
581  for (int i=0; i<=dc; i++)
582  cc[i] = ca[i] - cb[i];
583  }
584  else if (da > 0) {
585  cc[0] = ca[0] - cb[0];
586  for (int i=1; i<=dc; i++)
587  cc[i] = ca[i];
588  }
589  else if (db >= 0) {
590  cc[0] = ca[0] - cb[0];
591  for (int i=1; i<=dc; i++)
592  cc[i] = -cb[i];
593  }
594 
595  return c;
596 }
597 
598 template <typename T>
599 Taylor<T>
601  const Base< Taylor<T> >& bb)
602 {
603  const Taylor<T>& b = bb.derived();
604 
605  int dc = b.degree();
606 
607  Taylor<T> c(dc);
608  const T* cb = b.coeff();
609  T* cc = c.coeff();
610 
611  cc[0] = a - cb[0];
612  for (int i=1; i<=dc; i++)
613  cc[i] = -cb[i];
614 
615  return c;
616 }
617 
618 template <typename T>
619 Taylor<T>
620 operator-(const Base< Taylor<T> >& aa,
621  const typename Taylor<T>::value_type& b)
622 {
623  const Taylor<T>& a = aa.derived();
624 
625  int dc = a.degree();
626 
627  Taylor<T> c(dc);
628  const T* ca = a.coeff();
629  T* cc = c.coeff();
630 
631  cc[0] = ca[0] - b;
632  for (int i=1; i<=dc; i++)
633  cc[i] = ca[i];
634 
635  return c;
636 }
637 
638 template <typename T>
639 Taylor<T>
640 operator*(const Base< Taylor<T> >& aa,
641  const Base< Taylor<T> >& bb)
642 {
643  const Taylor<T>& a = aa.derived();
644  const Taylor<T>& b = bb.derived();
645 
646  int da = a.degree();
647  int db = b.degree();
648  int dc = da > db ? da : db;
649 
650 #ifdef SACADO_DEBUG
651  if ((da != db) && (da != 0) && (db != 0))
652  throw "operator+(): Arguments have incompatible degrees!";
653 #endif
654 
655  Taylor<T> c(dc);
656  const T* ca = a.coeff();
657  const T* cb = b.coeff();
658  T* cc = c.coeff();
659 
660  if (da > 0 && db > 0) {
661  T tmp;
662  for (int i=0; i<=dc; i++) {
663  tmp = T(0.0);
664  for (int k=0; k<=i; k++)
665  tmp += ca[k]*cb[i-k];
666  cc[i] = tmp;
667  }
668  }
669  else if (da > 0) {
670  for (int i=0; i<=dc; i++)
671  cc[i] = ca[i]*cb[0];
672  }
673  else if (db >= 0) {
674  for (int i=0; i<=dc; i++)
675  cc[i] = ca[0]*cb[i];
676  }
677 
678  return c;
679 }
680 
681 template <typename T>
682 Taylor<T>
684  const Base< Taylor<T> >& bb)
685 {
686  const Taylor<T>& b = bb.derived();
687 
688  int dc = b.degree();
689 
690  Taylor<T> c(dc);
691  const T* cb = b.coeff();
692  T* cc = c.coeff();
693 
694  for (int i=0; i<=dc; i++)
695  cc[i] = a*cb[i];
696 
697  return c;
698 }
699 
700 template <typename T>
701 Taylor<T>
702 operator*(const Base< Taylor<T> >& aa,
703  const typename Taylor<T>::value_type& b)
704 {
705  const Taylor<T>& a = aa.derived();
706 
707  int dc = a.degree();
708 
709  Taylor<T> c(dc);
710  const T* ca = a.coeff();
711  T* cc = c.coeff();
712 
713  for (int i=0; i<=dc; i++)
714  cc[i] = ca[i]*b;
715 
716  return c;
717 }
718 
719 template <typename T>
720 Taylor<T>
721 operator/(const Base< Taylor<T> >& aa,
722  const Base< Taylor<T> >& bb)
723 {
724  const Taylor<T>& a = aa.derived();
725  const Taylor<T>& b = bb.derived();
726 
727  int da = a.degree();
728  int db = b.degree();
729  int dc = da > db ? da : db;
730 
731 #ifdef SACADO_DEBUG
732  if ((da != db) && (da != 0) && (db != 0))
733  throw "operator+(): Arguments have incompatible degrees!";
734 #endif
735 
736  Taylor<T> c(dc);
737  const T* ca = a.coeff();
738  const T* cb = b.coeff();
739  T* cc = c.coeff();
740 
741  if (da > 0 && db > 0) {
742  T tmp;
743  for (int i=0; i<=dc; i++) {
744  tmp = ca[i];
745  for (int k=0; k<=i; k++)
746  tmp -= cb[k]*cc[i-k];
747  cc[i] = tmp / cb[0];
748  }
749  }
750  else if (da > 0) {
751  for (int i=0; i<=dc; i++)
752  cc[i] = ca[i]/cb[0];
753  }
754  else if (db >= 0) {
755  T tmp;
756  cc[0] = ca[0] / cb[0];
757  for (int i=1; i<=dc; i++) {
758  tmp = T(0.0);
759  for (int k=0; k<=i; k++)
760  tmp -= cb[k]*cc[i-k];
761  cc[i] = tmp / cb[0];
762  }
763  }
764 
765  return c;
766 }
767 
768 template <typename T>
769 Taylor<T>
771  const Base< Taylor<T> >& bb)
772 {
773  const Taylor<T>& b = bb.derived();
774 
775  int dc = b.degree();
776 
777  Taylor<T> c(dc);
778  const T* cb = b.coeff();
779  T* cc = c.coeff();
780 
781  T tmp;
782  cc[0] = a / cb[0];
783  for (int i=1; i<=dc; i++) {
784  tmp = T(0.0);
785  for (int k=0; k<=i; k++)
786  tmp -= cb[k]*cc[i-k];
787  cc[i] = tmp / cb[0];
788  }
789 
790  return c;
791 }
792 
793 template <typename T>
794 Taylor<T>
795 operator/(const Base< Taylor<T> >& aa,
796  const typename Taylor<T>::value_type& b)
797 {
798  const Taylor<T>& a = aa.derived();
799 
800  int dc = a.degree();
801 
802  Taylor<T> c(dc);
803  const T* ca = a.coeff();
804  T* cc = c.coeff();
805 
806  for (int i=0; i<=dc; i++)
807  cc[i] = ca[i]/b;
808 
809  return c;
810 }
811 
812 template <typename T>
813 Taylor<T>
814 exp(const Base< Taylor<T> >& aa)
815 {
816  const Taylor<T>& a = aa.derived();
817  int dc = a.degree();
818 
819  Taylor<T> c(dc);
820  const T* ca = a.coeff();
821  T* cc = c.coeff();
822 
823  T tmp;
824  cc[0] = std::exp(ca[0]);
825  for (int i=1; i<=dc; i++) {
826  tmp = T(0.0);
827  for (int k=1; k<=i; k++)
828  tmp += k*cc[i-k]*ca[k];
829  cc[i] = tmp / i;
830  }
831 
832  return c;
833 }
834 
835 template <typename T>
836 Taylor<T>
837 log(const Base< Taylor<T> >& aa)
838 {
839  const Taylor<T>& a = aa.derived();
840  int dc = a.degree();
841 
842  Taylor<T> c(dc);
843  const T* ca = a.coeff();
844  T* cc = c.coeff();
845 
846  T tmp;
847  cc[0] = std::log(ca[0]);
848  for (int i=1; i<=dc; i++) {
849  tmp = i*ca[i];
850  for (int k=1; k<=i-1; k++)
851  tmp -= k*ca[i-k]*cc[k];
852  cc[i] = tmp / (i*ca[0]);
853  }
854 
855  return c;
856 }
857 
858 template <typename T>
859 Taylor<T>
860 log10(const Base< Taylor<T> >& aa)
861 {
862  const Taylor<T>& a = aa.derived();
863  return log(a) / std::log(10.0);
864 }
865 
866 template <typename T>
867 Taylor<T>
868 sqrt(const Base< Taylor<T> >& aa)
869 {
870  const Taylor<T>& a = aa.derived();
871  int dc = a.degree();
872 
873  Taylor<T> c(dc);
874  const T* ca = a.coeff();
875  T* cc = c.coeff();
876 
877  T tmp;
878  cc[0] = std::sqrt(ca[0]);
879  for (int i=1; i<=dc; i++) {
880  tmp = ca[i];
881  for (int k=1; k<=i-1; k++)
882  tmp -= cc[k]*cc[i-k];
883  cc[i] = tmp / (2.0*cc[0]);
884  }
885 
886  return c;
887 }
888 
889 template <typename T>
890 Taylor<T>
891 cbrt(const Base< Taylor<T> >& aa)
892 {
893  return pow(aa, typename Taylor<T>::value_type(1.0/3.0));
894 }
895 
896 template <typename T>
897 Taylor<T>
898 pow(const Base< Taylor<T> >& aa,
899  const Base< Taylor<T> >& bb)
900 {
901  const Taylor<T>& a = aa.derived();
902  const Taylor<T>& b = bb.derived();
903  return exp(b*log(a));
904 }
905 
906 template <typename T>
907 Taylor<T>
908 pow(const typename Taylor<T>::value_type& a,
909  const Base< Taylor<T> >& bb)
910 {
911  const Taylor<T>& b = bb.derived();
912  return exp(b*std::log(a));
913 }
914 
915 template <typename T>
916 Taylor<T>
917 pow(const Base< Taylor<T> >& aa,
918  const typename Taylor<T>::value_type& b)
919 {
920  const Taylor<T>& a = aa.derived();
921  return exp(b*log(a));
922 }
923 
924 template <typename T>
925 void
926 sincos(const Base< Taylor<T> >& aa,
927  Taylor<T>& s,
928  Taylor<T>& c)
929 {
930  const Taylor<T>& a = aa.derived();
931  int dc = a.degree();
932  if (s.degree() != dc)
933  s.resize(dc, false);
934  if (c.degree() != dc)
935  c.resize(dc, false);
936 
937  const T* ca = a.coeff();
938  T* cs = s.coeff();
939  T* cc = c.coeff();
940 
941  T tmp1;
942  T tmp2;
943  cs[0] = std::sin(ca[0]);
944  cc[0] = std::cos(ca[0]);
945  for (int i=1; i<=dc; i++) {
946  tmp1 = T(0.0);
947  tmp2 = T(0.0);
948  for (int k=1; k<=i; k++) {
949  tmp1 += k*ca[k]*cc[i-k];
950  tmp2 -= k*ca[k]*cs[i-k];
951  }
952  cs[i] = tmp1 / i;
953  cc[i] = tmp2 / i;
954  }
955 }
956 
957 template <typename T>
958 Taylor<T>
959 sin(const Base< Taylor<T> >& aa)
960 {
961  const Taylor<T>& a = aa.derived();
962  int dc = a.degree();
963  Taylor<T> s(dc);
964  Taylor<T> c(dc);
965  sincos(a, s, c);
966 
967  return s;
968 }
969 
970 template <typename T>
971 Taylor<T>
972 cos(const Base< Taylor<T> >& aa)
973 {
974  const Taylor<T>& a = aa.derived();
975  int dc = a.degree();
976  Taylor<T> s(dc);
977  Taylor<T> c(dc);
978  sincos(a, s, c);
979 
980  return c;
981 }
982 
983 template <typename T>
984 Taylor<T>
985 tan(const Base< Taylor<T> >& aa)
986 {
987  const Taylor<T>& a = aa.derived();
988  int dc = a.degree();
989  Taylor<T> s(dc);
990  Taylor<T> c(dc);
991 
992  sincos(a, s, c);
993 
994  return s / c;
995 }
996 
997 template <typename T>
998 void
999 sinhcosh(const Base< Taylor<T> >& aa,
1000  Taylor<T>& s,
1001  Taylor<T>& c)
1002 {
1003  const Taylor<T>& a = aa.derived();
1004  int dc = a.degree();
1005  if (s.degree() != dc)
1006  s.resize(dc, false);
1007  if (c.degree() != dc)
1008  c.resize(dc, false);
1009 
1010  const T* ca = a.coeff();
1011  T* cs = s.coeff();
1012  T* cc = c.coeff();
1013 
1014  T tmp1;
1015  T tmp2;
1016  cs[0] = std::sinh(ca[0]);
1017  cc[0] = std::cosh(ca[0]);
1018  for (int i=1; i<=dc; i++) {
1019  tmp1 = T(0.0);
1020  tmp2 = T(0.0);
1021  for (int k=1; k<=i; k++) {
1022  tmp1 += k*ca[k]*cc[i-k];
1023  tmp2 += k*ca[k]*cs[i-k];
1024  }
1025  cs[i] = tmp1 / i;
1026  cc[i] = tmp2 / i;
1027  }
1028 }
1029 
1030 template <typename T>
1031 Taylor<T>
1032 sinh(const Base< Taylor<T> >& aa)
1033 {
1034  const Taylor<T>& a = aa.derived();
1035  int dc = a.degree();
1036  Taylor<T> s(dc);
1037  Taylor<T> c(dc);
1038  sinhcosh(a, s, c);
1039 
1040  return s;
1041 }
1042 
1043 template <typename T>
1044 Taylor<T>
1045 cosh(const Base< Taylor<T> >& aa)
1046 {
1047  const Taylor<T>& a = aa.derived();
1048  int dc = a.degree();
1049  Taylor<T> s(dc);
1050  Taylor<T> c(dc);
1051  sinhcosh(a, s, c);
1052 
1053  return c;
1054 }
1055 
1056 template <typename T>
1057 Taylor<T>
1058 tanh(const Base< Taylor<T> >& aa)
1059 {
1060  const Taylor<T>& a = aa.derived();
1061  int dc = a.degree();
1062  Taylor<T> s(dc);
1063  Taylor<T> c(dc);
1064 
1065  sinhcosh(a, s, c);
1066 
1067  return s / c;
1068 }
1069 
1070 template <typename T>
1071 Taylor<T>
1072 quad(const typename Taylor<T>::value_type& c0,
1073  const Base< Taylor<T> >& aa,
1074  const Base< Taylor<T> >& bb)
1075 {
1076  const Taylor<T>& a = aa.derived();
1077  const Taylor<T>& b = bb.derived();
1078  int dc = a.degree();
1079 
1080  Taylor<T> c(dc);
1081  const T* ca = a.coeff();
1082  const T* cb = b.coeff();
1083  T* cc = c.coeff();
1084 
1085  T tmp;
1086  cc[0] = c0;
1087  for (int i=1; i<=dc; i++) {
1088  tmp = T(0.0);
1089  for (int k=1; k<=i; k++)
1090  tmp += k*ca[k]*cb[i-k];
1091  cc[i] = tmp / i;
1092  }
1093 
1094  return c;
1095 }
1096 
1097 template <typename T>
1098 Taylor<T>
1099 acos(const Base< Taylor<T> >& aa)
1100 {
1101  const Taylor<T>& a = aa.derived();
1102  Taylor<T> b = -1.0 / sqrt(1.0 - a*a);
1103  return quad(std::acos(a.coeff(0)), a, b);
1104 }
1105 
1106 template <typename T>
1107 Taylor<T>
1108 asin(const Base< Taylor<T> >& aa)
1109 {
1110  const Taylor<T>& a = aa.derived();
1111  Taylor<T> b = 1.0 / sqrt(1.0 - a*a);
1112  return quad(std::asin(a.coeff(0)), a, b);
1113 }
1114 
1115 template <typename T>
1116 Taylor<T>
1117 atan(const Base< Taylor<T> >& aa)
1118 {
1119  const Taylor<T>& a = aa.derived();
1120  Taylor<T> b = 1.0 / (1.0 + a*a);
1121  return quad(std::atan(a.coeff(0)), a, b);
1122 }
1123 
1124 template <typename T>
1125 Taylor<T>
1126 atan2(const Base< Taylor<T> >& aa,
1127  const Base< Taylor<T> >& bb)
1128 {
1129  const Taylor<T>& a = aa.derived();
1130  const Taylor<T>& b = bb.derived();
1131 
1132  Taylor<T> c = atan(a/b);
1133  c.fastAccessCoeff(0) = atan2(a.coeff(0),b.coeff(0));
1134  return c;
1135 }
1136 
1137 template <typename T>
1138 Taylor<T>
1140  const Base< Taylor<T> >& bb)
1141 {
1142  const Taylor<T>& b = bb.derived();
1143 
1144  Taylor<T> c = atan(a/b);
1145  c.fastAccessCoeff(0) = atan2(a,b.coeff(0));
1146  return c;
1147 }
1148 
1149 template <typename T>
1150 Taylor<T>
1151 atan2(const Base< Taylor<T> >& aa,
1152  const typename Taylor<T>::value_type& b)
1153 {
1154  const Taylor<T>& a = aa.derived();
1155 
1156  Taylor<T> c = atan(a/b);
1157  c.fastAccessCoeff(0) = atan2(a.coeff(0),b);
1158  return c;
1159 }
1160 
1161 template <typename T>
1162 Taylor<T>
1163 acosh(const Base< Taylor<T> >& aa)
1164 {
1165  const Taylor<T>& a = aa.derived();
1166  Taylor<T> b = -1.0 / sqrt(1.0 - a*a);
1167  return quad(acosh(a.coeff(0)), a, b);
1168 }
1169 
1170 template <typename T>
1171 Taylor<T>
1172 asinh(const Base< Taylor<T> >& aa)
1173 {
1174  const Taylor<T>& a = aa.derived();
1175  Taylor<T> b = 1.0 / sqrt(a*a - 1.0);
1176  return quad(asinh(a.coeff(0)), a, b);
1177 }
1178 
1179 template <typename T>
1180 Taylor<T>
1181 atanh(const Base< Taylor<T> >& aa)
1182 {
1183  const Taylor<T>& a = aa.derived();
1184  Taylor<T> b = 1.0 / (1.0 - a*a);
1185  return quad(atanh(a.coeff(0)), a, b);
1186 }
1187 
1188 template <typename T>
1189 Taylor<T>
1190 fabs(const Base< Taylor<T> >& aa)
1191 {
1192  const Taylor<T>& a = aa.derived();
1193  if (a.coeff(0) >= 0)
1194  return a;
1195  else
1196  return -a;
1197 }
1198 
1199 template <typename T>
1200 Taylor<T>
1201 abs(const Base< Taylor<T> >& aa)
1202 {
1203  const Taylor<T>& a = aa.derived();
1204  if (a.coeff(0) >= 0)
1205  return a;
1206  else
1207  return -a;
1208 }
1209 
1210 template <typename T>
1211 Taylor<T>
1212 max(const Base< Taylor<T> >& aa,
1213  const Base< Taylor<T> >& bb)
1214 {
1215  const Taylor<T>& a = aa.derived();
1216  const Taylor<T>& b = bb.derived();
1217 
1218  if (a.coeff(0) >= b.coeff(0))
1219  return a;
1220  else
1221  return b;
1222 }
1223 
1224 template <typename T>
1225 Taylor<T>
1226 max(const typename Taylor<T>::value_type& a,
1227  const Base< Taylor<T> >& bb)
1228 {
1229  const Taylor<T>& b = bb.derived();
1230 
1231  if (a >= b.coeff(0))
1232  return Taylor<T>(b.degree(), a);
1233  else
1234  return b;
1235 }
1236 
1237 template <typename T>
1238 Taylor<T>
1239 max(const Base< Taylor<T> >& aa,
1240  const typename Taylor<T>::value_type& b)
1241 {
1242  const Taylor<T>& a = aa.derived();
1243 
1244  if (a.coeff(0) >= b)
1245  return a;
1246  else
1247  return Taylor<T>(a.degree(), b);
1248 }
1249 
1250 template <typename T>
1251 Taylor<T>
1252 min(const Base< Taylor<T> >& aa,
1253  const Base< Taylor<T> >& bb)
1254 {
1255  const Taylor<T>& a = aa.derived();
1256  const Taylor<T>& b = bb.derived();
1257 
1258  if (a.coeff(0) <= b.coeff(0))
1259  return a;
1260  else
1261  return b;
1262 }
1263 
1264 template <typename T>
1265 Taylor<T>
1266 min(const typename Taylor<T>::value_type& a,
1267  const Base< Taylor<T> >& bb)
1268 {
1269  const Taylor<T>& b = bb.derived();
1270 
1271  if (a <= b.coeff(0))
1272  return Taylor<T>(b.degree(), a);
1273  else
1274  return b;
1275 }
1276 
1277 template <typename T>
1278 Taylor<T>
1279 min(const Base< Taylor<T> >& aa,
1280  const typename Taylor<T>::value_type& b)
1281 {
1282  const Taylor<T>& a = aa.derived();
1283 
1284  if (a.coeff(0) <= b)
1285  return a;
1286  else
1287  return Taylor<T>(a.degree(), b);
1288 }
1289 
1290 template <typename T>
1291 bool
1293  const Base< Taylor<T> >& bb)
1294 {
1295  const Taylor<T>& a = aa.derived();
1296  const Taylor<T>& b = bb.derived();
1297  return a.coeff(0) == b.coeff(0);
1298 }
1299 
1300 template <typename T>
1301 bool
1303  const Base< Taylor<T> >& bb)
1304 {
1305  const Taylor<T>& b = bb.derived();
1306  return a == b.coeff(0);
1307 }
1308 
1309 template <typename T>
1310 bool
1312  const typename Taylor<T>::value_type& b)
1313 {
1314  const Taylor<T>& a = aa.derived();
1315  return a.coeff(0) == b;
1316 }
1317 
1318 template <typename T>
1319 bool
1321  const Base< Taylor<T> >& bb)
1322 {
1323  const Taylor<T>& a = aa.derived();
1324  const Taylor<T>& b = bb.derived();
1325  return a.coeff(0) != b.coeff(0);
1326 }
1327 
1328 template <typename T>
1329 bool
1331  const Base< Taylor<T> >& bb)
1332 {
1333  const Taylor<T>& b = bb.derived();
1334  return a != b.coeff(0);
1335 }
1336 
1337 template <typename T>
1338 bool
1340  const typename Taylor<T>::value_type& b)
1341 {
1342  const Taylor<T>& a = aa.derived();
1343  return a.coeff(0) != b;
1344 }
1345 
1346 template <typename T>
1347 bool
1348 operator<=(const Base< Taylor<T> >& aa,
1349  const Base< Taylor<T> >& bb)
1350 {
1351  const Taylor<T>& a = aa.derived();
1352  const Taylor<T>& b = bb.derived();
1353  return a.coeff(0) <= b.coeff(0);
1354 }
1355 
1356 template <typename T>
1357 bool
1358 operator<=(const typename Taylor<T>::value_type& a,
1359  const Base< Taylor<T> >& bb)
1360 {
1361  const Taylor<T>& b = bb.derived();
1362  return a <= b.coeff(0);
1363 }
1364 
1365 template <typename T>
1366 bool
1367 operator<=(const Base< Taylor<T> >& aa,
1368  const typename Taylor<T>::value_type& b)
1369 {
1370  const Taylor<T>& a = aa.derived();
1371  return a.coeff(0) <= b;
1372 }
1373 
1374 template <typename T>
1375 bool
1377  const Base< Taylor<T> >& bb)
1378 {
1379  const Taylor<T>& a = aa.derived();
1380  const Taylor<T>& b = bb.derived();
1381  return a.coeff(0) >= b.coeff(0);
1382 }
1383 
1384 template <typename T>
1385 bool
1387  const Base< Taylor<T> >& bb)
1388 {
1389  const Taylor<T>& b = bb.derived();
1390  return a >= b.coeff(0);
1391 }
1392 
1393 template <typename T>
1394 bool
1396  const typename Taylor<T>::value_type& b)
1397 {
1398  const Taylor<T>& a = aa.derived();
1399  return a.coeff(0) >= b;
1400 }
1401 
1402 template <typename T>
1403 bool
1404 operator<(const Base< Taylor<T> >& aa,
1405  const Base< Taylor<T> >& bb)
1406 {
1407  const Taylor<T>& a = aa.derived();
1408  const Taylor<T>& b = bb.derived();
1409  return a.coeff(0) < b.coeff(0);
1410 }
1411 
1412 template <typename T>
1413 bool
1414 operator<(const typename Taylor<T>::value_type& a,
1415  const Base< Taylor<T> >& bb)
1416 {
1417  const Taylor<T>& b = bb.derived();
1418  return a < b.coeff(0);
1419 }
1420 
1421 template <typename T>
1422 bool
1423 operator<(const Base< Taylor<T> >& aa,
1424  const typename Taylor<T>::value_type& b)
1425 {
1426  const Taylor<T>& a = aa.derived();
1427  return a.coeff(0) < b;
1428 }
1429 
1430 template <typename T>
1431 bool
1432 operator>(const Base< Taylor<T> >& aa,
1433  const Base< Taylor<T> >& bb)
1434 {
1435  const Taylor<T>& a = aa.derived();
1436  const Taylor<T>& b = bb.derived();
1437  return a.coeff(0) > b.coeff(0);
1438 }
1439 
1440 template <typename T>
1441 bool
1443  const Base< Taylor<T> >& bb)
1444 {
1445  const Taylor<T>& b = bb.derived();
1446  return a > b.coeff(0);
1447 }
1448 
1449 template <typename T>
1450 bool
1451 operator>(const Base< Taylor<T> >& aa,
1452  const typename Taylor<T>::value_type& b)
1453 {
1454  const Taylor<T>& a = aa.derived();
1455  return a.coeff(0) > b;
1456 }
1457 
1458 template <typename T>
1459 bool toBool(const Taylor<T>& x) {
1460  bool is_zero = true;
1461  for (int i=0; i<=x.degree(); i++)
1462  is_zero = is_zero && (x.coeff(i) == 0.0);
1463  return !is_zero;
1464 }
1465 
1466 template <typename T>
1467 inline bool
1468 operator && (const Base< Taylor<T> >& xx1, const Base< Taylor<T> >& xx2)
1469 {
1470  const Taylor<T>& x1 = xx1.derived();
1471  const Taylor<T>& x2 = xx2.derived();
1472  return toBool(x1) && toBool(x2);
1473 }
1474 
1475 template <typename T>
1476 inline bool
1478  const Base< Taylor<T> >& xx2)
1479 {
1480  const Taylor<T>& x2 = xx2.derived();
1481  return a && toBool(x2);
1482 }
1483 
1484 template <typename T>
1485 inline bool
1487  const typename Taylor<T>::value_type& b)
1488 {
1489  const Taylor<T>& x1 = xx1.derived();
1490  return toBool(x1) && b;
1491 }
1492 
1493 template <typename T>
1494 inline bool
1495 operator || (const Base< Taylor<T> >& xx1, const Base< Taylor<T> >& xx2)
1496 {
1497  const Taylor<T>& x1 = xx1.derived();
1498  const Taylor<T>& x2 = xx2.derived();
1499  return toBool(x1) || toBool(x2);
1500 }
1501 
1502 template <typename T>
1503 inline bool
1505  const Base< Taylor<T> >& xx2)
1506 {
1507  const Taylor<T>& x2 = xx2.derived();
1508  return a || toBool(x2);
1509 }
1510 
1511 template <typename T>
1512 inline bool
1514  const typename Taylor<T>::value_type& b)
1515 {
1516  const Taylor<T>& x1 = xx1.derived();
1517  return toBool(x1) || b;
1518 }
1519 
1520 template <typename T>
1521 std::ostream&
1522 operator << (std::ostream& os, const Base< Taylor<T> >& aa)
1523 {
1524  const Taylor<T>& a = aa.derived();
1525  os << "[ ";
1526 
1527  for (int i=0; i<=a.degree(); i++) {
1528  os << a.coeff(i) << " ";
1529  }
1530 
1531  os << "]";
1532  return os;
1533 }
1534 
1535 } // namespace Tay
1536 } // namespace Sacado
Taylor< T > operator+(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
A generic handle class.
ASinExprType< T >::expr_type asin(const Expr< T > &expr)
void resizeCoeffs(int len)
Resize coefficient array to new size.
void resize(int d, bool keep_coeffs)
Resize polynomial to degree d.
void sinhcosh(const Base< Taylor< T > > &a, Taylor< T > &s, Taylor< T > &c)
T value_type
Typename of values.
asin(expr.val())
cosh(expr.val())
int deg_
Degree of polynomial.
void makeOwnCopy()
Make handle have its own copy of rep.
Taylor< T > log(const Base< Taylor< T > > &a)
Taylor< T > & operator-=(const T &x)
Subtraction-assignment operator with constant right-hand-side.
Taylor< T > asinh(const Base< Taylor< T > > &a)
TanhExprType< T >::expr_type tanh(const Expr< T > &expr)
PowExprType< Expr< T1 >, Expr< T2 > >::expr_type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
TanExprType< T >::expr_type tan(const Expr< T > &expr)
static SACADO_INLINE_FUNCTION void zero(T *dest, int sz)
Zero out array dest of length sz.
bool operator>(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > operator+() const
Unary-plus operator.
int len_
Length of allocated polynomial array.
TaylorData & operator=(const TaylorData &x)
Assignment operator.
Taylor< T > & operator*=(const T &x)
Multiplication-assignment operator with constant right-hand-side.
atan(expr.val())
bool operator||(const Base< Taylor< T > > &xx1, const Base< Taylor< T > > &xx2)
Taylor< T > operator/(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > & operator=(const T &val)
Assignment operator with constant right-hand-side.
ACosExprType< T >::expr_type acos(const Expr< T > &expr)
bool operator>=(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Base class for Sacado types to control overload resolution.
Definition: Sacado_Base.hpp:26
Taylor< T > sin(const Base< Taylor< T > > &a)
bool operator&&(const Base< Taylor< T > > &xx1, const Base< Taylor< T > > &xx2)
#define T
Definition: Sacado_rad.hpp:553
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
Taylor< T > cos(const Base< Taylor< T > > &a)
T & fastAccessCoeff(int i)
Returns degree i term without bounds checking.
Taylor< T > operator-() const
Unary-minus operator.
Taylor< T > quad(const typename Taylor< T >::value_type &c0, const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > operator*(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > sinh(const Base< Taylor< T > > &a)
Taylor< T > sqrt(const Base< Taylor< T > > &a)
sqrt(expr.val())
Log10ExprType< T >::expr_type log10(const Expr< T > &expr)
Taylor< T > max(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
sinh(expr.val())
Taylor< T > cbrt(const Base< Taylor< T > > &a)
static SACADO_INLINE_FUNCTION void destroy_and_release(T *m, int sz)
Destroy array elements and release memory.
Taylor< T > fabs(const Base< Taylor< T > > &a)
const T * coeff() const
Returns Taylor coefficient array.
Sacado::Handle< TaylorData > th
ATanExprType< T >::expr_type atan(const Expr< T > &expr)
Taylor< T > abs(const Base< Taylor< T > > &a)
void reserve(int d)
Reserve space for a degree d polynomial.
sin(expr.val())
Taylor< T > atan2(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > exp(const Base< Taylor< T > > &a)
T * coeff_
Taylor polynomial coefficients.
static SACADO_INLINE_FUNCTION void copy(const T *src, T *dest, int sz)
Copy array from src to dest of length sz.
void sincos(const Base< Taylor< T > > &a, Taylor< T > &s, Taylor< T > &c)
Taylor< T > atanh(const Base< Taylor< T > > &a)
log(expr.val())
Taylor polynomial class.
Taylor< T > cosh(const Base< Taylor< T > > &a)
bool operator==(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
acos(expr.val())
const derived_type & derived() const
Definition: Sacado_Base.hpp:28
Taylor< T > min(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
Taylor< T > operator-(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
exp(expr.val())
Taylor< T > & operator/=(const T &x)
Division-assignment operator with constant right-hand-side.
Taylor< T > acosh(const Base< Taylor< T > > &a)
Taylor< T > & operator+=(const T &x)
Addition-assignment operator with constant right-hand-side.
int degree() const
Returns degree of polynomial.
bool operator!=(const Base< Taylor< T > > &a, const Base< Taylor< T > > &b)
cos(expr.val())
static SACADO_INLINE_FUNCTION T * get_and_fill(int sz)
Get memory for new array of length sz and fill with zeros.
Taylor()
Default constructor.
bool toBool(const Taylor< T > &x)