Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
sacado_math.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2017 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 
25 # include <deal.II/differentiation/ad/ad_number_traits.h>
26 # include <deal.II/differentiation/ad/sacado_number_types.h>
27 
28 
29 # ifndef DOXYGEN
30 
34 namespace std
35 {
39  template <
40  typename ADNumberType,
41  typename = typename std::enable_if<
44  ADNumberType>::value>::type>
45  inline ADNumberType
46  erf(ADNumberType x)
47  {
48  // Reference:
49  // Handbook of Mathematical Functions: with Formulas, Graphs, and
50  // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
51  // Mathematics. 1972.
52  //
53  // Current implementation source:
54  // https://www.johndcook.com/blog/cpp_erf/
55  // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
56  //
57  // Note: This implementation has a reported maximum round-off error
58  // of 1.5e-7.
59 
60  // Constants
61  const double a1 = 0.254829592;
62  const double a2 = -0.284496736;
63  const double a3 = 1.421413741;
64  const double a4 = -1.453152027;
65  const double a5 = 1.061405429;
66  const double p = 0.3275911;
67 
68  // Save the sign of x
69  const bool neg_val =
70  (x < ::internal::NumberType<ADNumberType>::value(0.0) ? true :
71  false);
72  x = std::fabs(x);
73 
74  // Abramowitz1972a equation 7.1.26
75  const ADNumberType t = 1.0 / (1.0 + p * x);
76  const ADNumberType y =
77  1.0 -
78  (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
79 
80  if (!neg_val)
81  return y;
82  else
83  return -y;
84  }
85 
86 
90  template <
91  typename ADNumberType,
92  typename = typename std::enable_if<
94  inline ADNumberType
95  erfc(const ADNumberType &x)
96  {
97  return 1.0 - std::erf(x);
98  }
99 
100 } // namespace std
101 
102 # endif // DOXYGEN
103 
104 #endif // DEAL_II_TRILINOS_WITH_SACADO
105 
106 #endif
STL namespace.
Expression erfc(const Expression &x)
Expression erf(const Expression &x)