1/* Implementation of the ERFC_SCALED intrinsic.
2   Copyright (C) 2008-2020 Free Software Foundation, Inc.
3
4This file is part of the GNU Fortran runtime library (libgfortran).
5
6Libgfortran is free software; you can redistribute it and/or
7modify it under the terms of the GNU General Public
8License as published by the Free Software Foundation; either
9version 3 of the License, or (at your option) any later version.
10
11Libgfortran is distributed in the hope that it will be useful,
12but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14GNU General Public License for more details.
15
16Under Section 7 of GPL version 3, you are granted additional
17permissions described in the GCC Runtime Library Exception, version
183.1, as published by the Free Software Foundation.
19
20You should have received a copy of the GNU General Public License and
21a copy of the GCC Runtime Library Exception along with this program;
22see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23<http://www.gnu.org/licenses/>.  */
24
25#include "libgfortran.h"
26
27/* This implementation of ERFC_SCALED is based on the netlib algorithm
28   available at http://www.netlib.org/specfun/erf  */
29
30#ifdef HAVE_GFC_REAL_4
31#undef KIND
32#define KIND 4
33#include "erfc_scaled_inc.c"
34#endif
35
36#ifdef HAVE_GFC_REAL_8
37#undef KIND
38#define KIND 8
39#include "erfc_scaled_inc.c"
40#endif
41
42#ifdef HAVE_GFC_REAL_10
43#undef KIND
44#define KIND 10
45#include "erfc_scaled_inc.c"
46#endif
47
48#ifdef HAVE_GFC_REAL_16
49
50/* For quadruple-precision, netlib's implementation is
51   not accurate enough.  We provide another one.  */
52
53#ifdef GFC_REAL_16_IS_FLOAT128
54
55# define _THRESH -106.566990228185312813205074546585730Q
56# define _M_2_SQRTPI M_2_SQRTPIq
57# define _INF __builtin_infq()
58# define _ERFC(x) erfcq(x)
59# define _EXP(x) expq(x)
60
61#else
62
63# define _THRESH -106.566990228185312813205074546585730L
64# ifndef M_2_SQRTPIl
65#  define M_2_SQRTPIl 1.128379167095512573896158903121545172L
66# endif
67# define _M_2_SQRTPI M_2_SQRTPIl
68# define _INF __builtin_infl()
69# ifdef HAVE_ERFCL
70#  define _ERFC(x) erfcl(x)
71# endif
72# ifdef HAVE_EXPL
73#  define _EXP(x) expl(x)
74# endif
75
76#endif
77
78#if defined(_ERFC) && defined(_EXP)
79
80extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
81export_proto(erfc_scaled_r16);
82
83GFC_REAL_16
84erfc_scaled_r16 (GFC_REAL_16 x)
85{
86  if (x < _THRESH)
87    {
88      return _INF;
89    }
90  if (x < 12)
91    {
92      /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
93	 This is not perfect, but much better than netlib.  */
94      return _ERFC(x) * _EXP(x * x);
95    }
96  else
97    {
98      /* Calculate ERFC_SCALED(x) using a power series in 1/x:
99	 ERFC_SCALED(x) = 1 / (x * sqrt(pi))
100			 * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
101					      / (2 * x**2)**n)
102       */
103      GFC_REAL_16 sum = 0, oldsum;
104      GFC_REAL_16 inv2x2 = 1 / (2 * x * x);
105      GFC_REAL_16 fac = 1;
106      int n = 1;
107
108      while (n < 200)
109	{
110	  fac *= - (2*n - 1) * inv2x2;
111	  oldsum = sum;
112	  sum += fac;
113
114	  if (sum == oldsum)
115	    break;
116
117	  n++;
118	}
119
120      return (1 + sum) / x * (_M_2_SQRTPI / 2);
121    }
122}
123
124#endif
125
126#endif
127