OSDN Git Service

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