OSDN Git Service

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