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