diff options
Diffstat (limited to 'math')
| -rw-r--r-- | math/Makefile | 8 | ||||
| -rw-r--r-- | math/libm-test.inc | 135 | ||||
| -rw-r--r-- | math/math_private.h | 392 | ||||
| -rw-r--r-- | math/s_cacosh.c | 12 | ||||
| -rw-r--r-- | math/s_cacoshf.c | 14 | ||||
| -rw-r--r-- | math/s_cacoshl.c | 14 | ||||
| -rw-r--r-- | math/s_casinh.c | 7 | ||||
| -rw-r--r-- | math/s_casinhf.c | 7 | ||||
| -rw-r--r-- | math/s_casinhl.c | 7 | ||||
| -rw-r--r-- | math/s_csqrt.c | 26 | ||||
| -rw-r--r-- | math/s_csqrtf.c | 26 | ||||
| -rw-r--r-- | math/s_csqrtl.c | 26 |
12 files changed, 251 insertions, 423 deletions
diff --git a/math/Makefile b/math/Makefile index 6a6688dd0f..3a671aac71 100644 --- a/math/Makefile +++ b/math/Makefile @@ -26,9 +26,6 @@ headers := math.h bits/mathcalls.h bits/mathinline.h bits/huge_val.h \ bits/fenv.h bits/fenvinline.h bits/mathdef.h tgmath.h \ bits/math-finite.h -# Internal header files. -distribute := math_ldbl.h math_private.h machine/asm.h - # FPU support code. aux := setfpucw fpu_control @@ -84,7 +81,6 @@ calls = s_isinf s_isnan s_finite s_copysign s_modf s_scalbn s_frexp s_ldexp \ generated += $(foreach s,.c .S l.c l.S f.c f.S,$(calls:s_%=m_%$s)) routines = $(calls) $(calls:=f) $(long-c-$(long-double-fcts)) long-c-yes = $(calls:=l) -distribute += $(filter-out $(generated),$(long-m-yes:=.c) $(long-c-yes:=.c)) # Rules for the test suite. tests = test-matherr test-fenv atest-exp atest-sincos atest-exp2 basic-test \ @@ -94,7 +90,6 @@ tests = test-matherr test-fenv atest-exp atest-sincos atest-exp2 basic-test \ # We do the `long double' tests only if this data type is available and # distinct from `double'. test-longdouble-yes = test-ldouble test-ildoubl -distribute += $(test-longdouble-yes:=.c) ifneq (no,$(PERL)) libm-tests = test-float test-double $(test-longdouble-$(long-double-fcts)) \ @@ -139,15 +134,12 @@ CPPFLAGS-test-idouble.c = -U__LIBC_INTERNAL_MATH_INLINES -D__FAST_MATH__ \ CPPFLAGS-test-ildoubl.c = -U__LIBC_INTERNAL_MATH_INLINES -D__FAST_MATH__ \ -DTEST_FAST_MATH -fno-builtin -distribute += libm-test.inc gen-libm-test.pl README.libm-test - # The -lieee module sets the _LIB_VERSION_ switch to IEEE mode # for error handling in the -lm functions. install-lib += libieee.a non-lib.a += libieee.a extra-objs += libieee.a ieee-math.o -distribute += ieee-math.c include ../Rules diff --git a/math/libm-test.inc b/math/libm-test.inc index a51f125198..aeba1ae7be 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -1065,6 +1065,32 @@ cacos_test (void) TEST_c_c (cacos, nan_value, nan_value, nan_value, nan_value); + TEST_c_c (cacos, plus_zero, -1.5L, M_PI_2l, 1.194763217287109304111930828519090523536L); + TEST_c_c (cacos, minus_zero, -1.5L, M_PI_2l, 1.194763217287109304111930828519090523536L); + TEST_c_c (cacos, plus_zero, -1.0L, M_PI_2l, 0.8813735870195430252326093249797923090282L); + TEST_c_c (cacos, minus_zero, -1.0L, M_PI_2l, 0.8813735870195430252326093249797923090282L); + TEST_c_c (cacos, plus_zero, -0.5L, M_PI_2l, 0.4812118250596034474977589134243684231352L); + TEST_c_c (cacos, minus_zero, -0.5L, M_PI_2l, 0.4812118250596034474977589134243684231352L); + TEST_c_c (cacos, plus_zero, 0.5L, M_PI_2l, -0.4812118250596034474977589134243684231352L); + TEST_c_c (cacos, minus_zero, 0.5L, M_PI_2l, -0.4812118250596034474977589134243684231352L); + TEST_c_c (cacos, plus_zero, 1.0L, M_PI_2l, -0.8813735870195430252326093249797923090282L); + TEST_c_c (cacos, minus_zero, 1.0L, M_PI_2l, -0.8813735870195430252326093249797923090282L); + TEST_c_c (cacos, plus_zero, 1.5L, M_PI_2l, -1.194763217287109304111930828519090523536L); + TEST_c_c (cacos, minus_zero, 1.5L, M_PI_2l, -1.194763217287109304111930828519090523536L); + + TEST_c_c (cacos, -1.5L, plus_zero, M_PIl, -0.9624236501192068949955178268487368462704L); + TEST_c_c (cacos, -1.5L, minus_zero, M_PIl, 0.9624236501192068949955178268487368462704L); + TEST_c_c (cacos, -1.0L, plus_zero, M_PIl, minus_zero); + TEST_c_c (cacos, -1.0L, minus_zero, M_PIl, plus_zero); + TEST_c_c (cacos, -0.5L, plus_zero, 2.094395102393195492308428922186335256131L, minus_zero); + TEST_c_c (cacos, -0.5L, minus_zero, 2.094395102393195492308428922186335256131L, plus_zero); + TEST_c_c (cacos, 0.5L, plus_zero, 1.047197551196597746154214461093167628066L, minus_zero); + TEST_c_c (cacos, 0.5L, minus_zero, 1.047197551196597746154214461093167628066L, plus_zero); + TEST_c_c (cacos, 1.0L, plus_zero, plus_zero, minus_zero); + TEST_c_c (cacos, 1.0L, minus_zero, plus_zero, plus_zero); + TEST_c_c (cacos, 1.5L, plus_zero, plus_zero, -0.9624236501192068949955178268487368462704L); + TEST_c_c (cacos, 1.5L, minus_zero, plus_zero, 0.9624236501192068949955178268487368462704L); + TEST_c_c (cacos, 0.75L, 1.25L, 1.11752014915610270578240049553777969L, -1.13239363160530819522266333696834467L); TEST_c_c (cacos, -2, -3, 2.1414491111159960199416055713254211L, 1.9833870299165354323470769028940395L); @@ -1127,6 +1153,32 @@ cacosh_test (void) TEST_c_c (cacosh, nan_value, nan_value, nan_value, nan_value); + TEST_c_c (cacosh, plus_zero, -1.5L, 1.194763217287109304111930828519090523536L, -M_PI_2l); + TEST_c_c (cacosh, minus_zero, -1.5L, 1.194763217287109304111930828519090523536L, -M_PI_2l); + TEST_c_c (cacosh, plus_zero, -1.0L, 0.8813735870195430252326093249797923090282L, -M_PI_2l); + TEST_c_c (cacosh, minus_zero, -1.0L, 0.8813735870195430252326093249797923090282L, -M_PI_2l); + TEST_c_c (cacosh, plus_zero, -0.5L, 0.4812118250596034474977589134243684231352L, -M_PI_2l); + TEST_c_c (cacosh, minus_zero, -0.5L, 0.4812118250596034474977589134243684231352L, -M_PI_2l); + TEST_c_c (cacosh, plus_zero, 0.5L, 0.4812118250596034474977589134243684231352L, M_PI_2l); + TEST_c_c (cacosh, minus_zero, 0.5L, 0.4812118250596034474977589134243684231352L, M_PI_2l); + TEST_c_c (cacosh, plus_zero, 1.0L, 0.8813735870195430252326093249797923090282L, M_PI_2l); + TEST_c_c (cacosh, minus_zero, 1.0L, 0.8813735870195430252326093249797923090282L, M_PI_2l); + TEST_c_c (cacosh, plus_zero, 1.5L, 1.194763217287109304111930828519090523536L, M_PI_2l); + TEST_c_c (cacosh, minus_zero, 1.5L, 1.194763217287109304111930828519090523536L, M_PI_2l); + + TEST_c_c (cacosh, -1.5L, plus_zero, 0.9624236501192068949955178268487368462704L, M_PIl); + TEST_c_c (cacosh, -1.5L, minus_zero, 0.9624236501192068949955178268487368462704L, -M_PIl); + TEST_c_c (cacosh, -1.0L, plus_zero, plus_zero, M_PIl); + TEST_c_c (cacosh, -1.0L, minus_zero, plus_zero, -M_PIl); + TEST_c_c (cacosh, -0.5L, plus_zero, plus_zero, 2.094395102393195492308428922186335256131L); + TEST_c_c (cacosh, -0.5L, minus_zero, plus_zero, -2.094395102393195492308428922186335256131L); + TEST_c_c (cacosh, 0.5L, plus_zero, plus_zero, 1.047197551196597746154214461093167628066L); + TEST_c_c (cacosh, 0.5L, minus_zero, plus_zero, -1.047197551196597746154214461093167628066L); + TEST_c_c (cacosh, 1.0L, plus_zero, plus_zero, plus_zero); + TEST_c_c (cacosh, 1.0L, minus_zero, plus_zero, minus_zero); + TEST_c_c (cacosh, 1.5L, plus_zero, 0.9624236501192068949955178268487368462704L, plus_zero); + TEST_c_c (cacosh, 1.5L, minus_zero, 0.9624236501192068949955178268487368462704L, minus_zero); + TEST_c_c (cacosh, 0.75L, 1.25L, 1.13239363160530819522266333696834467L, 1.11752014915610270578240049553777969L); TEST_c_c (cacosh, -2, -3, 1.9833870299165354323470769028940395L, -2.1414491111159960199416055713254211L); @@ -1258,6 +1310,32 @@ casin_test (void) TEST_c_c (casin, nan_value, nan_value, nan_value, nan_value); + TEST_c_c (casin, plus_zero, -1.5L, plus_zero, -1.194763217287109304111930828519090523536L); + TEST_c_c (casin, minus_zero, -1.5L, minus_zero, -1.194763217287109304111930828519090523536L); + TEST_c_c (casin, plus_zero, -1.0L, plus_zero, -0.8813735870195430252326093249797923090282L); + TEST_c_c (casin, minus_zero, -1.0L, minus_zero, -0.8813735870195430252326093249797923090282L); + TEST_c_c (casin, plus_zero, -0.5L, plus_zero, -0.4812118250596034474977589134243684231352L); + TEST_c_c (casin, minus_zero, -0.5L, minus_zero, -0.4812118250596034474977589134243684231352L); + TEST_c_c (casin, plus_zero, 0.5L, plus_zero, 0.4812118250596034474977589134243684231352L); + TEST_c_c (casin, minus_zero, 0.5L, minus_zero, 0.4812118250596034474977589134243684231352L); + TEST_c_c (casin, plus_zero, 1.0L, plus_zero, 0.8813735870195430252326093249797923090282L); + TEST_c_c (casin, minus_zero, 1.0L, minus_zero, 0.8813735870195430252326093249797923090282L); + TEST_c_c (casin, plus_zero, 1.5L, plus_zero, 1.194763217287109304111930828519090523536L); + TEST_c_c (casin, minus_zero, 1.5L, minus_zero, 1.194763217287109304111930828519090523536L); + + TEST_c_c (casin, -1.5L, plus_zero, -M_PI_2l, 0.9624236501192068949955178268487368462704L); + TEST_c_c (casin, -1.5L, minus_zero, -M_PI_2l, -0.9624236501192068949955178268487368462704L); + TEST_c_c (casin, -1.0L, plus_zero, -M_PI_2l, plus_zero); + TEST_c_c (casin, -1.0L, minus_zero, -M_PI_2l, minus_zero); + TEST_c_c (casin, -0.5L, plus_zero, -0.5235987755982988730771072305465838140329L, plus_zero); + TEST_c_c (casin, -0.5L, minus_zero, -0.5235987755982988730771072305465838140329L, minus_zero); + TEST_c_c (casin, 0.5L, plus_zero, 0.5235987755982988730771072305465838140329L, plus_zero); + TEST_c_c (casin, 0.5L, minus_zero, 0.5235987755982988730771072305465838140329L, minus_zero); + TEST_c_c (casin, 1.0L, plus_zero, M_PI_2l, plus_zero); + TEST_c_c (casin, 1.0L, minus_zero, M_PI_2l, minus_zero); + TEST_c_c (casin, 1.5L, plus_zero, M_PI_2l, 0.9624236501192068949955178268487368462704L); + TEST_c_c (casin, 1.5L, minus_zero, M_PI_2l, -0.9624236501192068949955178268487368462704L); + TEST_c_c (casin, 0.75L, 1.25L, 0.453276177638793913448921196101971749L, 1.13239363160530819522266333696834467L); TEST_c_c (casin, -2, -3, -0.57065278432109940071028387968566963L, -1.9833870299165354323470769028940395L); @@ -1322,6 +1400,32 @@ casinh_test (void) TEST_c_c (casinh, nan_value, nan_value, nan_value, nan_value); + TEST_c_c (casinh, plus_zero, -1.5L, 0.9624236501192068949955178268487368462704L, -M_PI_2l); + TEST_c_c (casinh, minus_zero, -1.5L, -0.9624236501192068949955178268487368462704L, -M_PI_2l); + TEST_c_c (casinh, plus_zero, -1.0L, plus_zero, -M_PI_2l); + TEST_c_c (casinh, minus_zero, -1.0L, minus_zero, -M_PI_2l); + TEST_c_c (casinh, plus_zero, -0.5L, plus_zero, -0.5235987755982988730771072305465838140329L); + TEST_c_c (casinh, minus_zero, -0.5L, minus_zero, -0.5235987755982988730771072305465838140329L); + TEST_c_c (casinh, plus_zero, 0.5L, plus_zero, 0.5235987755982988730771072305465838140329L); + TEST_c_c (casinh, minus_zero, 0.5L, minus_zero, 0.5235987755982988730771072305465838140329L); + TEST_c_c (casinh, plus_zero, 1.0L, plus_zero, M_PI_2l); + TEST_c_c (casinh, minus_zero, 1.0L, minus_zero, M_PI_2l); + TEST_c_c (casinh, plus_zero, 1.5L, 0.9624236501192068949955178268487368462704L, M_PI_2l); + TEST_c_c (casinh, minus_zero, 1.5L, -0.9624236501192068949955178268487368462704L, M_PI_2l); + + TEST_c_c (casinh, -1.5L, plus_zero, -1.194763217287109304111930828519090523536L, plus_zero); + TEST_c_c (casinh, -1.5L, minus_zero, -1.194763217287109304111930828519090523536L, minus_zero); + TEST_c_c (casinh, -1.0L, plus_zero, -0.8813735870195430252326093249797923090282L, plus_zero); + TEST_c_c (casinh, -1.0L, minus_zero, -0.8813735870195430252326093249797923090282L, minus_zero); + TEST_c_c (casinh, -0.5L, plus_zero, -0.4812118250596034474977589134243684231352L, plus_zero); + TEST_c_c (casinh, -0.5L, minus_zero, -0.4812118250596034474977589134243684231352L, minus_zero); + TEST_c_c (casinh, 0.5L, plus_zero, 0.4812118250596034474977589134243684231352L, plus_zero); + TEST_c_c (casinh, 0.5L, minus_zero, 0.4812118250596034474977589134243684231352L, minus_zero); + TEST_c_c (casinh, 1.0L, plus_zero, 0.8813735870195430252326093249797923090282L, plus_zero); + TEST_c_c (casinh, 1.0L, minus_zero, 0.8813735870195430252326093249797923090282L, minus_zero); + TEST_c_c (casinh, 1.5L, plus_zero, 1.194763217287109304111930828519090523536L, plus_zero); + TEST_c_c (casinh, 1.5L, minus_zero, 1.194763217287109304111930828519090523536L, minus_zero); + TEST_c_c (casinh, 0.75L, 1.25L, 1.03171853444778027336364058631006594L, 0.911738290968487636358489564316731207L); TEST_c_c (casinh, -2, -3, -1.9686379257930962917886650952454982L, -0.96465850440760279204541105949953237L); @@ -2554,6 +2658,24 @@ csqrt_test (void) part). */ TEST_c_c (csqrt, 0, -1, M_SQRT_2_2, -M_SQRT_2_2); + TEST_c_c (csqrt, 0x1.fffffep+127L, 0x1.fffffep+127L, 2.026714405498316804978751017492482558075e+19L, 8.394925938143272988211878516208015586281e+18L); + TEST_c_c (csqrt, 0x1.fffffep+127L, 1.0L, 1.844674352395372953599975585936590505260e+19L, 2.710505511993121390769065968615872097053e-20L); + TEST_c_c (csqrt, 0x1p-149L, 0x1p-149L, 4.112805464342778798097003462770175200803e-23L, 1.703579802732953750368659735601389709551e-23L); + TEST_c_c (csqrt, 0x1p-147L, 0x1p-147L, 8.225610928685557596194006925540350401606e-23L, 3.407159605465907500737319471202779419102e-23L); + +#ifndef TEST_FLOAT + TEST_c_c (csqrt, 0x1.fffffffffffffp+1023L, 0x1.fffffffffffffp+1023L, 1.473094556905565378990473658199034571917e+154L, 6.101757441282702188537080005372547713595e+153L); + TEST_c_c (csqrt, 0x1.fffffffffffffp+1023L, 0x1p+1023L, 1.379778091031440685006200821918878702861e+154L, 3.257214233483129514781233066898042490248e+153L); + TEST_c_c (csqrt, 0x1p-1074L, 0x1p-1074L, 2.442109726130830256743814843868934877597e-162L, 1.011554969366634726113090867589031782487e-162L); + TEST_c_c (csqrt, 0x1p-1073L, 0x1p-1073L, 3.453664695497464982856905711457966660085e-162L, 1.430554756764195530630723976279903095110e-162L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (csqrt, 0x1.fp+16383L, 0x1.fp+16383L, 1.179514222452201722651836720466795901016e+2466L, 4.885707879516577666702435054303191575148e+2465L); + TEST_c_c (csqrt, 0x1.fp+16383L, 0x1p+16383L, 1.106698967236475180613254276996359485630e+2466L, 2.687568007603946993388538156299100955642e+2465L); + TEST_c_c (csqrt, 0x1p-16440L, 0x1p-16441L, 3.514690655930285351254618340783294558136e-2475L, 8.297059146828716918029689466551384219370e-2476L); +#endif + END (csqrt, complex); } @@ -3454,6 +3576,19 @@ hypot_test (void) TEST_ff_f (hypot, 0.75L, 1.25L, 1.45773797371132511771853821938639577L); +#if !(defined TEST_FLOAT && defined TEST_INLINE) + TEST_ff_f (hypot, 0x3p125L, 0x4p125L, 0x5p125L); + TEST_ff_f (hypot, 0x1.234566p-126L, 0x1.234566p-126L, 1.891441686191081936598531534017449451173e-38L); +#endif + +#ifndef TEST_FLOAT + TEST_ff_f (hypot, 0x3p1021L, 0x4p1021L, 0x5p1021L); +#endif + +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 && !defined TEST_INLINE + TEST_ff_f (hypot, 0x3p16381L, 0x4p16381L, 0x5p16381L); +#endif + END (hypot); } diff --git a/math/math_private.h b/math/math_private.h deleted file mode 100644 index 777762dd33..0000000000 --- a/math/math_private.h +++ /dev/null @@ -1,392 +0,0 @@ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * from: @(#)fdlibm.h 5.1 93/09/24 - */ - -#ifndef _MATH_PRIVATE_H_ -#define _MATH_PRIVATE_H_ - -#include <endian.h> -#include <stdint.h> -#include <sys/types.h> - -/* The original fdlibm code used statements like: - n0 = ((*(int*)&one)>>29)^1; * index of high word * - ix0 = *(n0+(int*)&x); * high word of x * - ix1 = *((1-n0)+(int*)&x); * low word of x * - to dig two 32 bit words out of the 64 bit IEEE floating point - value. That is non-ANSI, and, moreover, the gcc instruction - scheduler gets it wrong. We instead use the following macros. - Unlike the original code, we determine the endianness at compile - time, not at run time; I don't see much benefit to selecting - endianness at run time. */ - -/* A union which permits us to convert between a double and two 32 bit - ints. */ - -#if __FLOAT_WORD_ORDER == BIG_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t msw; - u_int32_t lsw; - } parts; - uint64_t word; -} ieee_double_shape_type; - -#endif - -#if __FLOAT_WORD_ORDER == LITTLE_ENDIAN - -typedef union -{ - double value; - struct - { - u_int32_t lsw; - u_int32_t msw; - } parts; - uint64_t word; -} ieee_double_shape_type; - -#endif - -/* Get two 32 bit ints from a double. */ - -#define EXTRACT_WORDS(ix0,ix1,d) \ -do { \ - ieee_double_shape_type ew_u; \ - ew_u.value = (d); \ - (ix0) = ew_u.parts.msw; \ - (ix1) = ew_u.parts.lsw; \ -} while (0) - -/* Get the more significant 32 bit int from a double. */ - -#define GET_HIGH_WORD(i,d) \ -do { \ - ieee_double_shape_type gh_u; \ - gh_u.value = (d); \ - (i) = gh_u.parts.msw; \ -} while (0) - -/* Get the less significant 32 bit int from a double. */ - -#define GET_LOW_WORD(i,d) \ -do { \ - ieee_double_shape_type gl_u; \ - gl_u.value = (d); \ - (i) = gl_u.parts.lsw; \ -} while (0) - -/* Get all in one, efficient on 64-bit machines. */ -#define EXTRACT_WORDS64(i,d) \ -do { \ - ieee_double_shape_type gh_u; \ - gh_u.value = (d); \ - (i) = gh_u.word; \ -} while (0) - -/* Set a double from two 32 bit ints. */ - -#define INSERT_WORDS(d,ix0,ix1) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.parts.msw = (ix0); \ - iw_u.parts.lsw = (ix1); \ - (d) = iw_u.value; \ -} while (0) - -/* Get all in one, efficient on 64-bit machines. */ -#define INSERT_WORDS64(d,i) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.word = (i); \ - (d) = iw_u.value; \ -} while (0) - -/* Set the more significant 32 bits of a double from an int. */ - -#define SET_HIGH_WORD(d,v) \ -do { \ - ieee_double_shape_type sh_u; \ - sh_u.value = (d); \ - sh_u.parts.msw = (v); \ - (d) = sh_u.value; \ -} while (0) - -/* Set the less significant 32 bits of a double from an int. */ - -#define SET_LOW_WORD(d,v) \ -do { \ - ieee_double_shape_type sl_u; \ - sl_u.value = (d); \ - sl_u.parts.lsw = (v); \ - (d) = sl_u.value; \ -} while (0) - -/* A union which permits us to convert between a float and a 32 bit - int. */ - -typedef union -{ - float value; - u_int32_t word; -} ieee_float_shape_type; - -/* Get a 32 bit int from a float. */ - -#define GET_FLOAT_WORD(i,d) \ -do { \ - ieee_float_shape_type gf_u; \ - gf_u.value = (d); \ - (i) = gf_u.word; \ -} while (0) - -/* Set a float from a 32 bit int. */ - -#define SET_FLOAT_WORD(d,i) \ -do { \ - ieee_float_shape_type sf_u; \ - sf_u.word = (i); \ - (d) = sf_u.value; \ -} while (0) - -/* Get long double macros from a separate header. */ -#include <math_ldbl.h> - -/* ieee style elementary functions */ -extern double __ieee754_sqrt (double); -extern double __ieee754_acos (double); -extern double __ieee754_acosh (double); -extern double __ieee754_log (double); -extern double __ieee754_atanh (double); -extern double __ieee754_asin (double); -extern double __ieee754_atan2 (double,double); -extern double __ieee754_exp (double); -extern double __ieee754_exp2 (double); -extern double __ieee754_exp10 (double); -extern double __ieee754_cosh (double); -extern double __ieee754_fmod (double,double); -extern double __ieee754_pow (double,double); -extern double __ieee754_lgamma_r (double,int *); -extern double __ieee754_gamma_r (double,int *); -extern double __ieee754_lgamma (double); -extern double __ieee754_gamma (double); -extern double __ieee754_log10 (double); -extern double __ieee754_log2 (double); -extern double __ieee754_sinh (double); -extern double __ieee754_hypot (double,double); -extern double __ieee754_j0 (double); -extern double __ieee754_j1 (double); -extern double __ieee754_y0 (double); -extern double __ieee754_y1 (double); -extern double __ieee754_jn (int,double); -extern double __ieee754_yn (int,double); -extern double __ieee754_remainder (double,double); -extern int32_t __ieee754_rem_pio2 (double,double*); -extern double __ieee754_scalb (double,double); - -/* fdlibm kernel function */ -extern double __kernel_standard (double,double,int); -extern float __kernel_standard_f (float,float,int); -extern double __kernel_sin (double,double,int); -extern double __kernel_cos (double,double); -extern double __kernel_tan (double,double,int); -extern int __kernel_rem_pio2 (double*,double*,int,int,int, const int32_t*); - -/* internal functions. */ -extern double __copysign (double x, double __y); - -extern inline double __copysign (double x, double y) -{ return __builtin_copysign (x, y); } - -/* ieee style elementary float functions */ -extern float __ieee754_sqrtf (float); -extern float __ieee754_acosf (float); -extern float __ieee754_acoshf (float); -extern float __ieee754_logf (float); -extern float __ieee754_atanhf (float); -extern float __ieee754_asinf (float); -extern float __ieee754_atan2f (float,float); -extern float __ieee754_expf (float); -extern float __ieee754_exp2f (float); -extern float __ieee754_exp10f (float); -extern float __ieee754_coshf (float); -extern float __ieee754_fmodf (float,float); -extern float __ieee754_powf (float,float); -extern float __ieee754_lgammaf_r (float,int *); -extern float __ieee754_gammaf_r (float,int *); -extern float __ieee754_lgammaf (float); -extern float __ieee754_gammaf (float); -extern float __ieee754_log10f (float); -extern float __ieee754_log2f (float); -extern float __ieee754_sinhf (float); -extern float __ieee754_hypotf (float,float); -extern float __ieee754_j0f (float); -extern float __ieee754_j1f (float); -extern float __ieee754_y0f (float); -extern float __ieee754_y1f (float); -extern float __ieee754_jnf (int,float); -extern float __ieee754_ynf (int,float); -extern float __ieee754_remainderf (float,float); -extern int32_t __ieee754_rem_pio2f (float,float*); -extern float __ieee754_scalbf (float,float); - - -/* float versions of fdlibm kernel functions */ -extern float __kernel_sinf (float,float,int); -extern float __kernel_cosf (float,float); -extern float __kernel_tanf (float,float,int); -extern int __kernel_rem_pio2f (float*,float*,int,int,int, const int32_t*); - -/* internal functions. */ -extern float __copysignf (float x, float __y); - -extern inline float __copysignf (float x, float y) -{ return __builtin_copysignf (x, y); } - -/* ieee style elementary long double functions */ -extern long double __ieee754_sqrtl (long double); -extern long double __ieee754_acosl (long double); -extern long double __ieee754_acoshl (long double); -extern long double __ieee754_logl (long double); -extern long double __ieee754_atanhl (long double); -extern long double __ieee754_asinl (long double); -extern long double __ieee754_atan2l (long double,long double); -extern long double __ieee754_expl (long double); -extern long double __ieee754_exp2l (long double); -extern long double __ieee754_exp10l (long double); -extern long double __ieee754_coshl (long double); -extern long double __ieee754_fmodl (long double,long double); -extern long double __ieee754_powl (long double,long double); -extern long double __ieee754_lgammal_r (long double,int *); -extern long double __ieee754_gammal_r (long double,int *); -extern long double __ieee754_lgammal (long double); -extern long double __ieee754_gammal (long double); -extern long double __ieee754_log10l (long double); -extern long double __ieee754_log2l (long double); -extern long double __ieee754_sinhl (long double); -extern long double __ieee754_hypotl (long double,long double); -extern long double __ieee754_j0l (long double); -extern long double __ieee754_j1l (long double); -extern long double __ieee754_y0l (long double); -extern long double __ieee754_y1l (long double); -extern long double __ieee754_jnl (int,long double); -extern long double __ieee754_ynl (int,long double); -extern long double __ieee754_remainderl (long double,long double); -extern int __ieee754_rem_pio2l (long double,long double*); -extern long double __ieee754_scalbl (long double,long double); - -/* long double versions of fdlibm kernel functions */ -extern long double __kernel_sinl (long double,long double,int); -extern long double __kernel_cosl (long double,long double); -extern long double __kernel_tanl (long double,long double,int); -extern void __kernel_sincosl (long double,long double, - long double *,long double *, int); -extern int __kernel_rem_pio2l (long double*,long double*,int,int, - int,const int*); - -#ifndef NO_LONG_DOUBLE -/* prototypes required to compile the ldbl-96 support without warnings */ -extern int __finitel (long double); -extern int __ilogbl (long double); -extern int __isinfl (long double); -extern int __isnanl (long double); -extern long double __atanl (long double); -extern long double __copysignl (long double, long double); -extern long double __expm1l (long double); -extern long double __floorl (long double); -extern long double __frexpl (long double, int *); -extern long double __ldexpl (long double, int); -extern long double __log1pl (long double); -extern long double __nanl (const char *); -extern long double __rintl (long double); -extern long double __scalbnl (long double, int); -extern long double __sqrtl (long double x); -extern long double fabsl (long double x); -extern void __sincosl (long double, long double *, long double *); -extern long double __logbl (long double x); -extern long double __significandl (long double x); - -extern inline long double __copysignl (long double x, long double y) -{ return __builtin_copysignl (x, y); } - -#endif - -/* Prototypes for functions of the IBM Accurate Mathematical Library. */ -extern double __exp1 (double __x, double __xx, double __error); -extern double __sin (double __x); -extern double __cos (double __x); -extern int __branred (double __x, double *__a, double *__aa); -extern void __doasin (double __x, double __dx, double __v[]); -extern void __dubsin (double __x, double __dx, double __v[]); -extern void __dubcos (double __x, double __dx, double __v[]); -extern double __halfulp (double __x, double __y); -extern double __sin32 (double __x, double __res, double __res1); -extern double __cos32 (double __x, double __res, double __res1); -extern double __mpsin (double __x, double __dx); -extern double __mpcos (double __x, double __dx); -extern double __mpsin1 (double __x); -extern double __mpcos1 (double __x); -extern double __slowexp (double __x); -extern double __slowpow (double __x, double __y, double __z); -extern void __docos (double __x, double __dx, double __v[]); - -#ifndef math_opt_barrier -# define math_opt_barrier(x) \ -({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; }) -# define math_force_eval(x) \ -({ __typeof (x) __x = (x); __asm __volatile__ ("" : : "m" (__x)); }) -#endif - - -/* The standards only specify one variant of the fenv.h interfaces. - But at least for some architectures we can be more efficient if we - know what operations are going to be performed. Therefore we - define additional interfaces. By default they refer to the normal - interfaces. */ -#define libc_feholdexcept(e) (void) feholdexcept (e) -#define libc_feholdexceptf(e) (void) feholdexcept (e) -#define libc_feholdexceptl(e) (void) feholdexcept (e) - -#define libc_feholdexcept_setround(e, r) \ - do { feholdexcept (e); fesetround (r); } while (0) -#define libc_feholdexcept_setroundf(e, r) \ - do { feholdexcept (e); fesetround (r); } while (0) -#define libc_feholdexcept_setroundl(e, r) \ - do { feholdexcept (e); fesetround (r); } while (0) - -#define libc_fetestexcept(e) fetestexcept (e) -#define libc_fetestexceptf(e) fetestexcept (e) -#define libc_fetestexceptl(e) fetestexcept (e) - -#define libc_fesetenv(e) (void) fesetenv (e) -#define libc_fesetenvf(e) (void) fesetenv (e) -#define libc_fesetenvl(e) (void) fesetenv (e) - -#define libc_feupdateenv(e) (void) feupdateenv (e) -#define libc_feupdateenvf(e) (void) feupdateenv (e) -#define libc_feupdateenvl(e) (void) feupdateenv (e) - -#define __nan(str) \ - (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) -#define __nanf(str) \ - (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) -#define __nanl(str) \ - (__builtin_constant_p (str) && str[0] == '\0' ? NAN : __nan (str)) - -#endif /* _MATH_PRIVATE_H_ */ diff --git a/math/s_cacosh.c b/math/s_cacosh.c index 95bf005cd3..ef49088c3d 100644 --- a/math/s_cacosh.c +++ b/math/s_cacosh.c @@ -1,5 +1,5 @@ /* Return arc hyperbole cosine for double value. - Copyright (C) 1997, 2006, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -66,9 +66,13 @@ __cacosh (__complex__ double x) } /* The factor 16 is just a guess. */ else if (16.0 * fabs (__imag__ x) < fabs (__real__ x)) - /* Kahan's formula which avoid cancellation through subtraction in - some cases. */ - res = 2.0 * __clog (__csqrt ((x + 1.0) / 2.0) + __csqrt ((x - 1.0) / 2.0)); + { + /* Kahan's formula which avoid cancellation through subtraction in + some cases. */ + res = 2.0 * __clog (__csqrt ((x + 1.0) / 2.0) + __csqrt ((x - 1.0) / 2.0)); + if (signbit (__real__ res)) + __real__ res = 0.0; + } else { __complex__ double y; diff --git a/math/s_cacoshf.c b/math/s_cacoshf.c index d634a1b756..fc716ac519 100644 --- a/math/s_cacoshf.c +++ b/math/s_cacoshf.c @@ -1,5 +1,5 @@ /* Return arc hyperbole cosine for float value. - Copyright (C) 1997, 2006, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -66,10 +66,14 @@ __cacoshf (__complex__ float x) } /* The factor 16 is just a guess. */ else if (16.0 * fabsf (__imag__ x) < fabsf (__real__ x)) - /* Kahan's formula which avoid cancellation through subtraction in - some cases. */ - res = 2.0 * __clogf (__csqrtf ((x + 1.0) / 2.0) - + __csqrtf ((x - 1.0) / 2.0)); + { + /* Kahan's formula which avoid cancellation through subtraction in + some cases. */ + res = 2.0 * __clogf (__csqrtf ((x + 1.0) / 2.0) + + __csqrtf ((x - 1.0) / 2.0)); + if (signbit (__real__ res)) + __real__ res = 0.0f; + } else { __complex__ float y; diff --git a/math/s_cacoshl.c b/math/s_cacoshl.c index dcdb1d5a64..3bcab1d277 100644 --- a/math/s_cacoshl.c +++ b/math/s_cacoshl.c @@ -1,5 +1,5 @@ /* Return arc hyperbole cosine for long double value. - Copyright (C) 1997, 1998, 2006, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -66,10 +66,14 @@ __cacoshl (__complex__ long double x) } /* The factor 16 is just a guess. */ else if (16.0L * fabsl (__imag__ x) < fabsl (__real__ x)) - /* Kahan's formula which avoid cancellation through subtraction in - some cases. */ - res = 2.0L * __clogl (__csqrtl ((x + 1.0L) / 2.0L) - + __csqrtl ((x - 1.0L) / 2.0L)); + { + /* Kahan's formula which avoid cancellation through subtraction in + some cases. */ + res = 2.0L * __clogl (__csqrtl ((x + 1.0L) / 2.0L) + + __csqrtl ((x - 1.0L) / 2.0L)); + if (signbit (__real__ res)) + __real__ res = 0.0L; + } else { __complex__ long double y; diff --git a/math/s_casinh.c b/math/s_casinh.c index e06be1cdc6..24ae676a8f 100644 --- a/math/s_casinh.c +++ b/math/s_casinh.c @@ -1,5 +1,5 @@ /* Return arc hyperbole sine for double value. - Copyright (C) 1997, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -73,6 +73,11 @@ __casinh (__complex__ double x) __imag__ y += __imag__ x; res = __clog (y); + + /* Ensure zeros have correct sign and results are correct if + very close to branch cuts. */ + __real__ res = __copysign (__real__ res, __real__ x); + __imag__ res = __copysign (__imag__ res, __imag__ x); } return res; diff --git a/math/s_casinhf.c b/math/s_casinhf.c index 37c67401a3..c7098ef309 100644 --- a/math/s_casinhf.c +++ b/math/s_casinhf.c @@ -1,5 +1,5 @@ /* Return arc hyperbole sine for float value. - Copyright (C) 1997, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -73,6 +73,11 @@ __casinhf (__complex__ float x) __imag__ y += __imag__ x; res = __clogf (y); + + /* Ensure zeros have correct sign and results are correct if + very close to branch cuts. */ + __real__ res = __copysignf (__real__ res, __real__ x); + __imag__ res = __copysignf (__imag__ res, __imag__ x); } return res; diff --git a/math/s_casinhl.c b/math/s_casinhl.c index de41cb8999..80d7bf500d 100644 --- a/math/s_casinhl.c +++ b/math/s_casinhl.c @@ -1,5 +1,5 @@ /* Return arc hyperbole sine for long double value. - Copyright (C) 1997, 1998, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -73,6 +73,11 @@ __casinhl (__complex__ long double x) __imag__ y += __imag__ x; res = __clogl (y); + + /* Ensure zeros have correct sign and results are correct if + very close to branch cuts. */ + __real__ res = __copysignl (__real__ res, __real__ x); + __imag__ res = __copysignl (__imag__ res, __imag__ x); } return res; diff --git a/math/s_csqrt.c b/math/s_csqrt.c index 76585e889c..002ea5fdc2 100644 --- a/math/s_csqrt.c +++ b/math/s_csqrt.c @@ -1,5 +1,5 @@ /* Complex square root of double value. - Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,7 +21,7 @@ #include <complex.h> #include <math.h> #include <math_private.h> - +#include <float.h> __complex__ double __csqrt (__complex__ double x) @@ -83,6 +83,22 @@ __csqrt (__complex__ double x) else { double d, r, s; + int scale = 0; + + if (fabs (__real__ x) > DBL_MAX / 2.0 + || fabs (__imag__ x) > DBL_MAX / 2.0) + { + scale = 1; + __real__ x = __scalbn (__real__ x, -2 * scale); + __imag__ x = __scalbn (__imag__ x, -2 * scale); + } + else if (fabs (__real__ x) < DBL_MIN + && fabs (__imag__ x) < DBL_MIN) + { + scale = -(DBL_MANT_DIG / 2); + __real__ x = __scalbn (__real__ x, -2 * scale); + __imag__ x = __scalbn (__imag__ x, -2 * scale); + } d = __ieee754_hypot (__real__ x, __imag__ x); /* Use the identity 2 Re res Im res = Im x @@ -98,6 +114,12 @@ __csqrt (__complex__ double x) r = fabs ((0.5 * __imag__ x) / s); } + if (scale) + { + r = __scalbn (r, scale); + s = __scalbn (s, scale); + } + __real__ res = r; __imag__ res = __copysign (s, __imag__ x); } diff --git a/math/s_csqrtf.c b/math/s_csqrtf.c index d9949c685b..6539ba2249 100644 --- a/math/s_csqrtf.c +++ b/math/s_csqrtf.c @@ -1,5 +1,5 @@ /* Complex square root of float value. - Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,7 +21,7 @@ #include <complex.h> #include <math.h> #include <math_private.h> - +#include <float.h> __complex__ float __csqrtf (__complex__ float x) @@ -83,6 +83,22 @@ __csqrtf (__complex__ float x) else { float d, r, s; + int scale = 0; + + if (fabsf (__real__ x) > FLT_MAX / 2.0f + || fabsf (__imag__ x) > FLT_MAX / 2.0f) + { + scale = 1; + __real__ x = __scalbnf (__real__ x, -2 * scale); + __imag__ x = __scalbnf (__imag__ x, -2 * scale); + } + else if (fabsf (__real__ x) < FLT_MIN + && fabsf (__imag__ x) < FLT_MIN) + { + scale = -(FLT_MANT_DIG / 2); + __real__ x = __scalbnf (__real__ x, -2 * scale); + __imag__ x = __scalbnf (__imag__ x, -2 * scale); + } d = __ieee754_hypotf (__real__ x, __imag__ x); /* Use the identity 2 Re res Im res = Im x @@ -98,6 +114,12 @@ __csqrtf (__complex__ float x) r = fabsf ((0.5f * __imag__ x) / s); } + if (scale) + { + r = __scalbnf (r, scale); + s = __scalbnf (s, scale); + } + __real__ res = r; __imag__ res = __copysignf (s, __imag__ x); } diff --git a/math/s_csqrtl.c b/math/s_csqrtl.c index 0c624c7a73..64332f67b2 100644 --- a/math/s_csqrtl.c +++ b/math/s_csqrtl.c @@ -1,5 +1,5 @@ /* Complex square root of long double value. - Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc. + Copyright (C) 1997-2012 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,7 +21,7 @@ #include <complex.h> #include <math.h> #include <math_private.h> - +#include <float.h> __complex__ long double __csqrtl (__complex__ long double x) @@ -83,6 +83,22 @@ __csqrtl (__complex__ long double x) else { long double d, r, s; + int scale = 0; + + if (fabsl (__real__ x) > LDBL_MAX / 2.0L + || fabsl (__imag__ x) > LDBL_MAX / 2.0L) + { + scale = 1; + __real__ x = __scalbnl (__real__ x, -2 * scale); + __imag__ x = __scalbnl (__imag__ x, -2 * scale); + } + else if (fabsl (__real__ x) < LDBL_MIN + && fabsl (__imag__ x) < LDBL_MIN) + { + scale = -(LDBL_MANT_DIG / 2); + __real__ x = __scalbnl (__real__ x, -2 * scale); + __imag__ x = __scalbnl (__imag__ x, -2 * scale); + } d = __ieee754_hypotl (__real__ x, __imag__ x); /* Use the identity 2 Re res Im res = Im x @@ -98,6 +114,12 @@ __csqrtl (__complex__ long double x) r = fabsl ((0.5L * __imag__ x) / s); } + if (scale) + { + r = __scalbnl (r, scale); + s = __scalbnl (s, scale); + } + __real__ res = r; __imag__ res = __copysignl (s, __imag__ x); } |