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