1 /* Implementation of various C99 functions
2 Copyright (C) 2004 Free Software Foundation, Inc.
4 This file is part of the GNU Fortran 95 runtime library (libgfortran).
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.
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
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.
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. */
31 #include <sys/types.h>
34 #include "libgfortran.h"
41 return (float) acos(x);
49 return (float) asin(x);
55 atan2f(float y, float x)
57 return (float) atan2(y, x);
65 return (float) atan(x);
73 return (float) ceil(x);
77 #ifndef HAVE_COPYSIGNF
79 copysignf(float x, float y)
81 return (float) copysign(x, y);
89 return (float) cos(x);
97 return (float) cosh(x);
105 return (float) exp(x);
113 return (float) fabs(x);
121 return (float) floor(x);
127 frexpf(float x, int *exp)
129 return (float) frexp(x, exp);
135 hypotf(float x, float y)
137 return (float) hypot(x, y);
145 return (float) log(x);
153 return (float) log10(x);
159 scalbn(double x, int y)
161 return x * pow(FLT_RADIX, y);
167 scalbnf(float x, int y)
169 return (float) scalbn(x, y);
177 return (float) sin(x);
185 return (float) sinh(x);
193 return (float) sqrt(x);
201 return (float) tan(x);
209 return (float) tanh(x);
231 return (float) trunc (x);
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. */
240 nextafterf(float x, float y)
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;
247 if (isnan(x) || isnan(y))
252 return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
254 /* absx = fabsf (x); */
255 absx = (x < 0.0) ? -x : x;
257 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
258 if (__FLT_DENORM_MIN__ == 0.0f)
259 denorm_min = __FLT_MIN__;
261 denorm_min = __FLT_DENORM_MIN__;
263 if (absx < __FLT_MIN__)
270 /* Discard the fraction from x. */
271 frac = frexpf (absx, &exp);
272 delta = scalbnf (0.5f, exp);
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__;
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))
285 /* If that underflows to zero, then we're back to the minimum. */
300 powf(float x, float y)
302 return (float) pow(x, y);
306 /* Note that if fpclassify is not defined, then NaN is not handled */
308 /* Algorithm by Steven G. Kargl. */
311 /* Round to nearest integral value. If the argument is halfway between two
312 integral values then round away from zero. */
318 #if defined(fpclassify)
321 if (i == FP_INFINITE || i == FP_NAN)
343 /* Round to nearest integral value. If the argument is halfway between two
344 integral values then round away from zero. */
350 #if defined(fpclassify)
354 if (i == FP_INFINITE || i == FP_NAN)
376 /* log10 function for long double variables. The version provided here
377 reduces the argument until it fits into a double, then use log10. */
379 log10l(long double x)
381 #if LDBL_MAX_EXP > DBL_MAX_EXP
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);
395 #if LDBL_MIN_EXP < DBL_MIN_EXP
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);