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_FLOORF
 101 float
 102 floorf(float x)
 103 {
 104   return (float) floor(x);
 105 }
 106 #endif
 107
 108 #ifndef HAVE_FREXPF
 109 float
 110 frexpf(float x, int *exp)
 111 {
 112   return (float) frexp(x, exp);
 113 }
 114 #endif
 115
 116 #ifndef HAVE_HYPOTF
 117 float
 118 hypotf(float x, float y)
 119 {
 120   return (float) hypot(x, y);
 121 }
 122 #endif
 123
 124 #ifndef HAVE_LOGF
 125 float
 126 logf(float x)
 127 {
 128   return (float) log(x);
 129 }
 130 #endif
 131
 132 #ifndef HAVE_LOG10F
 133 float
 134 log10f(float x)
 135 {
 136   return (float) log10(x);
 137 }
 138 #endif
 139
 140 #ifndef HAVE_SCALBNF
 141 float
 142 scalbnf(float x, int y)
 143 {
 144   return (float) scalbn(x, y);
 145 }
 146 #endif
 147
 148 #ifndef HAVE_SINF
 149 float
 150 sinf(float x)
 151 {
 152   return (float) sin(x);
 153 }
 154 #endif
 155
 156 #ifndef HAVE_SINHF
 157 float
 158 sinhf(float x)
 159 {
 160   return (float) sinh(x);
 161 }
 162 #endif
 163
 164 #ifndef HAVE_SQRTF
 165 float
 166 sqrtf(float x)
 167 {
 168   return (float) sqrt(x);
 169 }
 170 #endif
 171
 172 #ifndef HAVE_TANF
 173 float
 174 tanf(float x)
 175 {
 176   return (float) tan(x);
 177 }
 178 #endif
 179
 180 #ifndef HAVE_TANHF
 181 float
 182 tanhf(float x)
 183 {
 184   return (float) tanh(x);
 185 }
 186 #endif
 187
 188 #ifndef HAVE_NEXTAFTERF
 189 /* This is a portable implementation of nextafterf that is intended to be
 190    independent of the floating point format or its in memory representation.
 191    This implementation works correctly with denormalized values.  */
 192 float
 193 nextafterf(float x, float y)
 194 {
 195   /* This variable is marked volatile to avoid excess precision problems
 196      on some platforms, including IA-32.  */
 197   volatile float delta;
 198   float absx, denorm_min;
 199
 200   if (isnan(x) || isnan(y))
 201     return x + y;
 202   if (x == y)
 203     return x;
 204
 205   /* absx = fabsf (x);  */
 206   absx = (x < 0.0) ? -x : x;
 207
 208   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 209   if (__FLT_DENORM_MIN__ == 0.0f)
 210     denorm_min = __FLT_MIN__;
 211   else
 212     denorm_min = __FLT_DENORM_MIN__;
 213
 214   if (absx < __FLT_MIN__)
 215     delta = denorm_min;
 216   else
 217     {
 218       float frac;
 219       int exp;
 220
 221       /* Discard the fraction from x.  */
 222       frac = frexpf (absx, &exp);
 223       delta = scalbnf (0.5f, exp);
 224
 225       /* Scale x by the epsilon of the representation.  By rights we should
 226          have been able to combine this with scalbnf, but some targets don't
 227          get that correct with denormals.  */
 228       delta *= __FLT_EPSILON__;
 229
 230       /* If we're going to be reducing the absolute value of X, and doing so
 231          would reduce the exponent of X, then the delta to be applied is
 232          one exponent smaller.  */
 233       if (frac == 0.5f && (y < x) == (x > 0))
 234         delta *= 0.5f;
 235
 236       /* If that underflows to zero, then we're back to the minimum.  */
 237       if (delta == 0.0f)
 238         delta = denorm_min;
 239     }
 240
 241   if (y < x)
 242     delta = -delta;
 243
 244   return x + delta;
 245 }
 246 #endif
 247
 248
 249 #ifndef HAVE_POWF
 250 float
 251 powf(float x, float y)
 252 {
 253   return (float) pow(x, y);
 254 }
 255 #endif
 256
 257 /* Note that if fpclassify is not defined, then NaN is not handled */
 258
 259 /* Algorithm by Steven G. Kargl.  */
 260
 261 #ifndef HAVE_ROUND
 262 /* Round to nearest integral value.  If the argument is halfway between two
 263    integral values then round away from zero.  */
 264
 265 double
 266 round(double x)
 267 {
 268    double t;
 269 #if defined(fpclassify)
 270    int i;
 271    i = fpclassify(x);
 272    if (i == FP_INFINITE || i == FP_NAN)
 273      return (x);
 274 #endif
 275
 276    if (x >= 0.0)
 277     {
 278       t = ceil(x);
 279       if (t - x > 0.5)
 280         t -= 1.0;
 281       return (t);
 282     }
 283    else
 284     {
 285       t = ceil(-x);
 286       if (t + x > 0.5)
 287         t -= 1.0;
 288       return (-t);
 289     }
 290 }
 291 #endif
 292
 293 #ifndef HAVE_ROUNDF
 294 /* Round to nearest integral value.  If the argument is halfway between two
 295    integral values then round away from zero.  */
 296
 297 float
 298 roundf(float x)
 299 {
 300    float t;
 301 #if defined(fpclassify)
 302    int i;
 303
 304    i = fpclassify(x);
 305    if (i == FP_INFINITE || i == FP_NAN)
 306      return (x);
 307 #endif
 308
 309    if (x >= 0.0)
 310     {
 311       t = ceilf(x);
 312       if (t - x > 0.5)
 313         t -= 1.0;
 314       return (t);
 315     }
 316    else
 317     {
 318       t = ceilf(-x);
 319       if (t + x > 0.5)
 320         t -= 1.0;
 321       return (-t);
 322     }
 323 }
 324 #endif
 325