libgfortran/intrinsics/c99_functions.c

   1 /* Implementation of various C99 functions
   2    Copyright (C) 2004 Free Software Foundation, Inc.
   3
   4 This file is part of the GNU Fortran 95 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 2 of the License, or (at your option) any later version.
  10
  11 In addition to the permissions in the GNU General Public License, the
  12 Free Software Foundation gives you unlimited permission to link the
  13 compiled version of this file into combinations with other programs,
  14 and to distribute those combinations without any restriction coming
  15 from the use of this file.  (The General Public License restrictions
  16 do apply in other respects; for example, they cover modification of
  17 the file, and distribution when not linked into a combine
  18 executable.)
  19
  20 Libgfortran is distributed in the hope that it will be useful,
  21 but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  23 GNU General Public License for more details.
  24
  25 You should have received a copy of the GNU General Public
  26 License along with libgfortran; see the file COPYING.  If not,
  27 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  28 Boston, MA 02111-1307, USA.  */
  29
  30 #include "config.h"
  31 #include <sys/types.h>
  32 #include <float.h>
  33 #include <math.h>
  34 #include "libgfortran.h"
  35
  36
  37 #ifndef HAVE_ACOSF
  38 float
  39 acosf(float x)
  40 {
  41   return (float) acos(x);
  42 }
  43 #endif
  44
  45 #ifndef HAVE_ASINF
  46 float
  47 asinf(float x)
  48 {
  49   return (float) asin(x);
  50 }
  51 #endif
  52
  53 #ifndef HAVE_ATAN2F
  54 float
  55 atan2f(float y, float x)
  56 {
  57   return (float) atan2(y, x);
  58 }
  59 #endif
  60
  61 #ifndef HAVE_ATANF
  62 float
  63 atanf(float x)
  64 {
  65   return (float) atan(x);
  66 }
  67 #endif
  68
  69 #ifndef HAVE_CEILF
  70 float
  71 ceilf(float x)
  72 {
  73   return (float) ceil(x);
  74 }
  75 #endif
  76
  77 #ifndef HAVE_COPYSIGNF
  78 float
  79 copysignf(float x, float y)
  80 {
  81   return (float) copysign(x, y);
  82 }
  83 #endif
  84
  85 #ifndef HAVE_COSF
  86 float
  87 cosf(float x)
  88 {
  89   return (float) cos(x);
  90 }
  91 #endif
  92
  93 #ifndef HAVE_COSHF
  94 float
  95 coshf(float x)
  96 {
  97   return (float) cosh(x);
  98 }
  99 #endif
 100
 101 #ifndef HAVE_EXPF
 102 float
 103 expf(float x)
 104 {
 105   return (float) exp(x);
 106 }
 107 #endif
 108
 109 #ifndef HAVE_FABSF
 110 float
 111 fabsf(float x)
 112 {
 113   return (float) fabs(x);
 114 }
 115 #endif
 116
 117 #ifndef HAVE_FLOORF
 118 float
 119 floorf(float x)
 120 {
 121   return (float) floor(x);
 122 }
 123 #endif
 124
 125 #ifndef HAVE_FREXPF
 126 float
 127 frexpf(float x, int *exp)
 128 {
 129   return (float) frexp(x, exp);
 130 }
 131 #endif
 132
 133 #ifndef HAVE_HYPOTF
 134 float
 135 hypotf(float x, float y)
 136 {
 137   return (float) hypot(x, y);
 138 }
 139 #endif
 140
 141 #ifndef HAVE_LOGF
 142 float
 143 logf(float x)
 144 {
 145   return (float) log(x);
 146 }
 147 #endif
 148
 149 #ifndef HAVE_LOG10F
 150 float
 151 log10f(float x)
 152 {
 153   return (float) log10(x);
 154 }
 155 #endif
 156
 157 #ifndef HAVE_SCALBNF
 158 float
 159 scalbnf(float x, int y)
 160 {
 161   return (float) scalbn(x, y);
 162 }
 163 #endif
 164
 165 #ifndef HAVE_SINF
 166 float
 167 sinf(float x)
 168 {
 169   return (float) sin(x);
 170 }
 171 #endif
 172
 173 #ifndef HAVE_SINHF
 174 float
 175 sinhf(float x)
 176 {
 177   return (float) sinh(x);
 178 }
 179 #endif
 180
 181 #ifndef HAVE_SQRTF
 182 float
 183 sqrtf(float x)
 184 {
 185   return (float) sqrt(x);
 186 }
 187 #endif
 188
 189 #ifndef HAVE_TANF
 190 float
 191 tanf(float x)
 192 {
 193   return (float) tan(x);
 194 }
 195 #endif
 196
 197 #ifndef HAVE_TANHF
 198 float
 199 tanhf(float x)
 200 {
 201   return (float) tanh(x);
 202 }
 203 #endif
 204
 205 #ifndef HAVE_NEXTAFTERF
 206 /* This is a portable implementation of nextafterf that is intended to be
 207    independent of the floating point format or its in memory representation.
 208    This implementation works correctly with denormalized values.  */
 209 float
 210 nextafterf(float x, float y)
 211 {
 212   /* This variable is marked volatile to avoid excess precision problems
 213      on some platforms, including IA-32.  */
 214   volatile float delta;
 215   float absx, denorm_min;
 216
 217   if (isnan(x) || isnan(y))
 218     return x + y;
 219   if (x == y)
 220     return x;
 221   if (!isfinite (x))
 222     return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
 223
 224   /* absx = fabsf (x);  */
 225   absx = (x < 0.0) ? -x : x;
 226
 227   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 228   if (__FLT_DENORM_MIN__ == 0.0f)
 229     denorm_min = __FLT_MIN__;
 230   else
 231     denorm_min = __FLT_DENORM_MIN__;
 232
 233   if (absx < __FLT_MIN__)
 234     delta = denorm_min;
 235   else
 236     {
 237       float frac;
 238       int exp;
 239
 240       /* Discard the fraction from x.  */
 241       frac = frexpf (absx, &exp);
 242       delta = scalbnf (0.5f, exp);
 243
 244       /* Scale x by the epsilon of the representation.  By rights we should
 245          have been able to combine this with scalbnf, but some targets don't
 246          get that correct with denormals.  */
 247       delta *= __FLT_EPSILON__;
 248
 249       /* If we're going to be reducing the absolute value of X, and doing so
 250          would reduce the exponent of X, then the delta to be applied is
 251          one exponent smaller.  */
 252       if (frac == 0.5f && (y < x) == (x > 0))
 253         delta *= 0.5f;
 254
 255       /* If that underflows to zero, then we're back to the minimum.  */
 256       if (delta == 0.0f)
 257         delta = denorm_min;
 258     }
 259
 260   if (y < x)
 261     delta = -delta;
 262
 263   return x + delta;
 264 }
 265 #endif
 266
 267
 268 #ifndef HAVE_POWF
 269 float
 270 powf(float x, float y)
 271 {
 272   return (float) pow(x, y);
 273 }
 274 #endif
 275
 276 /* Note that if fpclassify is not defined, then NaN is not handled */
 277
 278 /* Algorithm by Steven G. Kargl.  */
 279
 280 #ifndef HAVE_ROUND
 281 /* Round to nearest integral value.  If the argument is halfway between two
 282    integral values then round away from zero.  */
 283
 284 double
 285 round(double x)
 286 {
 287    double t;
 288 #if defined(fpclassify)
 289    int i;
 290    i = fpclassify(x);
 291    if (i == FP_INFINITE || i == FP_NAN)
 292      return (x);
 293 #endif
 294
 295    if (x >= 0.0)
 296     {
 297       t = ceil(x);
 298       if (t - x > 0.5)
 299         t -= 1.0;
 300       return (t);
 301     }
 302    else
 303     {
 304       t = ceil(-x);
 305       if (t + x > 0.5)
 306         t -= 1.0;
 307       return (-t);
 308     }
 309 }
 310 #endif
 311
 312 #ifndef HAVE_ROUNDF
 313 /* Round to nearest integral value.  If the argument is halfway between two
 314    integral values then round away from zero.  */
 315
 316 float
 317 roundf(float x)
 318 {
 319    float t;
 320 #if defined(fpclassify)
 321    int i;
 322
 323    i = fpclassify(x);
 324    if (i == FP_INFINITE || i == FP_NAN)
 325      return (x);
 326 #endif
 327
 328    if (x >= 0.0)
 329     {
 330       t = ceilf(x);
 331       if (t - x > 0.5)
 332         t -= 1.0;
 333       return (t);
 334     }
 335    else
 336     {
 337       t = ceilf(-x);
 338       if (t + x > 0.5)
 339         t -= 1.0;
 340       return (-t);
 341     }
 342 }
 343 #endif