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