OSDN Git Service

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