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