Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
sacado_math.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_differentiation_ad_sacado_math_h
17#define dealii_differentiation_ad_sacado_math_h
18
19#include <deal.II/base/config.h>
20
21#ifdef DEAL_II_TRILINOS_WITH_SACADO
22
23# include <deal.II/base/numbers.h>
24
27
28
29# ifndef DOXYGEN
30
35template <
36 typename ADNumberType,
37 typename = std::enable_if_t<
40inline ADNumberType
41erf(ADNumberType x)
42{
43 // Reference:
44 // Handbook of Mathematical Functions: with Formulas, Graphs, and
45 // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
46 // Mathematics. 1972.
47 //
48 // Current implementation source:
49 // https://www.johndcook.com/blog/cpp_erf/
50 // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
51 //
52 // Note: This implementation has a reported maximum round-off error
53 // of 1.5e-7.
54
55 // Constants
56 const double a1 = 0.254829592;
57 const double a2 = -0.284496736;
58 const double a3 = 1.421413741;
59 const double a4 = -1.453152027;
60 const double a5 = 1.061405429;
61 const double p = 0.3275911;
62
63 // Save the sign of x
64 const bool neg_val =
65 (x < ::internal::NumberType<ADNumberType>::value(0.0) ? true : false);
66 x = std::fabs(x);
67
68 // Abramowitz1972a equation 7.1.26
69 const ADNumberType t = 1.0 / (1.0 + p * x);
70 const ADNumberType y =
71 1.0 -
72 (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
73
74 if (!neg_val)
75 return y;
76 else
77 return -y;
78}
79
80
84template <typename ADNumberType,
85 typename = std::enable_if_t<
87inline ADNumberType
88erfc(const ADNumberType &x)
89{
90 return 1.0 - std::erf(x);
91}
92
94
95# endif // DOXYGEN
96
97#endif // DEAL_II_TRILINOS_WITH_SACADO
98
99#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Expression erfc(const Expression &x)
Expression erf(const Expression &x)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)