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., 51 Franklin Street, Fifth Floor,
  28 Boston, MA 02110-1301, 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_SCALBN
 158 double
 159 scalbn(double x, int y)
 160 {
 161   return x * pow(FLT_RADIX, y);
 162 }
 163 #endif
 164
 165 #ifndef HAVE_SCALBNF
 166 float
 167 scalbnf(float x, int y)
 168 {
 169   return (float) scalbn(x, y);
 170 }
 171 #endif
 172
 173 #ifndef HAVE_SINF
 174 float
 175 sinf(float x)
 176 {
 177   return (float) sin(x);
 178 }
 179 #endif
 180
 181 #ifndef HAVE_SINHF
 182 float
 183 sinhf(float x)
 184 {
 185   return (float) sinh(x);
 186 }
 187 #endif
 188
 189 #ifndef HAVE_SQRTF
 190 float
 191 sqrtf(float x)
 192 {
 193   return (float) sqrt(x);
 194 }
 195 #endif
 196
 197 #ifndef HAVE_TANF
 198 float
 199 tanf(float x)
 200 {
 201   return (float) tan(x);
 202 }
 203 #endif
 204
 205 #ifndef HAVE_TANHF
 206 float
 207 tanhf(float x)
 208 {
 209   return (float) tanh(x);
 210 }
 211 #endif
 212
 213 #ifndef HAVE_TRUNC
 214 double
 215 trunc(double x)
 216 {
 217   if (!isfinite (x))
 218     return x;
 219
 220   if (x < 0.0)
 221     return - floor (-x);
 222   else
 223     return floor (x);
 224 }
 225 #endif
 226
 227 #ifndef HAVE_TRUNCF
 228 float
 229 truncf(float x)
 230 {
 231   return (float) trunc (x);
 232 }
 233 #endif
 234
 235 #ifndef HAVE_NEXTAFTERF
 236 /* This is a portable implementation of nextafterf that is intended to be
 237    independent of the floating point format or its in memory representation.
 238    This implementation works correctly with denormalized values.  */
 239 float
 240 nextafterf(float x, float y)
 241 {
 242   /* This variable is marked volatile to avoid excess precision problems
 243      on some platforms, including IA-32.  */
 244   volatile float delta;
 245   float absx, denorm_min;
 246
 247   if (isnan(x) || isnan(y))
 248     return x + y;
 249   if (x == y)
 250     return x;
 251   if (!isfinite (x))
 252     return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
 253
 254   /* absx = fabsf (x);  */
 255   absx = (x < 0.0) ? -x : x;
 256
 257   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 258   if (__FLT_DENORM_MIN__ == 0.0f)
 259     denorm_min = __FLT_MIN__;
 260   else
 261     denorm_min = __FLT_DENORM_MIN__;
 262
 263   if (absx < __FLT_MIN__)
 264     delta = denorm_min;
 265   else
 266     {
 267       float frac;
 268       int exp;
 269
 270       /* Discard the fraction from x.  */
 271       frac = frexpf (absx, &exp);
 272       delta = scalbnf (0.5f, exp);
 273
 274       /* Scale x by the epsilon of the representation.  By rights we should
 275          have been able to combine this with scalbnf, but some targets don't
 276          get that correct with denormals.  */
 277       delta *= __FLT_EPSILON__;
 278
 279       /* If we're going to be reducing the absolute value of X, and doing so
 280          would reduce the exponent of X, then the delta to be applied is
 281          one exponent smaller.  */
 282       if (frac == 0.5f && (y < x) == (x > 0))
 283         delta *= 0.5f;
 284
 285       /* If that underflows to zero, then we're back to the minimum.  */
 286       if (delta == 0.0f)
 287         delta = denorm_min;
 288     }
 289
 290   if (y < x)
 291     delta = -delta;
 292
 293   return x + delta;
 294 }
 295 #endif
 296
 297
 298 #ifndef HAVE_POWF
 299 float
 300 powf(float x, float y)
 301 {
 302   return (float) pow(x, y);
 303 }
 304 #endif
 305
 306 /* Note that if fpclassify is not defined, then NaN is not handled */
 307
 308 /* Algorithm by Steven G. Kargl.  */
 309
 310 #ifndef HAVE_ROUND
 311 /* Round to nearest integral value.  If the argument is halfway between two
 312    integral values then round away from zero.  */
 313
 314 double
 315 round(double x)
 316 {
 317    double t;
 318 #if defined(fpclassify)
 319    int i;
 320    i = fpclassify(x);
 321    if (i == FP_INFINITE || i == FP_NAN)
 322      return (x);
 323 #endif
 324
 325    if (x >= 0.0)
 326     {
 327       t = ceil(x);
 328       if (t - x > 0.5)
 329         t -= 1.0;
 330       return (t);
 331     }
 332    else
 333     {
 334       t = ceil(-x);
 335       if (t + x > 0.5)
 336         t -= 1.0;
 337       return (-t);
 338     }
 339 }
 340 #endif
 341
 342 #ifndef HAVE_ROUNDF
 343 /* Round to nearest integral value.  If the argument is halfway between two
 344    integral values then round away from zero.  */
 345
 346 float
 347 roundf(float x)
 348 {
 349    float t;
 350 #if defined(fpclassify)
 351    int i;
 352
 353    i = fpclassify(x);
 354    if (i == FP_INFINITE || i == FP_NAN)
 355      return (x);
 356 #endif
 357
 358    if (x >= 0.0)
 359     {
 360       t = ceilf(x);
 361       if (t - x > 0.5)
 362         t -= 1.0;
 363       return (t);
 364     }
 365    else
 366     {
 367       t = ceilf(-x);
 368       if (t + x > 0.5)
 369         t -= 1.0;
 370       return (-t);
 371     }
 372 }
 373 #endif
 374
 375 #ifndef HAVE_LOG10L
 376 /* log10 function for long double variables. The version provided here
 377    reduces the argument until it fits into a double, then use log10.  */
 378 long double
 379 log10l(long double x)
 380 {
 381 #if LDBL_MAX_EXP > DBL_MAX_EXP
 382   if (x > DBL_MAX)
 383     {
 384       double val;
 385       int p2_result = 0;
 386       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 387       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 388       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 389       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 390       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 391       val = log10 ((double) x);
 392       return (val + p2_result * .30102999566398119521373889472449302L);
 393     }
 394 #endif
 395 #if LDBL_MIN_EXP < DBL_MIN_EXP
 396   if (x < DBL_MIN)
 397     {
 398       double val;
 399       int p2_result = 0;
 400       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 401       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 402       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 403       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 404       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 405       val = fabs(log10 ((double) x));
 406       return (- val - p2_result * .30102999566398119521373889472449302L);
 407     }
 408 #endif
 409     return log10 (x);
 410 }
 411 #endif