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