]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgfortran/intrinsics/c99_functions.c
projects / thirdparty / gcc.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Remove obsolete Tru64 UNIX V5.1B support
[thirdparty/gcc.git] / libgfortran / intrinsics / c99_functions.c
CommitLineData
a31bd8a4 1/* Implementation of various C99 functions
126b6848 2 Copyright (C) 2004, 2009, 2010, 2012 Free Software Foundation, Inc.
a31bd8a4 3
4This file is part of the GNU Fortran 95 runtime library (libgfortran).
5
6Libgfortran is free software; you can redistribute it and/or
b417ea8c 7modify it under the terms of the GNU General Public
a31bd8a4 8License as published by the Free Software Foundation; either
6bc9506f 9version 3 of the License, or (at your option) any later version.
a31bd8a4 10
11Libgfortran is distributed in the hope that it will be useful,
12but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
b417ea8c 14GNU General Public License for more details.
a31bd8a4 15
6bc9506f 16Under Section 7 of GPL version 3, you are granted additional
17permissions described in the GCC Runtime Library Exception, version
183.1, as published by the Free Software Foundation.
19
20You should have received a copy of the GNU General Public License and
21a copy of the GCC Runtime Library Exception along with this program;
22see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23<http://www.gnu.org/licenses/>. */
a31bd8a4 24
25#include "config.h"
d213114b 26
27#define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
a31bd8a4 28#include "libgfortran.h"
29
a160b48a 30/* IRIX's <math.h> declares a non-C99 compliant implementation of cabs,
31 which takes two floating point arguments instead of a single complex.
32 If <complex.h> is missing this prevents building of c99_functions.c.
33 To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */
34
35#if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
36#undef HAVE_CABS
37#undef HAVE_CABSF
38#undef HAVE_CABSL
39#define cabs __gfc_cabs
40#define cabsf __gfc_cabsf
41#define cabsl __gfc_cabsl
42#endif
43
c8f9aef9 44/* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
45 if not, we define a fallback version here. */
46#ifndef I
47# if defined(_Imaginary_I)
48# define I _Imaginary_I
49# elif defined(_Complex_I)
50# define I _Complex_I
51# else
52# define I (1.0fi)
53# endif
54#endif
7c715e7b 55
c8f9aef9 56/* Prototypes are included to silence -Wstrict-prototypes
57 -Wmissing-prototypes. */
7c715e7b 58
f29bc148 59
60/* Wrappers for systems without the various C99 single precision Bessel
61 functions. */
62
63#if defined(HAVE_J0) && ! defined(HAVE_J0F)
64#define HAVE_J0F 1
c8f9aef9 65float j0f (float);
f29bc148 66
67float
68j0f (float x)
69{
70 return (float) j0 ((double) x);
71}
72#endif
73
74#if defined(HAVE_J1) && !defined(HAVE_J1F)
75#define HAVE_J1F 1
c8f9aef9 76float j1f (float);
f29bc148 77
78float j1f (float x)
79{
80 return (float) j1 ((double) x);
81}
82#endif
83
84#if defined(HAVE_JN) && !defined(HAVE_JNF)
85#define HAVE_JNF 1
c8f9aef9 86float jnf (int, float);
f29bc148 87
88float
89jnf (int n, float x)
90{
91 return (float) jn (n, (double) x);
92}
93#endif
94
95#if defined(HAVE_Y0) && !defined(HAVE_Y0F)
96#define HAVE_Y0F 1
c8f9aef9 97float y0f (float);
f29bc148 98
99float
100y0f (float x)
101{
102 return (float) y0 ((double) x);
103}
104#endif
105
106#if defined(HAVE_Y1) && !defined(HAVE_Y1F)
107#define HAVE_Y1F 1
c8f9aef9 108float y1f (float);
f29bc148 109
110float
111y1f (float x)
112{
113 return (float) y1 ((double) x);
114}
115#endif
116
117#if defined(HAVE_YN) && !defined(HAVE_YNF)
118#define HAVE_YNF 1
c8f9aef9 119float ynf (int, float);
f29bc148 120
121float
122ynf (int n, float x)
123{
124 return (float) yn (n, (double) x);
125}
126#endif
127
128
129/* Wrappers for systems without the C99 erff() and erfcf() functions. */
130
131#if defined(HAVE_ERF) && !defined(HAVE_ERFF)
132#define HAVE_ERFF 1
c8f9aef9 133float erff (float);
f29bc148 134
135float
136erff (float x)
137{
138 return (float) erf ((double) x);
139}
140#endif
141
142#if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
143#define HAVE_ERFCF 1
c8f9aef9 144float erfcf (float);
f29bc148 145
146float
147erfcf (float x)
148{
149 return (float) erfc ((double) x);
150}
151#endif
152
153
8f838781 154#ifndef HAVE_ACOSF
d55715ee 155#define HAVE_ACOSF 1
c8f9aef9 156float acosf (float x);
157
8f838781 158float
c8f9aef9 159acosf (float x)
8f838781 160{
c8f9aef9 161 return (float) acos (x);
8f838781 162}
163#endif
164
0b4b7c0f 165#if HAVE_ACOSH && !HAVE_ACOSHF
c8f9aef9 166float acoshf (float x);
167
0b4b7c0f 168float
169acoshf (float x)
170{
171 return (float) acosh ((double) x);
172}
173#endif
174
8f838781 175#ifndef HAVE_ASINF
d55715ee 176#define HAVE_ASINF 1
c8f9aef9 177float asinf (float x);
178
8f838781 179float
c8f9aef9 180asinf (float x)
8f838781 181{
c8f9aef9 182 return (float) asin (x);
8f838781 183}
184#endif
185
0b4b7c0f 186#if HAVE_ASINH && !HAVE_ASINHF
c8f9aef9 187float asinhf (float x);
188
0b4b7c0f 189float
190asinhf (float x)
191{
192 return (float) asinh ((double) x);
193}
194#endif
195
8f838781 196#ifndef HAVE_ATAN2F
d55715ee 197#define HAVE_ATAN2F 1
c8f9aef9 198float atan2f (float y, float x);
199
8f838781 200float
c8f9aef9 201atan2f (float y, float x)
8f838781 202{
c8f9aef9 203 return (float) atan2 (y, x);
8f838781 204}
205#endif
206
207#ifndef HAVE_ATANF
d55715ee 208#define HAVE_ATANF 1
c8f9aef9 209float atanf (float x);
210
8f838781 211float
c8f9aef9 212atanf (float x)
8f838781 213{
c8f9aef9 214 return (float) atan (x);
8f838781 215}
216#endif
217
0b4b7c0f 218#if HAVE_ATANH && !HAVE_ATANHF
c8f9aef9 219float atanhf (float x);
220
0b4b7c0f 221float
222atanhf (float x)
223{
224 return (float) atanh ((double) x);
225}
226#endif
227
8f838781 228#ifndef HAVE_CEILF
d55715ee 229#define HAVE_CEILF 1
c8f9aef9 230float ceilf (float x);
231
8f838781 232float
c8f9aef9 233ceilf (float x)
8f838781 234{
c8f9aef9 235 return (float) ceil (x);
8f838781 236}
237#endif
238
239#ifndef HAVE_COPYSIGNF
d55715ee 240#define HAVE_COPYSIGNF 1
c8f9aef9 241float copysignf (float x, float y);
242
8f838781 243float
c8f9aef9 244copysignf (float x, float y)
8f838781 245{
c8f9aef9 246 return (float) copysign (x, y);
8f838781 247}
248#endif
249
250#ifndef HAVE_COSF
d55715ee 251#define HAVE_COSF 1
c8f9aef9 252float cosf (float x);
253
8f838781 254float
c8f9aef9 255cosf (float x)
8f838781 256{
c8f9aef9 257 return (float) cos (x);
8f838781 258}
259#endif
260
261#ifndef HAVE_COSHF
d55715ee 262#define HAVE_COSHF 1
c8f9aef9 263float coshf (float x);
264
8f838781 265float
c8f9aef9 266coshf (float x)
8f838781 267{
c8f9aef9 268 return (float) cosh (x);
8f838781 269}
270#endif
271
272#ifndef HAVE_EXPF
d55715ee 273#define HAVE_EXPF 1
c8f9aef9 274float expf (float x);
275
8f838781 276float
c8f9aef9 277expf (float x)
8f838781 278{
c8f9aef9 279 return (float) exp (x);
8f838781 280}
281#endif
282
334f03a1 283#ifndef HAVE_FABSF
d55715ee 284#define HAVE_FABSF 1
c8f9aef9 285float fabsf (float x);
286
334f03a1 287float
c8f9aef9 288fabsf (float x)
334f03a1 289{
c8f9aef9 290 return (float) fabs (x);
334f03a1 291}
292#endif
293
8f838781 294#ifndef HAVE_FLOORF
d55715ee 295#define HAVE_FLOORF 1
c8f9aef9 296float floorf (float x);
297
8f838781 298float
c8f9aef9 299floorf (float x)
8f838781 300{
c8f9aef9 301 return (float) floor (x);
8f838781 302}
303#endif
304
71d5b534 305#ifndef HAVE_FMODF
306#define HAVE_FMODF 1
c8f9aef9 307float fmodf (float x, float y);
308
71d5b534 309float
310fmodf (float x, float y)
311{
312 return (float) fmod (x, y);
313}
314#endif
315
8f838781 316#ifndef HAVE_FREXPF
d55715ee 317#define HAVE_FREXPF 1
c8f9aef9 318float frexpf (float x, int *exp);
319
8f838781 320float
c8f9aef9 321frexpf (float x, int *exp)
8f838781 322{
c8f9aef9 323 return (float) frexp (x, exp);
8f838781 324}
325#endif
326
327#ifndef HAVE_HYPOTF
d55715ee 328#define HAVE_HYPOTF 1
c8f9aef9 329float hypotf (float x, float y);
330
8f838781 331float
c8f9aef9 332hypotf (float x, float y)
8f838781 333{
c8f9aef9 334 return (float) hypot (x, y);
8f838781 335}
336#endif
337
338#ifndef HAVE_LOGF
d55715ee 339#define HAVE_LOGF 1
c8f9aef9 340float logf (float x);
341
8f838781 342float
c8f9aef9 343logf (float x)
8f838781 344{
c8f9aef9 345 return (float) log (x);
8f838781 346}
347#endif
348
349#ifndef HAVE_LOG10F
d55715ee 350#define HAVE_LOG10F 1
c8f9aef9 351float log10f (float x);
352
8f838781 353float
c8f9aef9 354log10f (float x)
8f838781 355{
c8f9aef9 356 return (float) log10 (x);
8f838781 357}
358#endif
359
97844b15 360#ifndef HAVE_SCALBN
d55715ee 361#define HAVE_SCALBN 1
c8f9aef9 362double scalbn (double x, int y);
363
97844b15 364double
c8f9aef9 365scalbn (double x, int y)
97844b15 366{
47b5bfdc 367#if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
368 return ldexp (x, y);
369#else
c8f9aef9 370 return x * pow (FLT_RADIX, y);
47b5bfdc 371#endif
97844b15 372}
373#endif
374
8f838781 375#ifndef HAVE_SCALBNF
d55715ee 376#define HAVE_SCALBNF 1
c8f9aef9 377float scalbnf (float x, int y);
378
8f838781 379float
c8f9aef9 380scalbnf (float x, int y)
8f838781 381{
c8f9aef9 382 return (float) scalbn (x, y);
8f838781 383}
384#endif
385
386#ifndef HAVE_SINF
d55715ee 387#define HAVE_SINF 1
c8f9aef9 388float sinf (float x);
389
8f838781 390float
c8f9aef9 391sinf (float x)
8f838781 392{
c8f9aef9 393 return (float) sin (x);
8f838781 394}
395#endif
396
397#ifndef HAVE_SINHF
d55715ee 398#define HAVE_SINHF 1
c8f9aef9 399float sinhf (float x);
400
8f838781 401float
c8f9aef9 402sinhf (float x)
8f838781 403{
c8f9aef9 404 return (float) sinh (x);
8f838781 405}
406#endif
407
408#ifndef HAVE_SQRTF
d55715ee 409#define HAVE_SQRTF 1
c8f9aef9 410float sqrtf (float x);
411
8f838781 412float
c8f9aef9 413sqrtf (float x)
8f838781 414{
c8f9aef9 415 return (float) sqrt (x);
8f838781 416}
417#endif
418
419#ifndef HAVE_TANF
d55715ee 420#define HAVE_TANF 1
c8f9aef9 421float tanf (float x);
422
8f838781 423float
c8f9aef9 424tanf (float x)
8f838781 425{
c8f9aef9 426 return (float) tan (x);
8f838781 427}
428#endif
429
430#ifndef HAVE_TANHF
d55715ee 431#define HAVE_TANHF 1
c8f9aef9 432float tanhf (float x);
433
8f838781 434float
c8f9aef9 435tanhf (float x)
8f838781 436{
c8f9aef9 437 return (float) tanh (x);
8f838781 438}
439#endif
440
50f0ca0a 441#ifndef HAVE_TRUNC
d55715ee 442#define HAVE_TRUNC 1
c8f9aef9 443double trunc (double x);
444
50f0ca0a 445double
c8f9aef9 446trunc (double x)
50f0ca0a 447{
448 if (!isfinite (x))
449 return x;
450
451 if (x < 0.0)
452 return - floor (-x);
453 else
454 return floor (x);
455}
456#endif
457
458#ifndef HAVE_TRUNCF
d55715ee 459#define HAVE_TRUNCF 1
c8f9aef9 460float truncf (float x);
461
50f0ca0a 462float
c8f9aef9 463truncf (float x)
50f0ca0a 464{
465 return (float) trunc (x);
466}
467#endif
468
8f838781 469#ifndef HAVE_NEXTAFTERF
d55715ee 470#define HAVE_NEXTAFTERF 1
8f838781 471/* This is a portable implementation of nextafterf that is intended to be
472 independent of the floating point format or its in memory representation.
b0cb0fad 473 This implementation works correctly with denormalized values. */
c8f9aef9 474float nextafterf (float x, float y);
475
8f838781 476float
c8f9aef9 477nextafterf (float x, float y)
8f838781 478{
b0cb0fad 479 /* This variable is marked volatile to avoid excess precision problems
480 on some platforms, including IA-32. */
481 volatile float delta;
482 float absx, denorm_min;
8f838781 483
c8f9aef9 484 if (isnan (x) || isnan (y))
b0cb0fad 485 return x + y;
8f838781 486 if (x == y)
487 return x;
65ce19bc 488 if (!isfinite (x))
489 return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
8f838781 490
b0cb0fad 491 /* absx = fabsf (x); */
492 absx = (x < 0.0) ? -x : x;
8f838781 493
b0cb0fad 494 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
495 if (__FLT_DENORM_MIN__ == 0.0f)
496 denorm_min = __FLT_MIN__;
497 else
498 denorm_min = __FLT_DENORM_MIN__;
499
500 if (absx < __FLT_MIN__)
501 delta = denorm_min;
8f838781 502 else
503 {
b0cb0fad 504 float frac;
505 int exp;
506
507 /* Discard the fraction from x. */
508 frac = frexpf (absx, &exp);
509 delta = scalbnf (0.5f, exp);
510
511 /* Scale x by the epsilon of the representation. By rights we should
512 have been able to combine this with scalbnf, but some targets don't
513 get that correct with denormals. */
514 delta *= __FLT_EPSILON__;
515
516 /* If we're going to be reducing the absolute value of X, and doing so
517 would reduce the exponent of X, then the delta to be applied is
518 one exponent smaller. */
519 if (frac == 0.5f && (y < x) == (x > 0))
520 delta *= 0.5f;
521
522 /* If that underflows to zero, then we're back to the minimum. */
523 if (delta == 0.0f)
524 delta = denorm_min;
8f838781 525 }
b0cb0fad 526
527 if (y < x)
528 delta = -delta;
529
530 return x + delta;
8f838781 531}
532#endif
533
a2d7a3ff 534
0dc2aff5 535#if !defined(HAVE_POWF) || defined(HAVE_BROKEN_POWF)
a2d7a3ff 536#ifndef HAVE_POWF
d55715ee 537#define HAVE_POWF 1
0dc2aff5 538#endif
c8f9aef9 539float powf (float x, float y);
540
a2d7a3ff 541float
c8f9aef9 542powf (float x, float y)
a2d7a3ff 543{
c8f9aef9 544 return (float) pow (x, y);
a2d7a3ff 545}
546#endif
547
79202f2e 548
bf09288e 549#ifndef HAVE_ROUND
550#define HAVE_ROUND 1
551/* Round to nearest integral value. If the argument is halfway between two
552 integral values then round away from zero. */
553double round (double x);
554
555double
556round (double x)
557{
558 double t;
559 if (!isfinite (x))
560 return (x);
561
562 if (x >= 0.0)
563 {
564 t = floor (x);
565 if (t - x <= -0.5)
566 t += 1.0;
567 return (t);
568 }
569 else
570 {
571 t = floor (-x);
572 if (t + x <= -0.5)
573 t += 1.0;
574 return (-t);
575 }
576}
577#endif
578
579
a31bd8a4 580/* Algorithm by Steven G. Kargl. */
581
a55a9c5d 582#if !defined(HAVE_ROUNDL)
ef080b63 583#define HAVE_ROUNDL 1
c8f9aef9 584long double roundl (long double x);
585
a55a9c5d 586#if defined(HAVE_CEILL)
ef080b63 587/* Round to nearest integral value. If the argument is halfway between two
588 integral values then round away from zero. */
589
590long double
c8f9aef9 591roundl (long double x)
ef080b63 592{
593 long double t;
594 if (!isfinite (x))
595 return (x);
596
597 if (x >= 0.0)
598 {
c8f9aef9 599 t = ceill (x);
ef080b63 600 if (t - x > 0.5)
601 t -= 1.0;
602 return (t);
603 }
604 else
605 {
c8f9aef9 606 t = ceill (-x);
ef080b63 607 if (t + x > 0.5)
608 t -= 1.0;
609 return (-t);
610 }
611}
a55a9c5d 612#else
613
614/* Poor version of roundl for system that don't have ceill. */
615long double
c8f9aef9 616roundl (long double x)
a55a9c5d 617{
618 if (x > DBL_MAX || x < -DBL_MAX)
619 {
620#ifdef HAVE_NEXTAFTERL
9bc69ba0 621 long double prechalf = nextafterl (0.5L, LDBL_MAX);
a55a9c5d 622#else
623 static long double prechalf = 0.5L;
624#endif
625 return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
626 }
627 else
628 /* Use round(). */
c8f9aef9 629 return round ((double) x);
a55a9c5d 630}
631
632#endif
ef080b63 633#endif
634
a31bd8a4 635#ifndef HAVE_ROUNDF
d55715ee 636#define HAVE_ROUNDF 1
a31bd8a4 637/* Round to nearest integral value. If the argument is halfway between two
638 integral values then round away from zero. */
c8f9aef9 639float roundf (float x);
a31bd8a4 640
641float
c8f9aef9 642roundf (float x)
a31bd8a4 643{
644 float t;
839904a0 645 if (!isfinite (x))
a31bd8a4 646 return (x);
647
648 if (x >= 0.0)
649 {
c8f9aef9 650 t = floorf (x);
1a7debab 651 if (t - x <= -0.5)
652 t += 1.0;
a31bd8a4 653 return (t);
654 }
655 else
656 {
c8f9aef9 657 t = floorf (-x);
1a7debab 658 if (t + x <= -0.5)
659 t += 1.0;
a31bd8a4 660 return (-t);
661 }
662}
663#endif
14c3c235 664
ef080b63 665
666/* lround{f,,l} and llround{f,,l} functions. */
667
668#if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
669#define HAVE_LROUNDF 1
c8f9aef9 670long int lroundf (float x);
671
ef080b63 672long int
673lroundf (float x)
674{
675 return (long int) roundf (x);
676}
677#endif
678
679#if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
680#define HAVE_LROUND 1
c8f9aef9 681long int lround (double x);
682
ef080b63 683long int
684lround (double x)
685{
686 return (long int) round (x);
687}
688#endif
689
690#if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
691#define HAVE_LROUNDL 1
c8f9aef9 692long int lroundl (long double x);
693
ef080b63 694long int
695lroundl (long double x)
696{
697 return (long long int) roundl (x);
698}
699#endif
700
701#if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
702#define HAVE_LLROUNDF 1
c8f9aef9 703long long int llroundf (float x);
704
ef080b63 705long long int
706llroundf (float x)
707{
708 return (long long int) roundf (x);
709}
710#endif
711
712#if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
713#define HAVE_LLROUND 1
c8f9aef9 714long long int llround (double x);
715
ef080b63 716long long int
717llround (double x)
718{
719 return (long long int) round (x);
720}
721#endif
722
723#if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
724#define HAVE_LLROUNDL 1
c8f9aef9 725long long int llroundl (long double x);
726
ef080b63 727long long int
728llroundl (long double x)
729{
730 return (long long int) roundl (x);
731}
732#endif
733
734
14c3c235 735#ifndef HAVE_LOG10L
d55715ee 736#define HAVE_LOG10L 1
14c3c235 737/* log10 function for long double variables. The version provided here
738 reduces the argument until it fits into a double, then use log10. */
c8f9aef9 739long double log10l (long double x);
740
14c3c235 741long double
c8f9aef9 742log10l (long double x)
14c3c235 743{
744#if LDBL_MAX_EXP > DBL_MAX_EXP
745 if (x > DBL_MAX)
746 {
747 double val;
748 int p2_result = 0;
749 if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
750 if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
751 if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
752 if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
753 if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
754 val = log10 ((double) x);
755 return (val + p2_result * .30102999566398119521373889472449302L);
756 }
757#endif
758#if LDBL_MIN_EXP < DBL_MIN_EXP
759 if (x < DBL_MIN)
760 {
761 double val;
762 int p2_result = 0;
763 if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
764 if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
765 if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
766 if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
767 if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
c8f9aef9 768 val = fabs (log10 ((double) x));
14c3c235 769 return (- val - p2_result * .30102999566398119521373889472449302L);
770 }
771#endif
772 return log10 (x);
773}
774#endif
d213114b 775
776
71d5b534 777#ifndef HAVE_FLOORL
778#define HAVE_FLOORL 1
c8f9aef9 779long double floorl (long double x);
780
71d5b534 781long double
782floorl (long double x)
783{
784 /* Zero, possibly signed. */
785 if (x == 0)
786 return x;
787
788 /* Large magnitude. */
789 if (x > DBL_MAX || x < (-DBL_MAX))
790 return x;
791
792 /* Small positive values. */
793 if (x >= 0 && x < DBL_MIN)
794 return 0;
795
796 /* Small negative values. */
797 if (x < 0 && x > (-DBL_MIN))
798 return -1;
799
800 return floor (x);
801}
802#endif
803
804
805#ifndef HAVE_FMODL
806#define HAVE_FMODL 1
c8f9aef9 807long double fmodl (long double x, long double y);
808
71d5b534 809long double
810fmodl (long double x, long double y)
811{
812 if (y == 0.0L)
813 return 0.0L;
814
815 /* Need to check that the result has the same sign as x and magnitude
816 less than the magnitude of y. */
817 return x - floorl (x / y) * y;
818}
819#endif
820
821
d213114b 822#if !defined(HAVE_CABSF)
d55715ee 823#define HAVE_CABSF 1
c8f9aef9 824float cabsf (float complex z);
825
d213114b 826float
827cabsf (float complex z)
828{
829 return hypotf (REALPART (z), IMAGPART (z));
830}
831#endif
832
833#if !defined(HAVE_CABS)
d55715ee 834#define HAVE_CABS 1
c8f9aef9 835double cabs (double complex z);
836
d213114b 837double
838cabs (double complex z)
839{
840 return hypot (REALPART (z), IMAGPART (z));
841}
842#endif
843
844#if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
d55715ee 845#define HAVE_CABSL 1
c8f9aef9 846long double cabsl (long double complex z);
847
d213114b 848long double
849cabsl (long double complex z)
850{
851 return hypotl (REALPART (z), IMAGPART (z));
852}
853#endif
854
855
856#if !defined(HAVE_CARGF)
d55715ee 857#define HAVE_CARGF 1
c8f9aef9 858float cargf (float complex z);
859
d213114b 860float
861cargf (float complex z)
862{
863 return atan2f (IMAGPART (z), REALPART (z));
864}
865#endif
866
867#if !defined(HAVE_CARG)
d55715ee 868#define HAVE_CARG 1
c8f9aef9 869double carg (double complex z);
870
d213114b 871double
872carg (double complex z)
873{
874 return atan2 (IMAGPART (z), REALPART (z));
875}
876#endif
877
878#if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
d55715ee 879#define HAVE_CARGL 1
c8f9aef9 880long double cargl (long double complex z);
881
d213114b 882long double
883cargl (long double complex z)
884{
885 return atan2l (IMAGPART (z), REALPART (z));
886}
887#endif
888
889
890/* exp(z) = exp(a)*(cos(b) + i sin(b)) */
891#if !defined(HAVE_CEXPF)
d55715ee 892#define HAVE_CEXPF 1
c8f9aef9 893float complex cexpf (float complex z);
894
d213114b 895float complex
896cexpf (float complex z)
897{
898 float a, b;
899 float complex v;
900
901 a = REALPART (z);
902 b = IMAGPART (z);
903 COMPLEX_ASSIGN (v, cosf (b), sinf (b));
904 return expf (a) * v;
905}
906#endif
907
908#if !defined(HAVE_CEXP)
d55715ee 909#define HAVE_CEXP 1
c8f9aef9 910double complex cexp (double complex z);
911
d213114b 912double complex
913cexp (double complex z)
914{
915 double a, b;
916 double complex v;
917
918 a = REALPART (z);
919 b = IMAGPART (z);
920 COMPLEX_ASSIGN (v, cos (b), sin (b));
921 return exp (a) * v;
922}
923#endif
924
925#if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
d55715ee 926#define HAVE_CEXPL 1
c8f9aef9 927long double complex cexpl (long double complex z);
928
d213114b 929long double complex
930cexpl (long double complex z)
931{
932 long double a, b;
933 long double complex v;
934
935 a = REALPART (z);
936 b = IMAGPART (z);
937 COMPLEX_ASSIGN (v, cosl (b), sinl (b));
938 return expl (a) * v;
939}
940#endif
941
942
943/* log(z) = log (cabs(z)) + i*carg(z) */
944#if !defined(HAVE_CLOGF)
d55715ee 945#define HAVE_CLOGF 1
c8f9aef9 946float complex clogf (float complex z);
947
d213114b 948float complex
949clogf (float complex z)
950{
951 float complex v;
952
953 COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
954 return v;
955}
956#endif
957
958#if !defined(HAVE_CLOG)
d55715ee 959#define HAVE_CLOG 1
c8f9aef9 960double complex clog (double complex z);
961
d213114b 962double complex
963clog (double complex z)
964{
965 double complex v;
966
967 COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
968 return v;
969}
970#endif
971
972#if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
d55715ee 973#define HAVE_CLOGL 1
c8f9aef9 974long double complex clogl (long double complex z);
975
d213114b 976long double complex
977clogl (long double complex z)
978{
979 long double complex v;
980
981 COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
982 return v;
983}
984#endif
985
986
987/* log10(z) = log10 (cabs(z)) + i*carg(z) */
988#if !defined(HAVE_CLOG10F)
d55715ee 989#define HAVE_CLOG10F 1
c8f9aef9 990float complex clog10f (float complex z);
991
d213114b 992float complex
993clog10f (float complex z)
994{
995 float complex v;
996
997 COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
998 return v;
999}
1000#endif
1001
1002#if !defined(HAVE_CLOG10)
d55715ee 1003#define HAVE_CLOG10 1
c8f9aef9 1004double complex clog10 (double complex z);
1005
d213114b 1006double complex
1007clog10 (double complex z)
1008{
1009 double complex v;
1010
1011 COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
1012 return v;
1013}
1014#endif
1015
1016#if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
d55715ee 1017#define HAVE_CLOG10L 1
c8f9aef9 1018long double complex clog10l (long double complex z);
1019
d213114b 1020long double complex
1021clog10l (long double complex z)
1022{
1023 long double complex v;
1024
1025 COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1026 return v;
1027}
1028#endif
1029
1030
1031/* pow(base, power) = cexp (power * clog (base)) */
1032#if !defined(HAVE_CPOWF)
d55715ee 1033#define HAVE_CPOWF 1
c8f9aef9 1034float complex cpowf (float complex base, float complex power);
1035
d213114b 1036float complex
1037cpowf (float complex base, float complex power)
1038{
1039 return cexpf (power * clogf (base));
1040}
1041#endif
1042
1043#if !defined(HAVE_CPOW)
d55715ee 1044#define HAVE_CPOW 1
c8f9aef9 1045double complex cpow (double complex base, double complex power);
1046
d213114b 1047double complex
1048cpow (double complex base, double complex power)
1049{
1050 return cexp (power * clog (base));
1051}
1052#endif
1053
1054#if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
d55715ee 1055#define HAVE_CPOWL 1
c8f9aef9 1056long double complex cpowl (long double complex base, long double complex power);
1057
d213114b 1058long double complex
1059cpowl (long double complex base, long double complex power)
1060{
1061 return cexpl (power * clogl (base));
1062}
1063#endif
1064
1065
1066/* sqrt(z). Algorithm pulled from glibc. */
1067#if !defined(HAVE_CSQRTF)
d55715ee 1068#define HAVE_CSQRTF 1
c8f9aef9 1069float complex csqrtf (float complex z);
1070
d213114b 1071float complex
1072csqrtf (float complex z)
1073{
1074 float re, im;
1075 float complex v;
1076
1077 re = REALPART (z);
1078 im = IMAGPART (z);
1079 if (im == 0)
1080 {
1081 if (re < 0)
1082 {
1083 COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1084 }
1085 else
1086 {
1087 COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1088 }
1089 }
1090 else if (re == 0)
1091 {
1092 float r;
1093
1094 r = sqrtf (0.5 * fabsf (im));
1095
e3126e9b 1096 COMPLEX_ASSIGN (v, r, copysignf (r, im));
d213114b 1097 }
1098 else
1099 {
1100 float d, r, s;
1101
1102 d = hypotf (re, im);
1103 /* Use the identity 2 Re res Im res = Im x
1104 to avoid cancellation error in d +/- Re x. */
1105 if (re > 0)
1106 {
1107 r = sqrtf (0.5 * d + 0.5 * re);
1108 s = (0.5 * im) / r;
1109 }
1110 else
1111 {
1112 s = sqrtf (0.5 * d - 0.5 * re);
1113 r = fabsf ((0.5 * im) / s);
1114 }
1115
1116 COMPLEX_ASSIGN (v, r, copysignf (s, im));
1117 }
1118 return v;
1119}
1120#endif
1121
1122#if !defined(HAVE_CSQRT)
d55715ee 1123#define HAVE_CSQRT 1
c8f9aef9 1124double complex csqrt (double complex z);
1125
d213114b 1126double complex
1127csqrt (double complex z)
1128{
1129 double re, im;
1130 double complex v;
1131
1132 re = REALPART (z);
1133 im = IMAGPART (z);
1134 if (im == 0)
1135 {
1136 if (re < 0)
1137 {
1138 COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1139 }
1140 else
1141 {
1142 COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1143 }
1144 }
1145 else if (re == 0)
1146 {
1147 double r;
1148
1149 r = sqrt (0.5 * fabs (im));
1150
e3126e9b 1151 COMPLEX_ASSIGN (v, r, copysign (r, im));
d213114b 1152 }
1153 else
1154 {
1155 double d, r, s;
1156
1157 d = hypot (re, im);
1158 /* Use the identity 2 Re res Im res = Im x
1159 to avoid cancellation error in d +/- Re x. */
1160 if (re > 0)
1161 {
1162 r = sqrt (0.5 * d + 0.5 * re);
1163 s = (0.5 * im) / r;
1164 }
1165 else
1166 {
1167 s = sqrt (0.5 * d - 0.5 * re);
1168 r = fabs ((0.5 * im) / s);
1169 }
1170
1171 COMPLEX_ASSIGN (v, r, copysign (s, im));
1172 }
1173 return v;
1174}
1175#endif
1176
1177#if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
d55715ee 1178#define HAVE_CSQRTL 1
c8f9aef9 1179long double complex csqrtl (long double complex z);
1180
d213114b 1181long double complex
1182csqrtl (long double complex z)
1183{
1184 long double re, im;
1185 long double complex v;
1186
1187 re = REALPART (z);
1188 im = IMAGPART (z);
1189 if (im == 0)
1190 {
1191 if (re < 0)
1192 {
1193 COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1194 }
1195 else
1196 {
1197 COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1198 }
1199 }
1200 else if (re == 0)
1201 {
1202 long double r;
1203
1204 r = sqrtl (0.5 * fabsl (im));
1205
1206 COMPLEX_ASSIGN (v, copysignl (r, im), r);
1207 }
1208 else
1209 {
1210 long double d, r, s;
1211
1212 d = hypotl (re, im);
1213 /* Use the identity 2 Re res Im res = Im x
1214 to avoid cancellation error in d +/- Re x. */
1215 if (re > 0)
1216 {
1217 r = sqrtl (0.5 * d + 0.5 * re);
1218 s = (0.5 * im) / r;
1219 }
1220 else
1221 {
1222 s = sqrtl (0.5 * d - 0.5 * re);
1223 r = fabsl ((0.5 * im) / s);
1224 }
1225
1226 COMPLEX_ASSIGN (v, r, copysignl (s, im));
1227 }
1228 return v;
1229}
1230#endif
1231
1232
1233/* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
1234#if !defined(HAVE_CSINHF)
d55715ee 1235#define HAVE_CSINHF 1
c8f9aef9 1236float complex csinhf (float complex a);
1237
d213114b 1238float complex
1239csinhf (float complex a)
1240{
1241 float r, i;
1242 float complex v;
1243
1244 r = REALPART (a);
1245 i = IMAGPART (a);
1246 COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1247 return v;
1248}
1249#endif
1250
1251#if !defined(HAVE_CSINH)
d55715ee 1252#define HAVE_CSINH 1
c8f9aef9 1253double complex csinh (double complex a);
1254
d213114b 1255double complex
1256csinh (double complex a)
1257{
1258 double r, i;
1259 double complex v;
1260
1261 r = REALPART (a);
1262 i = IMAGPART (a);
1263 COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1264 return v;
1265}
1266#endif
1267
1268#if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
d55715ee 1269#define HAVE_CSINHL 1
c8f9aef9 1270long double complex csinhl (long double complex a);
1271
d213114b 1272long double complex
1273csinhl (long double complex a)
1274{
1275 long double r, i;
1276 long double complex v;
1277
1278 r = REALPART (a);
1279 i = IMAGPART (a);
1280 COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1281 return v;
1282}
1283#endif
1284
1285
dda29e01 1286/* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */
d213114b 1287#if !defined(HAVE_CCOSHF)
d55715ee 1288#define HAVE_CCOSHF 1
c8f9aef9 1289float complex ccoshf (float complex a);
1290
d213114b 1291float complex
1292ccoshf (float complex a)
1293{
1294 float r, i;
1295 float complex v;
1296
1297 r = REALPART (a);
1298 i = IMAGPART (a);
dda29e01 1299 COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
d213114b 1300 return v;
1301}
1302#endif
1303
1304#if !defined(HAVE_CCOSH)
d55715ee 1305#define HAVE_CCOSH 1
c8f9aef9 1306double complex ccosh (double complex a);
1307
d213114b 1308double complex
1309ccosh (double complex a)
1310{
1311 double r, i;
1312 double complex v;
1313
1314 r = REALPART (a);
1315 i = IMAGPART (a);
dda29e01 1316 COMPLEX_ASSIGN (v, cosh (r) * cos (i), sinh (r) * sin (i));
d213114b 1317 return v;
1318}
1319#endif
1320
1321#if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
d55715ee 1322#define HAVE_CCOSHL 1
c8f9aef9 1323long double complex ccoshl (long double complex a);
1324
d213114b 1325long double complex
1326ccoshl (long double complex a)
1327{
1328 long double r, i;
1329 long double complex v;
1330
1331 r = REALPART (a);
1332 i = IMAGPART (a);
dda29e01 1333 COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
d213114b 1334 return v;
1335}
1336#endif
1337
1338
dda29e01 1339/* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */
d213114b 1340#if !defined(HAVE_CTANHF)
d55715ee 1341#define HAVE_CTANHF 1
c8f9aef9 1342float complex ctanhf (float complex a);
1343
d213114b 1344float complex
1345ctanhf (float complex a)
1346{
1347 float rt, it;
1348 float complex n, d;
1349
1350 rt = tanhf (REALPART (a));
1351 it = tanf (IMAGPART (a));
1352 COMPLEX_ASSIGN (n, rt, it);
dda29e01 1353 COMPLEX_ASSIGN (d, 1, rt * it);
d213114b 1354
1355 return n / d;
1356}
1357#endif
1358
1359#if !defined(HAVE_CTANH)
d55715ee 1360#define HAVE_CTANH 1
c8f9aef9 1361double complex ctanh (double complex a);
d213114b 1362double complex
1363ctanh (double complex a)
1364{
1365 double rt, it;
1366 double complex n, d;
1367
1368 rt = tanh (REALPART (a));
1369 it = tan (IMAGPART (a));
1370 COMPLEX_ASSIGN (n, rt, it);
dda29e01 1371 COMPLEX_ASSIGN (d, 1, rt * it);
d213114b 1372
1373 return n / d;
1374}
1375#endif
1376
1377#if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
d55715ee 1378#define HAVE_CTANHL 1
c8f9aef9 1379long double complex ctanhl (long double complex a);
1380
d213114b 1381long double complex
1382ctanhl (long double complex a)
1383{
1384 long double rt, it;
1385 long double complex n, d;
1386
1387 rt = tanhl (REALPART (a));
1388 it = tanl (IMAGPART (a));
1389 COMPLEX_ASSIGN (n, rt, it);
dda29e01 1390 COMPLEX_ASSIGN (d, 1, rt * it);
d213114b 1391
1392 return n / d;
1393}
1394#endif
1395
1396
1397/* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
1398#if !defined(HAVE_CSINF)
d55715ee 1399#define HAVE_CSINF 1
c8f9aef9 1400float complex csinf (float complex a);
1401
d213114b 1402float complex
1403csinf (float complex a)
1404{
1405 float r, i;
1406 float complex v;
1407
1408 r = REALPART (a);
1409 i = IMAGPART (a);
1410 COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1411 return v;
1412}
1413#endif
1414
1415#if !defined(HAVE_CSIN)
d55715ee 1416#define HAVE_CSIN 1
c8f9aef9 1417double complex csin (double complex a);
1418
d213114b 1419double complex
1420csin (double complex a)
1421{
1422 double r, i;
1423 double complex v;
1424
1425 r = REALPART (a);
1426 i = IMAGPART (a);
1427 COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1428 return v;
1429}
1430#endif
1431
1432#if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
d55715ee 1433#define HAVE_CSINL 1
c8f9aef9 1434long double complex csinl (long double complex a);
1435
d213114b 1436long double complex
1437csinl (long double complex a)
1438{
1439 long double r, i;
1440 long double complex v;
1441
1442 r = REALPART (a);
1443 i = IMAGPART (a);
1444 COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1445 return v;
1446}
1447#endif
1448
1449
1450/* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1451#if !defined(HAVE_CCOSF)
d55715ee 1452#define HAVE_CCOSF 1
c8f9aef9 1453float complex ccosf (float complex a);
1454
d213114b 1455float complex
1456ccosf (float complex a)
1457{
1458 float r, i;
1459 float complex v;
1460
1461 r = REALPART (a);
1462 i = IMAGPART (a);
1463 COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1464 return v;
1465}
1466#endif
1467
1468#if !defined(HAVE_CCOS)
d55715ee 1469#define HAVE_CCOS 1
c8f9aef9 1470double complex ccos (double complex a);
1471
d213114b 1472double complex
1473ccos (double complex a)
1474{
1475 double r, i;
1476 double complex v;
1477
1478 r = REALPART (a);
1479 i = IMAGPART (a);
1480 COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1481 return v;
1482}
1483#endif
1484
1485#if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
d55715ee 1486#define HAVE_CCOSL 1
c8f9aef9 1487long double complex ccosl (long double complex a);
1488
d213114b 1489long double complex
1490ccosl (long double complex a)
1491{
1492 long double r, i;
1493 long double complex v;
1494
1495 r = REALPART (a);
1496 i = IMAGPART (a);
1497 COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1498 return v;
1499}
1500#endif
1501
1502
1503/* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1504#if !defined(HAVE_CTANF)
d55715ee 1505#define HAVE_CTANF 1
c8f9aef9 1506float complex ctanf (float complex a);
1507
d213114b 1508float complex
1509ctanf (float complex a)
1510{
1511 float rt, it;
1512 float complex n, d;
1513
1514 rt = tanf (REALPART (a));
1515 it = tanhf (IMAGPART (a));
1516 COMPLEX_ASSIGN (n, rt, it);
1517 COMPLEX_ASSIGN (d, 1, - (rt * it));
1518
1519 return n / d;
1520}
1521#endif
1522
1523#if !defined(HAVE_CTAN)
d55715ee 1524#define HAVE_CTAN 1
c8f9aef9 1525double complex ctan (double complex a);
1526
d213114b 1527double complex
1528ctan (double complex a)
1529{
1530 double rt, it;
1531 double complex n, d;
1532
1533 rt = tan (REALPART (a));
1534 it = tanh (IMAGPART (a));
1535 COMPLEX_ASSIGN (n, rt, it);
1536 COMPLEX_ASSIGN (d, 1, - (rt * it));
1537
1538 return n / d;
1539}
1540#endif
1541
1542#if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
d55715ee 1543#define HAVE_CTANL 1
c8f9aef9 1544long double complex ctanl (long double complex a);
1545
d213114b 1546long double complex
1547ctanl (long double complex a)
1548{
1549 long double rt, it;
1550 long double complex n, d;
1551
1552 rt = tanl (REALPART (a));
1553 it = tanhl (IMAGPART (a));
1554 COMPLEX_ASSIGN (n, rt, it);
1555 COMPLEX_ASSIGN (d, 1, - (rt * it));
1556
1557 return n / d;
1558}
1559#endif
1560
6d516d29 1561
6f4274f9 1562/* Complex ASIN. Returns wrongly NaN for infinite arguments.
1563 Algorithm taken from Abramowitz & Stegun. */
1564
1565#if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1566#define HAVE_CASINF 1
c8f9aef9 1567complex float casinf (complex float z);
1568
6f4274f9 1569complex float
1570casinf (complex float z)
1571{
1572 return -I*clogf (I*z + csqrtf (1.0f-z*z));
1573}
1574#endif
1575
1576
1577#if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1578#define HAVE_CASIN 1
c8f9aef9 1579complex double casin (complex double z);
1580
6f4274f9 1581complex double
1582casin (complex double z)
1583{
1584 return -I*clog (I*z + csqrt (1.0-z*z));
1585}
1586#endif
1587
1588
1589#if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1590#define HAVE_CASINL 1
c8f9aef9 1591complex long double casinl (complex long double z);
1592
6f4274f9 1593complex long double
1594casinl (complex long double z)
1595{
1596 return -I*clogl (I*z + csqrtl (1.0L-z*z));
1597}
1598#endif
1599
1600
1601/* Complex ACOS. Returns wrongly NaN for infinite arguments.
1602 Algorithm taken from Abramowitz & Stegun. */
1603
1604#if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1605#define HAVE_CACOSF 1
c8f9aef9 1606complex float cacosf (complex float z);
1607
6f4274f9 1608complex float
1609cacosf (complex float z)
1610{
c8f9aef9 1611 return -I*clogf (z + I*csqrtf (1.0f-z*z));
6f4274f9 1612}
1613#endif
1614
1615
6f4274f9 1616#if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1617#define HAVE_CACOS 1
c8f9aef9 1618complex double cacos (complex double z);
1619
1620complex double
6f4274f9 1621cacos (complex double z)
1622{
1623 return -I*clog (z + I*csqrt (1.0-z*z));
1624}
1625#endif
1626
1627
1628#if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1629#define HAVE_CACOSL 1
c8f9aef9 1630complex long double cacosl (complex long double z);
1631
6f4274f9 1632complex long double
1633cacosl (complex long double z)
1634{
1635 return -I*clogl (z + I*csqrtl (1.0L-z*z));
1636}
1637#endif
1638
1639
1640/* Complex ATAN. Returns wrongly NaN for infinite arguments.
1641 Algorithm taken from Abramowitz & Stegun. */
1642
1643#if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1644#define HAVE_CACOSF 1
c8f9aef9 1645complex float catanf (complex float z);
1646
6f4274f9 1647complex float
1648catanf (complex float z)
1649{
1650 return I*clogf ((I+z)/(I-z))/2.0f;
1651}
1652#endif
1653
1654
1655#if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1656#define HAVE_CACOS 1
c8f9aef9 1657complex double catan (complex double z);
1658
6f4274f9 1659complex double
1660catan (complex double z)
1661{
1662 return I*clog ((I+z)/(I-z))/2.0;
1663}
1664#endif
1665
1666
1667#if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1668#define HAVE_CACOSL 1
c8f9aef9 1669complex long double catanl (complex long double z);
1670
6f4274f9 1671complex long double
1672catanl (complex long double z)
1673{
1674 return I*clogl ((I+z)/(I-z))/2.0L;
1675}
1676#endif
1677
1678
1679/* Complex ASINH. Returns wrongly NaN for infinite arguments.
1680 Algorithm taken from Abramowitz & Stegun. */
1681
1682#if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1683#define HAVE_CASINHF 1
c8f9aef9 1684complex float casinhf (complex float z);
1685
6f4274f9 1686complex float
1687casinhf (complex float z)
1688{
1689 return clogf (z + csqrtf (z*z+1.0f));
1690}
1691#endif
1692
1693
1694#if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1695#define HAVE_CASINH 1
c8f9aef9 1696complex double casinh (complex double z);
1697
6f4274f9 1698complex double
1699casinh (complex double z)
1700{
1701 return clog (z + csqrt (z*z+1.0));
1702}
1703#endif
1704
1705
1706#if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1707#define HAVE_CASINHL 1
c8f9aef9 1708complex long double casinhl (complex long double z);
1709
6f4274f9 1710complex long double
1711casinhl (complex long double z)
1712{
1713 return clogl (z + csqrtl (z*z+1.0L));
1714}
1715#endif
1716
1717
1718/* Complex ACOSH. Returns wrongly NaN for infinite arguments.
1719 Algorithm taken from Abramowitz & Stegun. */
1720
1721#if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1722#define HAVE_CACOSHF 1
c8f9aef9 1723complex float cacoshf (complex float z);
1724
6f4274f9 1725complex float
1726cacoshf (complex float z)
1727{
1728 return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1729}
1730#endif
1731
1732
1733#if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1734#define HAVE_CACOSH 1
c8f9aef9 1735complex double cacosh (complex double z);
1736
6f4274f9 1737complex double
1738cacosh (complex double z)
1739{
1740 return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1741}
1742#endif
1743
1744
1745#if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1746#define HAVE_CACOSHL 1
c8f9aef9 1747complex long double cacoshl (complex long double z);
1748
6f4274f9 1749complex long double
1750cacoshl (complex long double z)
1751{
1752 return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1753}
1754#endif
1755
1756
1757/* Complex ATANH. Returns wrongly NaN for infinite arguments.
1758 Algorithm taken from Abramowitz & Stegun. */
1759
1760#if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1761#define HAVE_CATANHF 1
c8f9aef9 1762complex float catanhf (complex float z);
1763
6f4274f9 1764complex float
1765catanhf (complex float z)
1766{
1767 return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1768}
1769#endif
1770
1771
1772#if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1773#define HAVE_CATANH 1
c8f9aef9 1774complex double catanh (complex double z);
1775
6f4274f9 1776complex double
1777catanh (complex double z)
1778{
1779 return clog ((1.0+z)/(1.0-z))/2.0;
1780}
1781#endif
1782
1783#if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1784#define HAVE_CATANHL 1
c8f9aef9 1785complex long double catanhl (complex long double z);
1786
6f4274f9 1787complex long double
1788catanhl (complex long double z)
1789{
1790 return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1791}
1792#endif
1793
1794
6d516d29 1795#if !defined(HAVE_TGAMMA)
1796#define HAVE_TGAMMA 1
c8f9aef9 1797double tgamma (double);
6d516d29 1798
1799/* Fallback tgamma() function. Uses the algorithm from
1800 http://www.netlib.org/specfun/gamma and references therein. */
1801
1802#undef SQRTPI
1803#define SQRTPI 0.9189385332046727417803297
1804
1805#undef PI
1806#define PI 3.1415926535897932384626434
1807
1808double
1809tgamma (double x)
1810{
1811 int i, n, parity;
1812 double fact, res, sum, xden, xnum, y, y1, ysq, z;
1813
1814 static double p[8] = {
1815 -1.71618513886549492533811e0, 2.47656508055759199108314e1,
1816 -3.79804256470945635097577e2, 6.29331155312818442661052e2,
1817 8.66966202790413211295064e2, -3.14512729688483675254357e4,
1818 -3.61444134186911729807069e4, 6.64561438202405440627855e4 };
1819
1820 static double q[8] = {
1821 -3.08402300119738975254353e1, 3.15350626979604161529144e2,
1822 -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1823 2.25381184209801510330112e4, 4.75584627752788110767815e3,
1824 -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1825
1826 static double c[7] = { -1.910444077728e-03,
1827 8.4171387781295e-04, -5.952379913043012e-04,
1828 7.93650793500350248e-04, -2.777777777777681622553e-03,
1829 8.333333333333333331554247e-02, 5.7083835261e-03 };
1830
1831 static const double xminin = 2.23e-308;
1832 static const double xbig = 171.624;
1833 static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1834 static double eps = 0;
1835
1836 if (eps == 0)
c8f9aef9 1837 eps = nextafter (1., 2.) - 1.;
6d516d29 1838
1839 parity = 0;
1840 fact = 1;
1841 n = 0;
1842 y = x;
1843
7e65c779 1844 if (isnan (x))
6d516d29 1845 return x;
1846
1847 if (y <= 0)
1848 {
1849 y = -x;
c8f9aef9 1850 y1 = trunc (y);
6d516d29 1851 res = y - y1;
1852
1853 if (res != 0)
1854 {
c8f9aef9 1855 if (y1 != trunc (y1*0.5l)*2)
6d516d29 1856 parity = 1;
c8f9aef9 1857 fact = -PI / sin (PI*res);
6d516d29 1858 y = y + 1;
1859 }
1860 else
1861 return x == 0 ? copysign (xinf, x) : xnan;
1862 }
1863
1864 if (y < eps)
1865 {
1866 if (y >= xminin)
1867 res = 1 / y;
1868 else
1869 return xinf;
1870 }
1871 else if (y < 13)
1872 {
1873 y1 = y;
1874 if (y < 1)
1875 {
1876 z = y;
1877 y = y + 1;
1878 }
1879 else
1880 {
1881 n = (int)y - 1;
1882 y = y - n;
1883 z = y - 1;
1884 }
1885
1886 xnum = 0;
1887 xden = 1;
1888 for (i = 0; i < 8; i++)
1889 {
1890 xnum = (xnum + p[i]) * z;
1891 xden = xden * z + q[i];
1892 }
1893
1894 res = xnum / xden + 1;
1895
1896 if (y1 < y)
1897 res = res / y1;
1898 else if (y1 > y)
1899 for (i = 1; i <= n; i++)
1900 {
1901 res = res * y;
1902 y = y + 1;
1903 }
1904 }
1905 else
1906 {
1907 if (y < xbig)
1908 {
1909 ysq = y * y;
1910 sum = c[6];
1911 for (i = 0; i < 6; i++)
1912 sum = sum / ysq + c[i];
1913
1914 sum = sum/y - y + SQRTPI;
c8f9aef9 1915 sum = sum + (y - 0.5) * log (y);
1916 res = exp (sum);
6d516d29 1917 }
1918 else
1919 return x < 0 ? xnan : xinf;
1920 }
1921
1922 if (parity)
1923 res = -res;
1924 if (fact != 1)
1925 res = fact / res;
1926
1927 return res;
1928}
1929#endif
1930
1931
1932
1933#if !defined(HAVE_LGAMMA)
1934#define HAVE_LGAMMA 1
c8f9aef9 1935double lgamma (double);
6d516d29 1936
1937/* Fallback lgamma() function. Uses the algorithm from
1938 http://www.netlib.org/specfun/algama and references therein,
1939 except for negative arguments (where netlib would return +Inf)
1940 where we use the following identity:
1941 lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1942 */
1943
1944double
1945lgamma (double y)
1946{
1947
1948#undef SQRTPI
1949#define SQRTPI 0.9189385332046727417803297
1950
1951#undef PI
1952#define PI 3.1415926535897932384626434
1953
1954#define PNT68 0.6796875
1955#define D1 -0.5772156649015328605195174
1956#define D2 0.4227843350984671393993777
1957#define D4 1.791759469228055000094023
1958
1959 static double p1[8] = {
1960 4.945235359296727046734888e0, 2.018112620856775083915565e2,
1961 2.290838373831346393026739e3, 1.131967205903380828685045e4,
1962 2.855724635671635335736389e4, 3.848496228443793359990269e4,
1963 2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1964 static double q1[8] = {
1965 6.748212550303777196073036e1, 1.113332393857199323513008e3,
1966 7.738757056935398733233834e3, 2.763987074403340708898585e4,
1967 5.499310206226157329794414e4, 6.161122180066002127833352e4,
1968 3.635127591501940507276287e4, 8.785536302431013170870835e3 };
1969 static double p2[8] = {
1970 4.974607845568932035012064e0, 5.424138599891070494101986e2,
1971 1.550693864978364947665077e4, 1.847932904445632425417223e5,
1972 1.088204769468828767498470e6, 3.338152967987029735917223e6,
1973 5.106661678927352456275255e6, 3.074109054850539556250927e6 };
1974 static double q2[8] = {
1975 1.830328399370592604055942e2, 7.765049321445005871323047e3,
1976 1.331903827966074194402448e5, 1.136705821321969608938755e6,
1977 5.267964117437946917577538e6, 1.346701454311101692290052e7,
1978 1.782736530353274213975932e7, 9.533095591844353613395747e6 };
1979 static double p4[8] = {
1980 1.474502166059939948905062e4, 2.426813369486704502836312e6,
1981 1.214755574045093227939592e8, 2.663432449630976949898078e9,
1982 2.940378956634553899906876e10, 1.702665737765398868392998e11,
1983 4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1984 static double q4[8] = {
1985 2.690530175870899333379843e3, 6.393885654300092398984238e5,
1986 4.135599930241388052042842e7, 1.120872109616147941376570e9,
1987 1.488613728678813811542398e10, 1.016803586272438228077304e11,
1988 3.417476345507377132798597e11, 4.463158187419713286462081e11 };
1989 static double c[7] = {
1990 -1.910444077728e-03, 8.4171387781295e-04,
1991 -5.952379913043012e-04, 7.93650793500350248e-04,
1992 -2.777777777777681622553e-03, 8.333333333333333331554247e-02,
1993 5.7083835261e-03 };
1994
1995 static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
1996 frtbig = 2.25e76;
1997
1998 int i;
1999 double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
2000
2001 if (eps == 0)
c8f9aef9 2002 eps = __builtin_nextafter (1., 2.) - 1.;
6d516d29 2003
2004 if ((y > 0) && (y <= xbig))
2005 {
2006 if (y <= eps)
c8f9aef9 2007 res = -log (y);
6d516d29 2008 else if (y <= 1.5)
2009 {
2010 if (y < PNT68)
2011 {
c8f9aef9 2012 corr = -log (y);
6d516d29 2013 xm1 = y;
2014 }
2015 else
2016 {
2017 corr = 0;
2018 xm1 = (y - 0.5) - 0.5;
2019 }
2020
2021 if ((y <= 0.5) || (y >= PNT68))
2022 {
2023 xden = 1;
2024 xnum = 0;
2025 for (i = 0; i < 8; i++)
2026 {
2027 xnum = xnum*xm1 + p1[i];
2028 xden = xden*xm1 + q1[i];
2029 }
2030 res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
2031 }
2032 else
2033 {
2034 xm2 = (y - 0.5) - 0.5;
2035 xden = 1;
2036 xnum = 0;
2037 for (i = 0; i < 8; i++)
2038 {
2039 xnum = xnum*xm2 + p2[i];
2040 xden = xden*xm2 + q2[i];
2041 }
2042 res = corr + xm2 * (D2 + xm2*(xnum/xden));
2043 }
2044 }
2045 else if (y <= 4)
2046 {
2047 xm2 = y - 2;
2048 xden = 1;
2049 xnum = 0;
2050 for (i = 0; i < 8; i++)
2051 {
2052 xnum = xnum*xm2 + p2[i];
2053 xden = xden*xm2 + q2[i];
2054 }
2055 res = xm2 * (D2 + xm2*(xnum/xden));
2056 }
2057 else if (y <= 12)
2058 {
2059 xm4 = y - 4;
2060 xden = -1;
2061 xnum = 0;
2062 for (i = 0; i < 8; i++)
2063 {
2064 xnum = xnum*xm4 + p4[i];
2065 xden = xden*xm4 + q4[i];
2066 }
2067 res = D4 + xm4*(xnum/xden);
2068 }
2069 else
2070 {
2071 res = 0;
2072 if (y <= frtbig)
2073 {
2074 res = c[6];
2075 ysq = y * y;
2076 for (i = 0; i < 6; i++)
2077 res = res / ysq + c[i];
2078 }
2079 res = res/y;
c8f9aef9 2080 corr = log (y);
6d516d29 2081 res = res + SQRTPI - 0.5*corr;
2082 res = res + y*(corr-1);
2083 }
2084 }
2085 else if (y < 0 && __builtin_floor (y) != y)
2086 {
2087 /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2088 For abs(y) very close to zero, we use a series expansion to
2089 the first order in y to avoid overflow. */
2090 if (y > -1.e-100)
2091 res = -2 * log (fabs (y)) - lgamma (-y);
2092 else
2093 res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2094 }
2095 else
2096 res = xinf;
2097
2098 return res;
2099}
2100#endif
2101
2102
2103#if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2104#define HAVE_TGAMMAF 1
c8f9aef9 2105float tgammaf (float);
6d516d29 2106
2107float
2108tgammaf (float x)
2109{
2110 return (float) tgamma ((double) x);
2111}
2112#endif
2113
2114#if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2115#define HAVE_LGAMMAF 1
c8f9aef9 2116float lgammaf (float);
6d516d29 2117
2118float
2119lgammaf (float x)
2120{
2121 return (float) lgamma ((double) x);
2122}
2123#endif
Mirror of https://gcc.gnu.org/git/gcc.git