]> 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
re PR libfortran/33595 (FAIL: gfortran.dg/nint_2.f90 -O0 execution test)
[thirdparty/gcc.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
Mirror of https://gcc.gnu.org/git/gcc.git