Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
numbers.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2006 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_numbers_h
16#define dealii_numbers_h
17
18
19#include <deal.II/base/config.h>
20
22#include <deal.II/base/types.h>
23
24#include <cmath>
25#include <complex>
26#include <cstddef>
27#include <type_traits>
28
29// Forward-declare the automatic differentiation types so we can add prototypes
30// for our own wrappers.
31#ifdef DEAL_II_WITH_ADOLC
32class adouble;
33namespace adtl
34{
35 class adouble;
36}
37#endif
38
40
41namespace internal
42{
59 template <typename Number>
61 {
65 constexpr static unsigned int max_width = 1;
66 };
67
74 template <>
76 {
80 constexpr static unsigned int max_width =
81#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
82 8;
83#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
84 4;
85#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
86 2;
87#else
88 1;
89#endif
90 };
91
98 template <>
100 {
104 constexpr static unsigned int max_width =
105#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
106 4;
107#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
108 16;
109#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
110 8;
111#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
112 4;
113#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
114 4;
115#else
116 1;
117#endif
118 };
119
120
121} // namespace internal
122
123// forward declarations to support abs or sqrt operations on VectorizedArray
124#ifndef DOXYGEN
125template <typename Number,
126 std::size_t width =
128class VectorizedArray;
129template <typename T>
130struct EnableIfScalar;
131#endif
132
133#ifdef DEAL_II_WITH_ADOLC
134# ifndef DOXYGEN
135// Prototype some inline functions present in adolc_math.h for use in
136// NumberTraits.
137//
138// ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
139// though, we don't want to include ADOL-C headers in this header since
140// numbers.h is in everything. To get around this: use C++ rules which permit
141// the use of forward-declared classes in function prototypes to declare some
142// functions which are defined in adolc_math.h. This permits us to write "using
143// ::abs;" in NumberTraits which will allow us to select the correct
144// overload (the one in ::) when instantiating NumberTraits for ADOL-C
145// types.
146
147adouble
148abs(const adouble &x);
149
150adtl::adouble
151abs(const adtl::adouble &x);
152# endif
153#endif
154
156
157namespace std
158{
159 template <typename Number, std::size_t width>
160 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
161 sqrt(const ::VectorizedArray<Number, width> &);
162 template <typename Number, std::size_t width>
163 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
164 abs(const ::VectorizedArray<Number, width> &);
165 template <typename Number, std::size_t width>
166 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
167 max(const ::VectorizedArray<Number, width> &,
168 const ::VectorizedArray<Number, width> &);
169 template <typename Number, std::size_t width>
170 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
171 min(const ::VectorizedArray<Number, width> &,
172 const ::VectorizedArray<Number, width> &);
173 template <typename Number, size_t width>
175 pow(const ::VectorizedArray<Number, width> &, const Number p);
176 template <typename Number, size_t width>
178 sin(const ::VectorizedArray<Number, width> &);
179 template <typename Number, size_t width>
181 cos(const ::VectorizedArray<Number, width> &);
182 template <typename Number, size_t width>
184 tan(const ::VectorizedArray<Number, width> &);
185 template <typename Number, size_t width>
187 exp(const ::VectorizedArray<Number, width> &);
188 template <typename Number, size_t width>
190 log(const ::VectorizedArray<Number, width> &);
191} // namespace std
192
194
210namespace numbers
211{
215 constexpr double E = 2.7182818284590452354;
216
220 constexpr double LOG2E = 1.4426950408889634074;
221
225 constexpr double LOG10E = 0.43429448190325182765;
226
230 constexpr double LN2 = 0.69314718055994530942;
231
235 constexpr double LN10 = 2.30258509299404568402;
236
240 constexpr double PI = 3.14159265358979323846;
241
245 constexpr double PI_2 = 1.57079632679489661923;
246
250 constexpr double PI_4 = 0.78539816339744830962;
251
255 constexpr double SQRT2 = 1.41421356237309504880;
256
260 constexpr double SQRT1_2 = 0.70710678118654752440;
261
271 bool
272 is_finite(const double x);
273
278 bool
279 is_finite(const std::complex<double> &x);
280
285 bool
286 is_finite(const std::complex<float> &x);
287
296 bool
297 is_finite(const std::complex<long double> &x);
298
309 template <typename Number1, typename Number2>
310 constexpr DEAL_II_HOST_DEVICE bool
311 values_are_equal(const Number1 &value_1, const Number2 &value_2);
312
323 template <typename Number1, typename Number2>
324 constexpr bool
325 values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
326
334 template <typename Number>
335 constexpr DEAL_II_HOST_DEVICE bool
336 value_is_zero(const Number &value);
337
348 template <typename Number1, typename Number2>
349 bool
350 value_is_less_than(const Number1 &value_1, const Number2 &value_2);
351
362 template <typename Number1, typename Number2>
363 bool
364 value_is_less_than_or_equal_to(const Number1 &value_1,
365 const Number2 &value_2);
366
367
368
379 template <typename Number1, typename Number2>
380 bool
381 value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
382
393 template <typename Number1, typename Number2>
394 bool
395 value_is_greater_than_or_equal_to(const Number1 &value_1,
396 const Number2 &value_2);
397
406 template <typename number>
408 {
414 static constexpr bool is_complex = false;
415
422 using real_type = number;
423
427 using double_type = double;
428
436 static constexpr DEAL_II_HOST_DEVICE const number &
437 conjugate(const number &x);
438
447 static constexpr DEAL_II_HOST_DEVICE real_type
448 abs_square(const number &x);
449
453 static real_type
454 abs(const number &x);
455 };
456
457
462 template <typename number>
463 struct NumberTraits<std::complex<number>>
464 {
470 static constexpr bool is_complex = true;
471
478 using real_type = number;
479
483 using double_type = std::complex<double>;
484
488 static constexpr std::complex<number>
489 conjugate(const std::complex<number> &x);
490
497 static constexpr real_type
498 abs_square(const std::complex<number> &x);
499
500
504 static real_type
505 abs(const std::complex<number> &x);
506 };
507
508 // --------------- inline and template functions ---------------- //
509
510 inline bool
511 is_nan(const double x)
512 {
513 return std::isnan(x);
514 }
515
516
517
518 inline bool
519 is_finite(const double x)
520 {
521 return std::isfinite(x);
522 }
523
524
525
526 inline bool
527 is_finite(const std::complex<double> &x)
528 {
529 // Check complex numbers for infinity
530 // by testing real and imaginary part
531 return (is_finite(x.real()) && is_finite(x.imag()));
532 }
533
534
535
536 inline bool
537 is_finite(const std::complex<float> &x)
538 {
539 // Check complex numbers for infinity
540 // by testing real and imaginary part
541 return (is_finite(x.real()) && is_finite(x.imag()));
542 }
543
544
545
546 inline bool
547 is_finite(const std::complex<long double> &x)
548 {
549 // Same for std::complex<long double>
550 return (is_finite(x.real()) && is_finite(x.imag()));
551 }
552
553
554 template <typename number>
555 constexpr DEAL_II_HOST_DEVICE const number &
557 {
558 return x;
559 }
560
561
562
563 template <typename number>
566 {
567 return x * x;
568 }
569
570
571
572 template <typename number>
575 {
576 // Make things work with AD types
577 using std::abs;
578#ifdef DEAL_II_WITH_ADOLC
579 // This one is a little tricky - we have our own abs function in ::,
580 // prototyped with forward-declared types in this file, but it only exists
581 // if we have ADOL-C: hence we only add this using statement in that
582 // situation
583 using ::abs;
584#endif
585 return abs(x);
586 }
587
588
589
590 template <typename number>
591 constexpr std::complex<number>
592 NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
593 {
594 return std::conj(x);
595 }
596
597
598
599 template <typename number>
600 typename NumberTraits<std::complex<number>>::real_type
601 NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
602 {
603 // Make things work with AD types
604 using std::abs;
605#ifdef DEAL_II_WITH_ADOLC
606 // Same comment as the non-complex case holds here
607 using ::abs;
608#endif
609 return abs(x);
610 }
611
612
613
614 template <typename number>
615 constexpr typename NumberTraits<std::complex<number>>::real_type
616 NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
617 {
618 return std::norm(x);
619 }
620
621} // namespace numbers
622
623
624// Forward declarations
626{
627 namespace AD
628 {
629 namespace internal
630 {
631 // Defined in differentiation/ad/ad_number_traits.h
632 template <typename T>
634 } // namespace internal
635
636 // Defined in differentiation/ad/ad_number_traits.h
637 template <typename NumberType>
639 } // namespace AD
640} // namespace Differentiation
641
642
643namespace internal
644{
649 template <typename From, typename To>
651 {
652 // Source: https://stackoverflow.com/a/16944130
653 private:
654 template <typename T>
655 static void f(T);
656
657 template <typename F, typename T>
658 static constexpr auto
659 test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
660 {
661 return true;
662 }
663
664 template <typename F, typename T>
665 static constexpr auto
666 test(...) -> bool
667 {
668 return false;
669 }
670
671 public:
672 static const bool value = test<From, To>(0);
673 };
674
675 /*
676 * The structs below are needed to convert between some special number types.
677 * Also see tensor.h for another specialization.
678 */
679 template <typename T>
681 {
682 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
683 value(const T &t)
684 {
685 return t;
686 }
687
688 // Below are generic functions that allows an overload for any
689 // type U that is transformable to type T. This is particularly
690 // useful when needing to cast exotic number types
691 // (e.g. auto-differentiable or symbolic numbers) to a floating
692 // point one, such as might happen when converting between tensor
693 // types.
694
695 // Type T is constructible from F.
696 template <typename F>
697 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
698 value(const F &f,
699 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
700 std::is_constructible_v<T, F>> * = nullptr)
701 {
702 return T(f);
703 }
704
705 // Type T is explicitly convertible (but not constructible) from F.
706 template <typename F>
707 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
708 value(const F &f,
709 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
710 !std::is_constructible_v<T, F> &&
712 nullptr)
713 {
714 return static_cast<T>(f);
715 }
716
717 // Sacado doesn't provide any conversion operators, so we have
718 // to extract the value and perform further conversions from there.
719 // To be safe, we extend this to other possible AD numbers that
720 // might fall into the same category.
721 template <typename F>
722 static T
724 const F &f,
725 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
726 !std::is_constructible_v<T, F> &&
729 {
731 }
732 };
733
734 template <typename T>
735 struct NumberType<std::complex<T>>
736 {
737 static constexpr const std::complex<T> &
738 value(const std::complex<T> &t)
739 {
740 return t;
741 }
742
743 static constexpr std::complex<T>
744 value(const T &t)
745 {
746 return std::complex<T>(t);
747 }
748
749 // Facilitate cast from complex<double> to complex<float>
750 template <typename U>
751 static constexpr std::complex<T>
752 value(const std::complex<U> &t)
753 {
754 return std::complex<T>(NumberType<T>::value(t.real()),
755 NumberType<T>::value(t.imag()));
756 }
757 };
758
759} // namespace internal
760
761namespace numbers
762{
763#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
764
775 // Defined in differentiation/ad/adolc_number_types.cc
776 bool
777 values_are_equal(const adouble &value_1, const adouble &value_2);
778
779
790 template <typename Number>
791 bool
792 values_are_equal(const adouble &value_1, const Number &value_2)
793 {
794 // Use the specialized definition for two ADOL-C taped types
795 return values_are_equal(
797 }
798
799
810 template <typename Number>
811 bool
812 values_are_equal(const Number &value_1, const adouble &value_2)
813 {
814 // Use the above definition
815 return values_are_equal(value_2, value_1);
816 }
817
829 // Defined in differentiation/ad/adolc_number_types.cc
830 bool
831 value_is_less_than(const adouble &value_1, const adouble &value_2);
832
833
845 template <typename Number>
846 bool
847 value_is_less_than(const adouble &value_1, const Number &value_2)
848 {
849 // Use the specialized definition for two ADOL-C taped types
850 return value_is_less_than(
852 }
853
854
866 template <typename Number>
867 bool
868 value_is_less_than(const Number &value_1, const adouble &value_2)
869 {
870 // Use the specialized definition for two ADOL-C taped types
871 return value_is_less_than(
873 }
874
875#endif
876
877
878 template <typename Number1, typename Number2>
879 constexpr DEAL_II_HOST_DEVICE bool
880 values_are_equal(const Number1 &value_1, const Number2 &value_2)
881 {
882 return (value_1 == ::internal::NumberType<Number1>::value(value_2));
883 }
884
885
886 template <typename Number1, typename Number2>
887 inline constexpr bool
888 values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
889 {
890 return !(values_are_equal(value_1, value_2));
891 }
892
893
894 template <typename Number>
895 constexpr DEAL_II_HOST_DEVICE bool
896 value_is_zero(const Number &value)
897 {
898 return values_are_equal(value, 0.0);
899 }
900
901
902 template <typename Number1, typename Number2>
903 inline bool
904 value_is_less_than(const Number1 &value_1, const Number2 &value_2)
905 {
906 return (value_1 < ::internal::NumberType<Number1>::value(value_2));
907 }
908
909
910 template <typename Number1, typename Number2>
911 inline bool
912 value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
913 {
914 return (value_is_less_than(value_1, value_2) ||
915 values_are_equal(value_1, value_2));
916 }
917
918
919 template <typename Number1, typename Number2>
920 bool
921 value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
922 {
923 return !(value_is_less_than_or_equal_to(value_1, value_2));
924 }
925
926
927 template <typename Number1, typename Number2>
928 inline bool
929 value_is_greater_than_or_equal_to(const Number1 &value_1,
930 const Number2 &value_2)
931 {
932 return !(value_is_less_than(value_1, value_2));
933 }
934} // namespace numbers
935
937
938#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition config.h:167
Definition numbers.h:34
constexpr double LOG10E
Definition numbers.h:225
constexpr double PI_2
Definition numbers.h:245
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:912
constexpr double E
Definition numbers.h:215
constexpr double PI
Definition numbers.h:240
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:921
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:888
constexpr double SQRT2
Definition numbers.h:255
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:896
constexpr double SQRT1_2
Definition numbers.h:260
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:880
constexpr double PI_4
Definition numbers.h:250
constexpr double LN10
Definition numbers.h:235
constexpr double LN2
Definition numbers.h:230
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:904
bool is_finite(const double x)
Definition numbers.h:519
constexpr double LOG2E
Definition numbers.h:220
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:929
bool is_nan(const double x)
Definition numbers.h:511
STL namespace.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition numbers.h:752
static constexpr std::complex< T > value(const T &t)
Definition numbers.h:744
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition numbers.h:738
static constexpr T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&std::is_constructible_v< T, F > > *=nullptr)
Definition numbers.h:698
static constexpr const T & value(const T &t)
Definition numbers.h:683
static constexpr T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&is_explicitly_convertible< const F, T >::value > *=nullptr)
Definition numbers.h:708
static T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition numbers.h:723
static constexpr unsigned int max_width
Definition numbers.h:65
static constexpr auto test(...) -> bool
Definition numbers.h:666
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition numbers.h:659
static constexpr const number & conjugate(const number &x)
Definition numbers.h:556
static constexpr bool is_complex
Definition numbers.h:414
static real_type abs(const number &x)
Definition numbers.h:574
static constexpr real_type abs_square(const number &x)
Definition numbers.h:565