Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2941-gbd3c83c843 2025-03-25 23:50: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
21#include <deal.II/base/types.h>
22
23#include <cmath>
24#include <complex>
25#include <cstddef>
26#include <type_traits>
27
28// Forward-declare the automatic differentiation types so we can add prototypes
29// for our own wrappers.
30#ifdef DEAL_II_WITH_ADOLC
31class adouble;
32namespace adtl
33{
34 class adouble;
35}
36#endif
37
39
40namespace internal
41{
58 template <typename Number>
60 {
64 constexpr static unsigned int max_width = 1;
65 };
66
73 template <>
75 {
79 constexpr static unsigned int max_width =
80#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
81 8;
82#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
83 4;
84#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
85 2;
86#else
87 1;
88#endif
89 };
90
97 template <>
99 {
103 constexpr static unsigned int max_width =
104#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
105 4;
106#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
107 16;
108#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
109 8;
110#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
111 4;
112#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
113 4;
114#else
115 1;
116#endif
117 };
118
119
120} // namespace internal
121
122// forward declarations to support abs or sqrt operations on VectorizedArray
123#ifndef DOXYGEN
124template <typename Number,
125 std::size_t width =
127class VectorizedArray;
128template <typename T>
129struct EnableIfScalar;
130#endif
131
132#ifdef DEAL_II_WITH_ADOLC
133# ifndef DOXYGEN
134// Prototype some inline functions present in adolc_math.h for use in
135// NumberTraits.
136//
137// ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
138// though, we don't want to include ADOL-C headers in this header since
139// numbers.h is in everything. To get around this: use C++ rules which permit
140// the use of forward-declared classes in function prototypes to declare some
141// functions which are defined in adolc_math.h. This permits us to write "using
142// ::abs;" in NumberTraits which will allow us to select the correct
143// overload (the one in ::) when instantiating NumberTraits for ADOL-C
144// types.
145
146adouble
147abs(const adouble &x);
148
149adtl::adouble
150abs(const adtl::adouble &x);
151# endif
152#endif
153
155
156namespace std
157{
158 template <typename Number, std::size_t width>
159 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
160 sqrt(const ::VectorizedArray<Number, width> &);
161 template <typename Number, std::size_t width>
162 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
163 abs(const ::VectorizedArray<Number, width> &);
164 template <typename Number, std::size_t width>
165 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
166 max(const ::VectorizedArray<Number, width> &,
167 const ::VectorizedArray<Number, width> &);
168 template <typename Number, std::size_t width>
169 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
170 min(const ::VectorizedArray<Number, width> &,
171 const ::VectorizedArray<Number, width> &);
172 template <typename Number, size_t width>
174 pow(const ::VectorizedArray<Number, width> &, const Number p);
175 template <typename Number, size_t width>
177 sin(const ::VectorizedArray<Number, width> &);
178 template <typename Number, size_t width>
180 cos(const ::VectorizedArray<Number, width> &);
181 template <typename Number, size_t width>
183 tan(const ::VectorizedArray<Number, width> &);
184 template <typename Number, size_t width>
186 exp(const ::VectorizedArray<Number, width> &);
187 template <typename Number, size_t width>
189 log(const ::VectorizedArray<Number, width> &);
190} // namespace std
191
193
209namespace numbers
210{
214 constexpr double E = 2.7182818284590452354;
215
219 constexpr double LOG2E = 1.4426950408889634074;
220
224 constexpr double LOG10E = 0.43429448190325182765;
225
229 constexpr double LN2 = 0.69314718055994530942;
230
234 constexpr double LN10 = 2.30258509299404568402;
235
239 constexpr double PI = 3.14159265358979323846;
240
244 constexpr double PI_2 = 1.57079632679489661923;
245
249 constexpr double PI_4 = 0.78539816339744830962;
250
254 constexpr double SQRT2 = 1.41421356237309504880;
255
259 constexpr double SQRT1_2 = 0.70710678118654752440;
260
270 bool
271 is_finite(const double x);
272
277 bool
278 is_finite(const std::complex<double> &x);
279
284 bool
285 is_finite(const std::complex<float> &x);
286
295 bool
296 is_finite(const std::complex<long double> &x);
297
308 template <typename Number1, typename Number2>
309 constexpr DEAL_II_HOST_DEVICE bool
310 values_are_equal(const Number1 &value_1, const Number2 &value_2);
311
322 template <typename Number1, typename Number2>
323 constexpr bool
324 values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
325
333 template <typename Number>
334 constexpr DEAL_II_HOST_DEVICE bool
335 value_is_zero(const Number &value);
336
347 template <typename Number1, typename Number2>
348 bool
349 value_is_less_than(const Number1 &value_1, const Number2 &value_2);
350
361 template <typename Number1, typename Number2>
362 bool
363 value_is_less_than_or_equal_to(const Number1 &value_1,
364 const Number2 &value_2);
365
366
367
378 template <typename Number1, typename Number2>
379 bool
380 value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
381
392 template <typename Number1, typename Number2>
393 bool
394 value_is_greater_than_or_equal_to(const Number1 &value_1,
395 const Number2 &value_2);
396
405 template <typename number>
407 {
413 static constexpr bool is_complex = false;
414
421 using real_type = number;
422
426 using double_type = double;
427
435 static constexpr DEAL_II_HOST_DEVICE const number &
436 conjugate(const number &x);
437
446 static constexpr DEAL_II_HOST_DEVICE real_type
447 abs_square(const number &x);
448
452 static real_type
453 abs(const number &x);
454 };
455
456
461 template <typename number>
462 struct NumberTraits<std::complex<number>>
463 {
469 static constexpr bool is_complex = true;
470
477 using real_type = number;
478
482 using double_type = std::complex<double>;
483
487 static constexpr std::complex<number>
488 conjugate(const std::complex<number> &x);
489
496 static constexpr real_type
497 abs_square(const std::complex<number> &x);
498
499
503 static real_type
504 abs(const std::complex<number> &x);
505 };
506
507 // --------------- inline and template functions ---------------- //
508
509 inline bool
510 is_nan(const double x)
511 {
512 return std::isnan(x);
513 }
514
515
516
517 inline bool
518 is_finite(const double x)
519 {
520 return std::isfinite(x);
521 }
522
523
524
525 inline bool
526 is_finite(const std::complex<double> &x)
527 {
528 // Check complex numbers for infinity
529 // by testing real and imaginary part
530 return (is_finite(x.real()) && is_finite(x.imag()));
531 }
532
533
534
535 inline bool
536 is_finite(const std::complex<float> &x)
537 {
538 // Check complex numbers for infinity
539 // by testing real and imaginary part
540 return (is_finite(x.real()) && is_finite(x.imag()));
541 }
542
543
544
545 inline bool
546 is_finite(const std::complex<long double> &x)
547 {
548 // Same for std::complex<long double>
549 return (is_finite(x.real()) && is_finite(x.imag()));
550 }
551
552
553 template <typename number>
554 constexpr DEAL_II_HOST_DEVICE const number &
556 {
557 return x;
558 }
559
560
561
562 template <typename number>
565 {
566 return x * x;
567 }
568
569
570
571 template <typename number>
574 {
575 // Make things work with AD types
576 using std::abs;
577#ifdef DEAL_II_WITH_ADOLC
578 // This one is a little tricky - we have our own abs function in ::,
579 // prototyped with forward-declared types in this file, but it only exists
580 // if we have ADOL-C: hence we only add this using statement in that
581 // situation
582 using ::abs;
583#endif
584 return abs(x);
585 }
586
587
588
589 template <typename number>
590 constexpr std::complex<number>
591 NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
592 {
593 return std::conj(x);
594 }
595
596
597
598 template <typename number>
599 typename NumberTraits<std::complex<number>>::real_type
600 NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
601 {
602 // Make things work with AD types
603 using std::abs;
604#ifdef DEAL_II_WITH_ADOLC
605 // Same comment as the non-complex case holds here
606 using ::abs;
607#endif
608 return abs(x);
609 }
610
611
612
613 template <typename number>
614 constexpr typename NumberTraits<std::complex<number>>::real_type
615 NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
616 {
617 return std::norm(x);
618 }
619
620} // namespace numbers
621
622
623// Forward declarations
625{
626 namespace AD
627 {
628 namespace internal
629 {
630 // Defined in differentiation/ad/ad_number_traits.h
631 template <typename T>
633 } // namespace internal
634
635 // Defined in differentiation/ad/ad_number_traits.h
636 template <typename NumberType>
638 } // namespace AD
639} // namespace Differentiation
640
641
642namespace internal
643{
648 template <typename From, typename To>
650 {
651 // Source: https://stackoverflow.com/a/16944130
652 private:
653 template <typename T>
654 static void f(T);
655
656 template <typename F, typename T>
657 static constexpr auto
658 test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
659 {
660 return true;
661 }
662
663 template <typename F, typename T>
664 static constexpr auto
665 test(...) -> bool
666 {
667 return false;
668 }
669
670 public:
671 static const bool value = test<From, To>(0);
672 };
673
674 /*
675 * The structs below are needed to convert between some special number types.
676 * Also see tensor.h for another specialization.
677 */
678 template <typename T>
680 {
681 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
682 value(const T &t)
683 {
684 return t;
685 }
686
687 // Below are generic functions that allows an overload for any
688 // type U that is transformable to type T. This is particularly
689 // useful when needing to cast exotic number types
690 // (e.g. auto-differentiable or symbolic numbers) to a floating
691 // point one, such as might happen when converting between tensor
692 // types.
693
694 // Type T is constructible from F.
695 template <typename F>
696 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
697 value(const F &f,
698 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
699 std::is_constructible_v<T, F>> * = nullptr)
700 {
701 return T(f);
702 }
703
704 // Type T is explicitly convertible (but not constructible) from F.
705 template <typename F>
706 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
707 value(const F &f,
708 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
709 !std::is_constructible_v<T, F> &&
711 nullptr)
712 {
713 return static_cast<T>(f);
714 }
715
716 // Sacado doesn't provide any conversion operators, so we have
717 // to extract the value and perform further conversions from there.
718 // To be safe, we extend this to other possible AD numbers that
719 // might fall into the same category.
720 template <typename F>
721 static T
723 const F &f,
724 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
725 !std::is_constructible_v<T, F> &&
728 {
730 }
731 };
732
733 template <typename T>
734 struct NumberType<std::complex<T>>
735 {
736 static constexpr const std::complex<T> &
737 value(const std::complex<T> &t)
738 {
739 return t;
740 }
741
742 static constexpr std::complex<T>
743 value(const T &t)
744 {
745 return std::complex<T>(t);
746 }
747
748 // Facilitate cast from complex<double> to complex<float>
749 template <typename U>
750 static constexpr std::complex<T>
751 value(const std::complex<U> &t)
752 {
753 return std::complex<T>(NumberType<T>::value(t.real()),
754 NumberType<T>::value(t.imag()));
755 }
756 };
757
758} // namespace internal
759
760namespace numbers
761{
762#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
763
774 // Defined in differentiation/ad/adolc_number_types.cc
775 bool
776 values_are_equal(const adouble &value_1, const adouble &value_2);
777
778
789 template <typename Number>
790 bool
791 values_are_equal(const adouble &value_1, const Number &value_2)
792 {
793 // Use the specialized definition for two ADOL-C taped types
794 return values_are_equal(
796 }
797
798
809 template <typename Number>
810 bool
811 values_are_equal(const Number &value_1, const adouble &value_2)
812 {
813 // Use the above definition
814 return values_are_equal(value_2, value_1);
815 }
816
828 // Defined in differentiation/ad/adolc_number_types.cc
829 bool
830 value_is_less_than(const adouble &value_1, const adouble &value_2);
831
832
844 template <typename Number>
845 bool
846 value_is_less_than(const adouble &value_1, const Number &value_2)
847 {
848 // Use the specialized definition for two ADOL-C taped types
849 return value_is_less_than(
851 }
852
853
865 template <typename Number>
866 bool
867 value_is_less_than(const Number &value_1, const adouble &value_2)
868 {
869 // Use the specialized definition for two ADOL-C taped types
870 return value_is_less_than(
872 }
873
874#endif
875
876
877 template <typename Number1, typename Number2>
878 constexpr DEAL_II_HOST_DEVICE bool
879 values_are_equal(const Number1 &value_1, const Number2 &value_2)
880 {
881 return (value_1 == ::internal::NumberType<Number1>::value(value_2));
882 }
883
884
885 template <typename Number1, typename Number2>
886 inline constexpr bool
887 values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
888 {
889 return !(values_are_equal(value_1, value_2));
890 }
891
892
893 template <typename Number>
894 constexpr DEAL_II_HOST_DEVICE bool
895 value_is_zero(const Number &value)
896 {
897 return values_are_equal(value, 0.0);
898 }
899
900
901 template <typename Number1, typename Number2>
902 inline bool
903 value_is_less_than(const Number1 &value_1, const Number2 &value_2)
904 {
905 return (value_1 < ::internal::NumberType<Number1>::value(value_2));
906 }
907
908
909 template <typename Number1, typename Number2>
910 inline bool
911 value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
912 {
913 return (value_is_less_than(value_1, value_2) ||
914 values_are_equal(value_1, value_2));
915 }
916
917
918 template <typename Number1, typename Number2>
919 bool
920 value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
921 {
922 return !(value_is_less_than_or_equal_to(value_1, value_2));
923 }
924
925
926 template <typename Number1, typename Number2>
927 inline bool
928 value_is_greater_than_or_equal_to(const Number1 &value_1,
929 const Number2 &value_2)
930 {
931 return !(value_is_less_than(value_1, value_2));
932 }
933} // namespace numbers
934
936
937#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_HOST_DEVICE
Definition config.h:169
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition config.h:170
Definition numbers.h:33
constexpr double LOG10E
Definition numbers.h:224
constexpr double PI_2
Definition numbers.h:244
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:911
constexpr double E
Definition numbers.h:214
constexpr double PI
Definition numbers.h:239
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:920
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:887
constexpr double SQRT2
Definition numbers.h:254
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:895
constexpr double SQRT1_2
Definition numbers.h:259
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:879
constexpr double PI_4
Definition numbers.h:249
constexpr double LN10
Definition numbers.h:234
constexpr double LN2
Definition numbers.h:229
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:903
bool is_finite(const double x)
Definition numbers.h:518
constexpr double LOG2E
Definition numbers.h:219
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:928
bool is_nan(const double x)
Definition numbers.h:510
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:751
static constexpr std::complex< T > value(const T &t)
Definition numbers.h:743
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition numbers.h:737
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:697
static constexpr const T & value(const T &t)
Definition numbers.h:682
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:707
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:722
static constexpr unsigned int max_width
Definition numbers.h:64
static constexpr auto test(...) -> bool
Definition numbers.h:665
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition numbers.h:658
static constexpr const number & conjugate(const number &x)
Definition numbers.h:555
static constexpr bool is_complex
Definition numbers.h:413
static real_type abs(const number &x)
Definition numbers.h:573
static constexpr real_type abs_square(const number &x)
Definition numbers.h:564