OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgfortran/] [intrinsics/] [erfc_scaled_inc.c] - Blame information for rev 733

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 733 jeremybenn
/* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c
2
   Copyright (c) 2008, 2010 Free Software Foundation, Inc.
3
 
4
This file is part of the GNU Fortran runtime library (libgfortran).
5
 
6
Libgfortran is free software; you can redistribute it and/or
7
modify it under the terms of the GNU General Public
8
License as published by the Free Software Foundation; either
9
version 3 of the License, or (at your option) any later version.
10
 
11
Libgfortran is distributed in the hope that it will be useful,
12
but WITHOUT ANY WARRANTY; without even the implied warranty of
13
MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE.  See the
14
GNU General Public License for more details.
15
 
16
Under Section 7 of GPL version 3, you are granted additional
17
permissions described in the GCC Runtime Library Exception, version
18
3.1, as published by the Free Software Foundation.
19
 
20
You should have received a copy of the GNU General Public License and
21
a copy of the GCC Runtime Library Exception along with this program;
22
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23
<http://www.gnu.org/licenses/>.  */
24
 
25
/* This implementation of ERFC_SCALED is based on the netlib algorithm
26
   available at http://www.netlib.org/specfun/erf  */
27
 
28
#define TYPE KIND_SUFFIX(GFC_REAL_,KIND)
29
#define CONCAT(x,y) x ## y
30
#define KIND_SUFFIX(x,y) CONCAT(x,y)
31
 
32
#if (KIND == 4)
33
 
34
# define EXP(x) expf(x)
35
# define TRUNC(x) truncf(x)
36
 
37
#elif (KIND == 8)
38
 
39
# define EXP(x) exp(x)
40
# define TRUNC(x) trunc(x)
41
 
42
#elif (KIND == 10) || (KIND == 16 && defined(GFC_REAL_16_IS_LONG_DOUBLE))
43
 
44
# ifdef HAVE_EXPL
45
#  define EXP(x) expl(x)
46
# endif
47
# ifdef HAVE_TRUNCL
48
#  define TRUNC(x) truncl(x)
49
# endif
50
 
51
#elif (KIND == 16 && defined(GFC_REAL_16_IS_FLOAT128))
52
 
53
#  define EXP(x) expq(x)
54
#  define TRUNC(x) truncq(x)
55
 
56
#else
57
 
58
# error "What exactly is it that you want me to do?"
59
 
60
#endif
61
 
62
#if defined(EXP) && defined(TRUNC)
63
 
64
extern TYPE KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE);
65
export_proto(KIND_SUFFIX(erfc_scaled_r,KIND));
66
 
67
TYPE
68
KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE x)
69
{
70
  /* The main computation evaluates near-minimax approximations
71
     from "Rational Chebyshev approximations for the error function"
72
     by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
73
     transportable program uses rational functions that theoretically
74
     approximate  erf(x)  and  erfc(x)  to at least 18 significant
75
     decimal digits.  The accuracy achieved depends on the arithmetic
76
     system, the compiler, the intrinsic functions, and proper
77
     selection of the machine-dependent constants.  */
78
 
79
  int i;
80
  TYPE del, res, xden, xnum, y, ysq;
81
 
82
#if (KIND == 4)
83
  static TYPE xneg = -9.382, xsmall = 5.96e-8,
84
              xbig = 9.194, xhuge = 2.90e+3, xmax = 4.79e+37;
85
#else
86
  static TYPE xneg = -26.628, xsmall = 1.11e-16,
87
              xbig = 26.543, xhuge = 6.71e+7, xmax = 2.53e+307;
88
#endif
89
 
90
#define SQRPI ((TYPE) 0.56418958354775628695L)
91
#define THRESH ((TYPE) 0.46875L)
92
 
93
  static TYPE a[5] = { 3.16112374387056560l, 113.864154151050156l,
94
    377.485237685302021l, 3209.37758913846947l, 0.185777706184603153l };
95
 
96
  static TYPE b[4] = { 23.6012909523441209l, 244.024637934444173l,
97
    1282.61652607737228l, 2844.23683343917062l };
98
 
99
  static TYPE c[9] = { 0.564188496988670089l, 8.88314979438837594l,
100
    66.1191906371416295l, 298.635138197400131l, 881.952221241769090l,
101
    1712.04761263407058l, 2051.07837782607147l, 1230.33935479799725l,
102
    2.15311535474403846e-8l };
103
 
104
  static TYPE d[8] = { 15.7449261107098347l, 117.693950891312499l,
105
    537.181101862009858l, 1621.38957456669019l, 3290.79923573345963l,
106
    4362.61909014324716l, 3439.36767414372164l, 1230.33935480374942l };
107
 
108
  static TYPE p[6] = { 0.305326634961232344l, 0.360344899949804439l,
109
    0.125781726111229246l, 0.0160837851487422766l,
110
    0.000658749161529837803l, 0.0163153871373020978l };
111
 
112
  static TYPE q[5] = { 2.56852019228982242l, 1.87295284992346047l,
113
    0.527905102951428412l, 0.0605183413124413191l,
114
    0.00233520497626869185l };
115
 
116
  y = (x > 0 ? x : -x);
117
  if (y <= THRESH)
118
    {
119
      ysq = 0;
120
      if (y > xsmall)
121
        ysq = y * y;
122
      xnum = a[4]*ysq;
123
      xden = ysq;
124
      for (i = 0; i <= 2; i++)
125
        {
126
          xnum = (xnum + a[i]) * ysq;
127
          xden = (xden + b[i]) * ysq;
128
        }
129
      res = x * (xnum + a[3]) / (xden + b[3]);
130
      res = 1 - res;
131
      res = EXP(ysq) * res;
132
      return res;
133
    }
134
  else if (y <= 4)
135
    {
136
      xnum = c[8]*y;
137
      xden = y;
138
      for (i = 0; i <= 6; i++)
139
        {
140
          xnum = (xnum + c[i]) * y;
141
          xden = (xden + d[i]) * y;
142
        }
143
      res = (xnum + c[7]) / (xden + d[7]);
144
    }
145
  else
146
    {
147
      res = 0;
148
      if (y >= xbig)
149
        {
150
          if (y >= xmax)
151
            goto finish;
152
          if (y >= xhuge)
153
            {
154
              res = SQRPI / y;
155
              goto finish;
156
            }
157
        }
158
      ysq = ((TYPE) 1) / (y * y);
159
      xnum = p[5]*ysq;
160
      xden = ysq;
161
      for (i = 0; i <= 3; i++)
162
        {
163
          xnum = (xnum + p[i]) * ysq;
164
          xden = (xden + q[i]) * ysq;
165
        }
166
      res = ysq *(xnum + p[4]) / (xden + q[4]);
167
      res = (SQRPI -  res) / y;
168
    }
169
 
170
finish:
171
  if (x < 0)
172
    {
173
      if (x < xneg)
174
        res = __builtin_inf ();
175
      else
176
        {
177
          ysq = TRUNC (x*((TYPE) 16))/((TYPE) 16);
178
          del = (x-ysq)*(x+ysq);
179
          y = EXP(ysq*ysq) * EXP(del);
180
          res = (y+y) - res;
181
        }
182
    }
183
  return res;
184
}
185
 
186
#endif
187
 
188
#undef EXP
189
#undef TRUNC
190
 
191
#undef CONCAT
192
#undef TYPE
193
#undef KIND_SUFFIX

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.