diff options
Diffstat (limited to 'sysdeps/libm-ieee754')
60 files changed, 3641 insertions, 242 deletions
diff --git a/sysdeps/libm-ieee754/s_cacos.c b/sysdeps/libm-ieee754/s_cacos.c new file mode 100644 index 0000000000..9b007598f0 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacos.c @@ -0,0 +1,41 @@ +/* Return cosine of complex double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +__complex__ double +__cacos (__complex__ double x) +{ + __complex__ double y; + __complex__ double res; + + y = __casin (x); + + __real__ res = (double) M_PI_2 - __real__ y; + __imag__ res = -__imag__ y; + + return res; +} +weak_alias (__cacos, cacos) +#ifdef NO_LONG_DOUBLE +strong_alias (__cacos, __cacosl) +weak_alias (__cacos, cacosl) +#endif diff --git a/sysdeps/libm-ieee754/s_cacosf.c b/sysdeps/libm-ieee754/s_cacosf.c new file mode 100644 index 0000000000..6fb132dc23 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacosf.c @@ -0,0 +1,37 @@ +/* Return cosine of complex float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +__complex__ float +__cacosf (__complex__ float x) +{ + __complex__ float y; + __complex__ float res; + + y = __casinf (x); + + __real__ res = (float) M_PI_2 - __real__ y; + __imag__ res = -__imag__ y; + + return res; +} +weak_alias (__cacosf, cacosf) diff --git a/sysdeps/libm-ieee754/s_cacosh.c b/sysdeps/libm-ieee754/s_cacosh.c new file mode 100644 index 0000000000..d938c64473 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacosh.c @@ -0,0 +1,88 @@ +/* Return arc hyperbole cosine for double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__cacosh (__complex__ double x) +{ + __complex__ double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VAL; + + if (rcls == FP_NAN) + __imag__ res = __nan (""); + else + __imag__ res = __copysign ((rcls == FP_INFINITE + ? (__real__ x < 0.0 + ? M_PI - M_PI_4 : M_PI_4) + : M_PI_2), __imag__ x); + } + else if (rcls == FP_INFINITE) + { + __real__ res = HUGE_VAL; + + if (icls >= FP_ZERO) + __imag__ res = __copysign (signbit (__real__ x) ? M_PI : 0.0, + __imag__ x); + else + __imag__ res = __nan (""); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + __real__ res = 0.0; + __imag__ res = __copysign (M_PI_2, __imag__ x); + } + else + { + __complex__ double y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) - 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrt (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clog (y); + } + + return res; +} +weak_alias (__cacosh, cacosh) +#ifdef NO_LONG_DOUBLE +strong_alias (__cacosh, __cacoshl) +weak_alias (__cacosh, cacoshl) +#endif diff --git a/sysdeps/libm-ieee754/s_cacoshf.c b/sysdeps/libm-ieee754/s_cacoshf.c new file mode 100644 index 0000000000..bcfebea123 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacoshf.c @@ -0,0 +1,84 @@ +/* Return arc hyperbole cosine for float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__cacoshf (__complex__ float x) +{ + __complex__ float res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VALF; + + if (rcls == FP_NAN) + __imag__ res = __nanf (""); + else + __imag__ res = __copysignf ((rcls == FP_INFINITE + ? (__real__ x < 0.0 + ? M_PI - M_PI_4 : M_PI_4) + : M_PI_2), __imag__ x); + } + else if (rcls == FP_INFINITE) + { + __real__ res = HUGE_VALF; + + if (icls >= FP_ZERO) + __imag__ res = __copysignf (signbit (__real__ x) ? M_PI : 0.0, + __imag__ x); + else + __imag__ res = __nanf (""); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + __real__ res = 0.0; + __imag__ res = __copysignf (M_PI_2, __imag__ x); + } + else + { + __complex__ float y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) - 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrtf (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clogf (y); + } + + return res; +} +weak_alias (__cacoshf, cacoshf) diff --git a/sysdeps/libm-ieee754/s_cacoshl.c b/sysdeps/libm-ieee754/s_cacoshl.c new file mode 100644 index 0000000000..ed5980c551 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacoshl.c @@ -0,0 +1,84 @@ +/* Return arc hyperbole cosine for long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__cacoshl (__complex__ long double x) +{ + __complex__ long double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VALL; + + if (rcls == FP_NAN) + __imag__ res = __nanl (""); + else + __imag__ res = __copysignl ((rcls == FP_INFINITE + ? (__real__ x < 0.0 + ? M_PI - M_PI_4 : M_PI_4) + : M_PI_2), __imag__ x); + } + else if (rcls == FP_INFINITE) + { + __real__ res = HUGE_VALL; + + if (icls >= FP_ZERO) + __imag__ res = __copysignl (signbit (__real__ x) ? M_PI : 0.0, + __imag__ x); + else + __imag__ res = __nanl (""); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + __real__ res = 0.0; + __imag__ res = __copysignl (M_PI_2, __imag__ x); + } + else + { + __complex__ long double y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) - 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrtl (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clogl (y); + } + + return res; +} +weak_alias (__cacoshl, cacoshl) diff --git a/sysdeps/libm-ieee754/s_cacosl.c b/sysdeps/libm-ieee754/s_cacosl.c new file mode 100644 index 0000000000..7a60de5a50 --- /dev/null +++ b/sysdeps/libm-ieee754/s_cacosl.c @@ -0,0 +1,37 @@ +/* Return cosine of complex long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +__complex__ long double +__cacosl (__complex__ long double x) +{ + __complex__ long double y; + __complex__ long double res; + + y = __casinl (x); + + __real__ res = M_PI_2 - __real__ y; + __imag__ res = -__imag__ y; + + return res; +} +weak_alias (__cacosl, cacosl) diff --git a/sysdeps/libm-ieee754/s_casin.c b/sysdeps/libm-ieee754/s_casin.c new file mode 100644 index 0000000000..516aea0464 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casin.c @@ -0,0 +1,66 @@ +/* Return arc sine of complex double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__casin (__complex__ double x) +{ + __complex__ double res; + + if (isnan (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0) + { + res = x; + } + else if (__isinf (__real__ x) || __isinf (__imag__ x)) + { + __real__ res = __nan (""); + __imag__ res = __copysign (HUGE_VAL, __imag__ x); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + __complex__ double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __casinh (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__casin, casin) +#ifdef NO_LONG_DOUBLE +strong_alias (__casin, __casinl) +weak_alias (__casin, casinl) +#endif diff --git a/sysdeps/libm-ieee754/s_casinf.c b/sysdeps/libm-ieee754/s_casinf.c new file mode 100644 index 0000000000..aaf0d66c28 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casinf.c @@ -0,0 +1,62 @@ +/* Return arc sine of complex float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__casinf (__complex__ float x) +{ + __complex__ float res; + + if (isnan (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0) + { + res = x; + } + else if (__isinff (__real__ x) || __isinff (__imag__ x)) + { + __real__ res = __nanf (""); + __imag__ res = __copysignf (HUGE_VALF, __imag__ x); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + __complex__ float y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __casinhf (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__casinf, casinf) diff --git a/sysdeps/libm-ieee754/s_casinh.c b/sysdeps/libm-ieee754/s_casinh.c new file mode 100644 index 0000000000..da7d1ed429 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casinh.c @@ -0,0 +1,84 @@ +/* Return arc hyperbole sine for double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__casinh (__complex__ double x) +{ + __complex__ double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysign (HUGE_VAL, __real__ x); + + if (rcls == FP_NAN) + __imag__ res = __nan (""); + else + __imag__ res = __copysign (rcls >= FP_ZERO ? M_PI_2 : M_PI_4, + __imag__ x); + } + else if (rcls <= FP_INFINITE) + { + __real__ res = __real__ x; + if ((rcls == FP_INFINITE && icls >= FP_ZERO) + || (rcls == FP_NAN && icls == FP_ZERO)) + __imag__ res = __copysign (0.0, __imag__ x); + else + __imag__ res = __nan (""); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + __complex__ double y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) + 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrt (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clog (y); + } + + return res; +} +weak_alias (__casinh, casinh) +#ifdef NO_LONG_DOUBLE +strong_alias (__casinh, __casinhl) +weak_alias (__casinh, casinhl) +#endif diff --git a/sysdeps/libm-ieee754/s_casinhf.c b/sysdeps/libm-ieee754/s_casinhf.c new file mode 100644 index 0000000000..e8441f4fa0 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casinhf.c @@ -0,0 +1,80 @@ +/* Return arc hyperbole sine for float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__casinhf (__complex__ float x) +{ + __complex__ float res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysignf (HUGE_VALF, __real__ x); + + if (rcls == FP_NAN) + __imag__ res = __nanf (""); + else + __imag__ res = __copysignf (rcls >= FP_ZERO ? M_PI_2 : M_PI_4, + __imag__ x); + } + else if (rcls <= FP_INFINITE) + { + __real__ res = __real__ x; + if ((rcls == FP_INFINITE && icls >= FP_ZERO) + || (rcls == FP_NAN && icls == FP_ZERO)) + __imag__ res = __copysignf (0.0, __imag__ x); + else + __imag__ res = __nanf (""); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + __complex__ float y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) + 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrtf (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clogf (y); + } + + return res; +} +weak_alias (__casinhf, casinhf) diff --git a/sysdeps/libm-ieee754/s_casinhl.c b/sysdeps/libm-ieee754/s_casinhl.c new file mode 100644 index 0000000000..cc6757b189 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casinhl.c @@ -0,0 +1,80 @@ +/* Return arc hyperbole sine for long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__casinhl (__complex__ long double x) +{ + __complex__ long double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysignl (HUGE_VALL, __real__ x); + + if (rcls == FP_NAN) + __imag__ res = __nanl (""); + else + __imag__ res = __copysignl (rcls >= FP_ZERO ? M_PI_2 : M_PI_4, + __imag__ x); + } + else if (rcls <= FP_INFINITE) + { + __real__ res = __real__ x; + if ((rcls == FP_INFINITE && icls >= FP_ZERO) + || (rcls == FP_NAN && icls == FP_ZERO)) + __imag__ res = __copysignl (0.0, __imag__ x); + else + __imag__ res = __nanl (""); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + __complex__ long double y; + + __real__ y = (__real__ x - __imag__ x) * (__real__ x + __imag__ x) + 1.0; + __imag__ y = 2.0 * __real__ x * __imag__ x; + + y = __csqrtl (y); + + __real__ y += __real__ x; + __imag__ y += __imag__ x; + + res = __clogl (y); + } + + return res; +} +weak_alias (__casinhl, casinhl) diff --git a/sysdeps/libm-ieee754/s_casinl.c b/sysdeps/libm-ieee754/s_casinl.c new file mode 100644 index 0000000000..cc750d4a11 --- /dev/null +++ b/sysdeps/libm-ieee754/s_casinl.c @@ -0,0 +1,62 @@ +/* Return arc sine of complex long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__casinl (__complex__ long double x) +{ + __complex__ long double res; + + if (isnan (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0) + { + res = x; + } + else if (__isinfl (__real__ x) || __isinfl (__imag__ x)) + { + __real__ res = __nanl (""); + __imag__ res = __copysignl (HUGE_VALL, __imag__ x); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + __complex__ long double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __casinhl (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__casinl, casinl) diff --git a/sysdeps/libm-ieee754/s_catan.c b/sysdeps/libm-ieee754/s_catan.c new file mode 100644 index 0000000000..bab87e95a8 --- /dev/null +++ b/sysdeps/libm-ieee754/s_catan.c @@ -0,0 +1,89 @@ +/* Return arc tangent of complex double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ double +__catan (__complex__ double x) +{ + __complex__ double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (rcls == FP_INFINITE) + { + __real__ res = __copysign (M_PI_2, __real__ x); + __imag__ res = __copysign (0.0, __imag__ x); + } + else if (icls == FP_INFINITE) + { + if (rcls >= FP_ZERO) + __real__ res = __copysign (M_PI_2, __real__ x); + else + __real__ res = __nan (""); + __imag__ res = __copysign (0.0, __imag__ x); + } + else if (icls == FP_ZERO || icls == FP_INFINITE) + { + __real__ res = __nan (""); + __imag__ res = __copysign (0.0, __imag__ x); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + double r2, num, den; + + r2 = __real__ x * __real__ x; + + den = 1 - r2 - __imag__ x * __imag__ x; + + __real__ res = 0.5 * __atan ((2.0 * __real__ x) / den); + + num = __imag__ x + 1.0; + num = r2 + num * num; + + den = __imag__ x - 1.0; + den = r2 + den * den; + + __imag__ res = 0.25 * __ieee754_log (num / den); + } + + return res; +} +weak_alias (__catan, catan) +#ifdef NO_LONG_DOUBLE +strong_alias (__catan, __catanl) +weak_alias (__catan, catanl) +#endif diff --git a/sysdeps/libm-ieee754/s_catanf.c b/sysdeps/libm-ieee754/s_catanf.c new file mode 100644 index 0000000000..92bdac9ec8 --- /dev/null +++ b/sysdeps/libm-ieee754/s_catanf.c @@ -0,0 +1,85 @@ +/* Return arc tangent of complex float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ float +__catanf (__complex__ float x) +{ + __complex__ float res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (rcls == FP_INFINITE) + { + __real__ res = __copysignf (M_PI_2, __real__ x); + __imag__ res = __copysignf (0.0, __imag__ x); + } + else if (icls == FP_INFINITE) + { + if (rcls >= FP_ZERO) + __real__ res = __copysignf (M_PI_2, __real__ x); + else + __real__ res = __nanf (""); + __imag__ res = __copysignf (0.0, __imag__ x); + } + else if (icls == FP_ZERO || icls == FP_INFINITE) + { + __real__ res = __nanf (""); + __imag__ res = __copysignf (0.0, __imag__ x); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + float r2, num, den; + + r2 = __real__ x * __real__ x; + + den = 1 - r2 - __imag__ x * __imag__ x; + + __real__ res = 0.5 * __atanf ((2.0 * __real__ x) / den); + + num = __imag__ x + 1.0; + num = r2 + num * num; + + den = __imag__ x - 1.0; + den = r2 + den * den; + + __imag__ res = 0.25 * __ieee754_logf (num / den); + } + + return res; +} +weak_alias (__catanf, catanf) diff --git a/sysdeps/libm-ieee754/s_catanh.c b/sysdeps/libm-ieee754/s_catanh.c new file mode 100644 index 0000000000..6c4b10e3db --- /dev/null +++ b/sysdeps/libm-ieee754/s_catanh.c @@ -0,0 +1,84 @@ +/* Return arc hyperbole tangent for double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ double +__catanh (__complex__ double x) +{ + __complex__ double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysign (0.0, __real__ x); + __imag__ res = __copysign (M_PI_2, __imag__ x); + } + else if (rcls == FP_INFINITE || rcls == FP_ZERO) + { + __real__ res = __copysign (0.0, __real__ x); + if (icls >= FP_ZERO) + __imag__ res = __copysign (M_PI_2, __imag__ x); + else + __imag__ res = __nan (""); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + double i2, num, den; + + i2 = __imag__ x * __imag__ x; + + num = 1.0 - __real__ x; + num = i2 + num * num; + + den = 1.0 + __real__ x; + den = i2 + den * den; + + __real__ res = 0.25 * __ieee754_log (num / den); + + den = 1 - __real__ x * __real__ x - i2; + + __imag__ res = 0.5 * __atan ((2.0 * __imag__ x) / den); + } + + return res; +} +weak_alias (__catanh, catanh) +#ifdef NO_LONG_DOUBLE +strong_alias (__catanh, __catanhl) +weak_alias (__catanh, catanhl) +#endif diff --git a/sysdeps/libm-ieee754/s_catanhf.c b/sysdeps/libm-ieee754/s_catanhf.c new file mode 100644 index 0000000000..5d195be905 --- /dev/null +++ b/sysdeps/libm-ieee754/s_catanhf.c @@ -0,0 +1,80 @@ +/* Return arc hyperbole tangent for float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ float +__catanhf (__complex__ float x) +{ + __complex__ float res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysignf (0.0, __real__ x); + __imag__ res = __copysignf (M_PI_2, __imag__ x); + } + else if (rcls == FP_INFINITE || rcls == FP_ZERO) + { + __real__ res = __copysignf (0.0, __real__ x); + if (icls >= FP_ZERO) + __imag__ res = __copysignf (M_PI_2, __imag__ x); + else + __imag__ res = __nanf (""); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + float i2, num, den; + + i2 = __imag__ x * __imag__ x; + + num = 1.0 - __real__ x; + num = i2 + num * num; + + den = 1.0 + __real__ x; + den = i2 + den * den; + + __real__ res = 0.25 * __ieee754_logf (num / den); + + den = 1 - __real__ x * __real__ x - i2; + + __imag__ res = 0.5 * __atanf ((2.0 * __imag__ x) / den); + } + + return res; +} +weak_alias (__catanhf, catanhf) diff --git a/sysdeps/libm-ieee754/s_catanhl.c b/sysdeps/libm-ieee754/s_catanhl.c new file mode 100644 index 0000000000..d8396a7961 --- /dev/null +++ b/sysdeps/libm-ieee754/s_catanhl.c @@ -0,0 +1,80 @@ +/* Return arc hyperbole tangent for long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ long double +__catanhl (__complex__ long double x) +{ + __complex__ long double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = __copysignl (0.0, __real__ x); + __imag__ res = __copysignl (M_PI_2, __imag__ x); + } + else if (rcls == FP_INFINITE || rcls == FP_ZERO) + { + __real__ res = __copysignl (0.0, __real__ x); + if (icls >= FP_ZERO) + __imag__ res = __copysignl (M_PI_2, __imag__ x); + else + __imag__ res = __nanl (""); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + long double i2, num, den; + + i2 = __imag__ x * __imag__ x; + + num = 1.0 - __real__ x; + num = i2 + num * num; + + den = 1.0 + __real__ x; + den = i2 + den * den; + + __real__ res = 0.25 * __ieee754_logl (num / den); + + den = 1 - __real__ x * __real__ x - i2; + + __imag__ res = 0.5 * __atanl ((2.0 * __imag__ x) / den); + } + + return res; +} +weak_alias (__catanhl, catanhl) diff --git a/sysdeps/libm-ieee754/s_catanl.c b/sysdeps/libm-ieee754/s_catanl.c new file mode 100644 index 0000000000..2fd8a4fa08 --- /dev/null +++ b/sysdeps/libm-ieee754/s_catanl.c @@ -0,0 +1,85 @@ +/* Return arc tangent of complex long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ long double +__catanl (__complex__ long double x) +{ + __complex__ long double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (rcls == FP_INFINITE) + { + __real__ res = __copysignl (M_PI_2, __real__ x); + __imag__ res = __copysignl (0.0, __imag__ x); + } + else if (icls == FP_INFINITE) + { + if (rcls >= FP_ZERO) + __real__ res = __copysignl (M_PI_2, __real__ x); + else + __real__ res = __nanl (""); + __imag__ res = __copysignl (0.0, __imag__ x); + } + else if (icls == FP_ZERO || icls == FP_INFINITE) + { + __real__ res = __nanl (""); + __imag__ res = __copysignl (0.0, __imag__ x); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else if (rcls == FP_ZERO && icls == FP_ZERO) + { + res = x; + } + else + { + long double r2, num, den; + + r2 = __real__ x * __real__ x; + + den = 1 - r2 - __imag__ x * __imag__ x; + + __real__ res = 0.5 * __atanl ((2.0 * __real__ x) / den); + + num = __imag__ x + 1.0; + num = r2 + num * num; + + den = __imag__ x - 1.0; + den = r2 + den * den; + + __imag__ res = 0.25 * __ieee754_logl (num / den); + } + + return res; +} +weak_alias (__catanl, catanl) diff --git a/sysdeps/libm-ieee754/s_ccos.c b/sysdeps/libm-ieee754/s_ccos.c new file mode 100644 index 0000000000..8a4b55dd99 --- /dev/null +++ b/sysdeps/libm-ieee754/s_ccos.c @@ -0,0 +1,64 @@ +/* Return cosine of complex double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__ccos (__complex__ double x) +{ + __complex__ double res; + + if (!isfinite (__real__ x) || __isnan (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nan (""); + __imag__ res = 0.0; + } + else if (__isinf (__imag__ x)) + { + __real__ res = HUGE_VAL; + __imag__ res = __nan (""); + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + __complex__ double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + res = __ccosh (y); + } + + return res; +} +weak_alias (__ccos, ccos) +#ifdef NO_LONG_DOUBLE +strong_alias (__ccos, __ccosl) +weak_alias (__ccos, ccosl) +#endif diff --git a/sysdeps/libm-ieee754/s_ccosf.c b/sysdeps/libm-ieee754/s_ccosf.c new file mode 100644 index 0000000000..9d1a97239c --- /dev/null +++ b/sysdeps/libm-ieee754/s_ccosf.c @@ -0,0 +1,60 @@ +/* Return cosine of complex float value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__ccosf (__complex__ float x) +{ + __complex__ float res; + + if (!isfinite (__real__ x) || __isnanf (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nanf (""); + __imag__ res = 0.0; + } + else if (__isinff (__imag__ x)) + { + __real__ res = HUGE_VALF; + __imag__ res = __nanf (""); + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + __complex__ float y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + res = __ccoshf (y); + } + + return res; +} +weak_alias (__ccosf, ccosf) diff --git a/sysdeps/libm-ieee754/s_ccosh.c b/sysdeps/libm-ieee754/s_ccosh.c index f01b245e77..b9d7b82f5e 100644 --- a/sysdeps/libm-ieee754/s_ccosh.c +++ b/sysdeps/libm-ieee754/s_ccosh.c @@ -21,47 +21,47 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ double __ccosh (__complex__ double x) { __complex__ double retval; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); __real__ x = fabs (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - double exp_val = __exp (__real__ x); - double rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + double cosh_val = __ieee754_cosh (__real__ x); - __real__ retval = 0.5 * (exp_val + rec_exp_val) * __cos (__imag__ x); - __imag__ retval = 0.5 * (exp_val + rec_exp_val) * __sin (__imag__ x); + __real__ retval = cosh_val * __cos (__imag__ x); + __imag__ retval = cosh_val * __sin (__imag__ x); } else { - if (__real__ x == 0) - { - __imag__ retval = 0.0; - __real__ retval = __nan ("") + __nan (""); - } - else - { - __real__ retval = __nan (""); - __imag__ retval = __nan ("") + __nan (""); - } + __imag__ retval = __real__ x == 0.0 ? 0.0 : __nan (""); + __real__ retval = __nan ("") + __nan (""); } } - else if (__isinf (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = HUGE_VAL; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysign (HUGE_VAL, __cos (__imag__ x)); __imag__ retval = __copysign (HUGE_VAL, __sin (__imag__ x)); } @@ -74,16 +74,8 @@ __ccosh (__complex__ double x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nan (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nan (""); - __imag__ retval = __nan (""); - } + __real__ retval = __nan (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nan (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_ccoshf.c b/sysdeps/libm-ieee754/s_ccoshf.c index 9f2774b6c7..758ec5b579 100644 --- a/sysdeps/libm-ieee754/s_ccoshf.c +++ b/sysdeps/libm-ieee754/s_ccoshf.c @@ -21,49 +21,47 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ float __ccoshf (__complex__ float x) { __complex__ float retval; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); __real__ x = fabsf (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - float exp_val = __expf (__real__ x); - float rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + float cosh_val = __ieee754_coshf (__real__ x); - __real__ retval = (0.5 * (exp_val + rec_exp_val) - * __cosf (__imag__ x)); - __imag__ retval = (0.5 * (exp_val + rec_exp_val) - * __sinf (__imag__ x)); + __real__ retval = cosh_val * __cosf (__imag__ x); + __imag__ retval = cosh_val * __sinf (__imag__ x); } else { - if (__real__ x == 0) - { - __imag__ retval = 0.0; - __real__ retval = __nanf ("") + __nanf (""); - } - else - { - __real__ retval = __nanf (""); - __imag__ retval = __nanf ("") + __nanf (""); - } + __imag__ retval = __real__ x == 0.0 ? 0.0 : __nanf (""); + __real__ retval = __nanf ("") + __nanf (""); } } - else if (__isinff (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = HUGE_VALF; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysignf (HUGE_VALF, __cosf (__imag__ x)); __imag__ retval = __copysignf (HUGE_VALF, __sinf (__imag__ x)); } @@ -76,16 +74,8 @@ __ccoshf (__complex__ float x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nanf (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nanf (""); - __imag__ retval = __nanf (""); - } + __real__ retval = __nanf (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nanf (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_ccoshl.c b/sysdeps/libm-ieee754/s_ccoshl.c index fd553418b4..5e8c399fb2 100644 --- a/sysdeps/libm-ieee754/s_ccoshl.c +++ b/sysdeps/libm-ieee754/s_ccoshl.c @@ -21,49 +21,47 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ long double __ccoshl (__complex__ long double x) { __complex__ long double retval; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); __real__ x = fabsl (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - long double exp_val = __expl (__real__ x); - long double rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + long double cosh_val = __ieee754_coshl (__real__ x); - __real__ retval = (0.5 * (exp_val + rec_exp_val) - * __cosl (__imag__ x)); - __imag__ retval = (0.5 * (exp_val + rec_exp_val) - * __sinl (__imag__ x)); + __real__ retval = cosh_val * __cosl (__imag__ x); + __imag__ retval = cosh_val * __sinl (__imag__ x); } else { - if (__real__ x == 0) - { - __imag__ retval = 0.0; - __real__ retval = __nanl ("") + __nanl (""); - } - else - { - __real__ retval = __nanl (""); - __imag__ retval = __nanl ("") + __nanl (""); - } + __imag__ retval = __real__ x == 0.0 ? 0.0 : __nanl (""); + __real__ retval = __nanl ("") + __nanl (""); } } - else if (__isinfl (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = HUGE_VALL; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysignl (HUGE_VALL, __cosl (__imag__ x)); __imag__ retval = __copysignl (HUGE_VALL, __sinl (__imag__ x)); } @@ -76,16 +74,8 @@ __ccoshl (__complex__ long double x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nanl (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nanl (""); - __imag__ retval = __nanl (""); - } + __real__ retval = __nanl (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nanl (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_ccosl.c b/sysdeps/libm-ieee754/s_ccosl.c new file mode 100644 index 0000000000..a41d48b970 --- /dev/null +++ b/sysdeps/libm-ieee754/s_ccosl.c @@ -0,0 +1,60 @@ +/* Return cosine of complex long double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__ccosl (__complex__ long double x) +{ + __complex__ long double res; + + if (!isfinite (__real__ x) || __isnanl (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nanl (""); + __imag__ res = 0.0; + } + else if (__isinfl (__imag__ x)) + { + __real__ res = HUGE_VALL; + __imag__ res = __nanl (""); + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + __complex__ long double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + res = __ccoshl (y); + } + + return res; +} +weak_alias (__ccosl, ccosl) diff --git a/sysdeps/libm-ieee754/s_ceill.c b/sysdeps/libm-ieee754/s_ceill.c index c5c86487ea..773be32995 100644 --- a/sysdeps/libm-ieee754/s_ceill.c +++ b/sysdeps/libm-ieee754/s_ceill.c @@ -48,34 +48,34 @@ static long double huge = 1.0e4930; GET_LDOUBLE_WORDS(se,i0,i1,x); sx = (se>>15)&1; j0 = (se&0x7fff)-0x3fff; - if(j0<32) { + if(j0<31) { if(j0<0) { /* raise inexact if x != 0 */ if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ if(sx) {es=0x8000;i0=0;i1=0;} else if((i0|i1)!=0) { es=0x3fff;i0=0;i1=0;} } } else { - i = (0xffffffff)>>j0; + i = (0x7fffffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(sx==0) { - if (j0>0) i0 += (0x80000000)>>(j0-1); + if (j0>0) i0 += (0x80000000)>>j0; else ++se; } i0 &= (~i); i1=0; } } - } else if (j0>63) { + } else if (j0>62) { if(j0==0x4000) return x+x; /* inf or NaN */ else return x; /* x is integral */ } else { - i = ((u_int32_t)(0xffffffff))>>(j0-32); + i = ((u_int32_t)(0xffffffff))>>(j0-31); if((i1&i)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(sx==0) { - if(j0==32) i0+=1; + if(j0==31) i0+=1; else { - j = i1 + (1<<(64-j0)); + j = i1 + (1<<(63-j0)); if(j<i1) i0+=1; /* got a carry */ i1 = j; } diff --git a/sysdeps/libm-ieee754/s_cexpf.c b/sysdeps/libm-ieee754/s_cexpf.c index c5d8f0cc07..0d4372103b 100644 --- a/sysdeps/libm-ieee754/s_cexpf.c +++ b/sysdeps/libm-ieee754/s_cexpf.c @@ -21,17 +21,23 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ float __cexpf (__complex__ float x) { __complex__ float retval; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - float exp_val = __expf (__real__ x); + /* Imaginary part is finite. */ + float exp_val = __ieee754_expf (__real__ x); if (isfinite (exp_val)) { @@ -52,14 +58,17 @@ __cexpf (__complex__ float x) __imag__ retval = __nanf (""); } } - else if (__isinff (__real__ x)) + else if (rcls == FP_INFINITE) { - if (isfinite (__imag__ x)) + /* Real part is infinite. */ + if (icls >= FP_ZERO) { + /* Imaginary part is finite. */ float value = signbit (__real__ x) ? 0.0 : HUGE_VALF; - if (__imag__ x == 0.0) + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = value; __imag__ retval = __imag__ x; } @@ -82,7 +91,7 @@ __cexpf (__complex__ float x) } else { - /* If the real part is NaN the result is NaN + iNan. */ + /* If the real part is NaN the result is NaN + iNaN. */ __real__ retval = __nanf (""); __imag__ retval = __nanf (""); } diff --git a/sysdeps/libm-ieee754/s_cexpl.c b/sysdeps/libm-ieee754/s_cexpl.c index f1cdf43ec8..ac27e1eeb8 100644 --- a/sysdeps/libm-ieee754/s_cexpl.c +++ b/sysdeps/libm-ieee754/s_cexpl.c @@ -1,4 +1,4 @@ -/* Return value of complex exponential function for float complex value. +/* Return value of complex exponential function for long double complex value. Copyright (C) 1997 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,17 +21,23 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ long double __cexpl (__complex__ long double x) { __complex__ long double retval; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - long double exp_val = __expl (__real__ x); + /* Imaginary part is finite. */ + long double exp_val = __ieee754_expl (__real__ x); if (isfinite (exp_val)) { @@ -52,14 +58,17 @@ __cexpl (__complex__ long double x) __imag__ retval = __nanl (""); } } - else if (__isinfl (__real__ x)) + else if (rcls == FP_INFINITE) { - if (isfinite (__imag__ x)) + /* Real part is infinite. */ + if (icls >= FP_ZERO) { + /* Imaginary part is finite. */ long double value = signbit (__real__ x) ? 0.0 : HUGE_VALL; - if (__imag__ x == 0.0) + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = value; __imag__ retval = __imag__ x; } @@ -89,4 +98,4 @@ __cexpl (__complex__ long double x) return retval; } -weak_alias (__cexp, cexp) +weak_alias (__cexpl, cexpl) diff --git a/sysdeps/libm-ieee754/s_clog.c b/sysdeps/libm-ieee754/s_clog.c index f00753b3bb..c14a734759 100644 --- a/sysdeps/libm-ieee754/s_clog.c +++ b/sysdeps/libm-ieee754/s_clog.c @@ -28,17 +28,20 @@ __complex__ double __clog (__complex__ double x) { __complex__ double result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - if (__real__ x == 0.0 && __imag__ x == 0.0) + if (rcls == FP_ZERO && icls == FP_ZERO) { + /* Real and imaginary part are 0.0. */ __imag__ result = signbit (__real__ x) ? M_PI : 0.0; - if (signbit (__imag__ x)) - __imag__ result = __copysign (__imag__ result, -1.0); + __imag__ result = __copysign (__imag__ result, __imag__ x); /* Yes, the following line raises an exception. */ __real__ result = -1.0 / fabs (__real__ x); } - else if (!__isnan (__real__ x) && !__isnan (__imag__ x)) + else if (rcls != FP_NAN && icls != FP_NAN) { + /* Neither real nor imaginary part is NaN. */ __real__ result = __ieee754_log (__ieee754_hypot (__real__ x, __imag__ x)); __imag__ result = __ieee754_atan2 (__imag__ x, __real__ x); @@ -46,7 +49,8 @@ __clog (__complex__ double x) else { __imag__ result = __nan (""); - if (__isinf (__real__ x) || __isinf (__imag__ x)) + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ __real__ result = HUGE_VAL; else __real__ result = __nan (""); diff --git a/sysdeps/libm-ieee754/s_clogf.c b/sysdeps/libm-ieee754/s_clogf.c index 4eafc82bf0..9c9aa83e33 100644 --- a/sysdeps/libm-ieee754/s_clogf.c +++ b/sysdeps/libm-ieee754/s_clogf.c @@ -28,17 +28,20 @@ __complex__ float __clogf (__complex__ float x) { __complex__ float result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - if (__real__ x == 0.0 && __imag__ x == 0.0) + if (rcls == FP_ZERO && icls == FP_ZERO) { + /* Real and imaginary part are 0.0. */ __imag__ result = signbit (__real__ x) ? M_PI : 0.0; - if (signbit (__imag__ x)) - __imag__ result = __copysignf (__imag__ result, -1.0); + __imag__ result = __copysignf (__imag__ result, __imag__ x); /* Yes, the following line raises an exception. */ __real__ result = -1.0 / fabsf (__real__ x); } - else if (!__isnanf (__real__ x) && !__isnanf (__imag__ x)) + else if (rcls != FP_NAN && icls != FP_NAN) { + /* Neither real nor imaginary part is NaN. */ __real__ result = __ieee754_logf (__ieee754_hypotf (__real__ x, __imag__ x)); __imag__ result = __ieee754_atan2f (__imag__ x, __real__ x); @@ -46,7 +49,8 @@ __clogf (__complex__ float x) else { __imag__ result = __nanf (""); - if (__isinff (__real__ x) || __isinff (__imag__ x)) + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ __real__ result = HUGE_VALF; else __real__ result = __nanf (""); diff --git a/sysdeps/libm-ieee754/s_clogl.c b/sysdeps/libm-ieee754/s_clogl.c index a299a95c03..51bee372a6 100644 --- a/sysdeps/libm-ieee754/s_clogl.c +++ b/sysdeps/libm-ieee754/s_clogl.c @@ -28,17 +28,20 @@ __complex__ long double __clogl (__complex__ long double x) { __complex__ long double result; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - if (__real__ x == 0.0 && __imag__ x == 0.0) + if (rcls == FP_ZERO && icls == FP_ZERO) { + /* Real and imaginary part are 0.0. */ __imag__ result = signbit (__real__ x) ? M_PI : 0.0; - if (signbit (__imag__ x)) - __imag__ result = __copysignl (__imag__ result, -1.0); + __imag__ result = __copysignl (__imag__ result, __imag__ x); /* Yes, the following line raises an exception. */ __real__ result = -1.0 / fabsl (__real__ x); } - else if (!__isnanl (__real__ x) && !__isnanl (__imag__ x)) + else if (rcls != FP_NAN && icls != FP_NAN) { + /* Neither real nor imaginary part is NaN. */ __real__ result = __ieee754_logl (__ieee754_hypotl (__real__ x, __imag__ x)); __imag__ result = __ieee754_atan2l (__imag__ x, __real__ x); @@ -46,7 +49,8 @@ __clogl (__complex__ long double x) else { __imag__ result = __nanl (""); - if (__isinfl (__real__ x) || __isinfl (__imag__ x)) + if (rcls == FP_INFINITE || icls == FP_INFINITE) + /* Real or imaginary part is infinite. */ __real__ result = HUGE_VALL; else __real__ result = __nanl (""); diff --git a/sysdeps/libm-ieee754/s_cpow.c b/sysdeps/libm-ieee754/s_cpow.c new file mode 100644 index 0000000000..074b38bd2d --- /dev/null +++ b/sysdeps/libm-ieee754/s_cpow.c @@ -0,0 +1,34 @@ +/* Complex power of double values. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__cpow (__complex__ double x, __complex__ double c) +{ + return __cexp (c * __clog (x)); +} +weak_alias (__cpow, cpow) +#ifdef NO_LONG_DOUBLE +strong_alias (__cpow, __cpowl) +weak_alias (__cpow, cpowl) +#endif diff --git a/sysdeps/libm-ieee754/s_cpowf.c b/sysdeps/libm-ieee754/s_cpowf.c new file mode 100644 index 0000000000..fa4541ca2b --- /dev/null +++ b/sysdeps/libm-ieee754/s_cpowf.c @@ -0,0 +1,30 @@ +/* Complex power of float values. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__cpowf (__complex__ float x, __complex__ float c) +{ + return __cexpf (c * __clogf (x)); +} +weak_alias (__cpowf, cpowf) diff --git a/sysdeps/libm-ieee754/s_cpowl.c b/sysdeps/libm-ieee754/s_cpowl.c new file mode 100644 index 0000000000..69097d5dcc --- /dev/null +++ b/sysdeps/libm-ieee754/s_cpowl.c @@ -0,0 +1,30 @@ +/* Complex power of long double values. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__cpowl (__complex__ long double x, __complex__ long double c) +{ + return __cexpl (c * __clogl (x)); +} +weak_alias (__cpowl, cpowl) diff --git a/sysdeps/libm-ieee754/s_csin.c b/sysdeps/libm-ieee754/s_csin.c new file mode 100644 index 0000000000..4639bcaaa6 --- /dev/null +++ b/sysdeps/libm-ieee754/s_csin.c @@ -0,0 +1,67 @@ +/* Complex sine function for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ double +__csin (__complex__ double x) +{ + __complex__ double res; + + if (!isfinite (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nan (""); + __imag__ res = 0.0; + } + else if (__isinf (__imag__ x)) + { + __real__ res = __nan (""); + __imag__ res = __imag__ x; + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + __complex__ double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __csinh (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__csin, csin) +#ifdef NO_LONG_DOUBLE +strong_alias (__csin, __csinl) +weak_alias (__csin, csinl) +#endif diff --git a/sysdeps/libm-ieee754/s_csinf.c b/sysdeps/libm-ieee754/s_csinf.c new file mode 100644 index 0000000000..f7f10e6b6f --- /dev/null +++ b/sysdeps/libm-ieee754/s_csinf.c @@ -0,0 +1,63 @@ +/* Complex sine function for float. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ float +__csinf (__complex__ float x) +{ + __complex__ float res; + + if (!isfinite (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nanf (""); + __imag__ res = 0.0; + } + else if (__isinff (__imag__ x)) + { + __real__ res = __nanf (""); + __imag__ res = __imag__ x; + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + __complex__ float y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __csinhf (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__csinf, csinf) diff --git a/sysdeps/libm-ieee754/s_csinh.c b/sysdeps/libm-ieee754/s_csinh.c index aab041bdf9..05cec85e7c 100644 --- a/sysdeps/libm-ieee754/s_csinh.c +++ b/sysdeps/libm-ieee754/s_csinh.c @@ -21,32 +21,38 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ double __csinh (__complex__ double x) { __complex__ double retval; int negate = signbit (__real__ x); + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); __real__ x = fabs (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - double exp_val = __exp (__real__ x); - double rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + double sinh_val = __ieee754_sinh (__real__ x); - __real__ retval = 0.5 * (exp_val - rec_exp_val) * __cos (__imag__ x); - __imag__ retval = 0.5 * (exp_val - rec_exp_val) * __sin (__imag__ x); + __real__ retval = sinh_val * __cos (__imag__ x); + __imag__ retval = sinh_val * __sin (__imag__ x); if (negate) __real__ retval = -__real__ retval; } else { - if (__real__ x == 0) + if (rcls == FP_ZERO) { + /* Real part is 0.0. */ __real__ retval = __copysign (0.0, negate ? -1.0 : 1.0); __imag__ retval = __nan ("") + __nan (""); } @@ -57,15 +63,18 @@ __csinh (__complex__ double x) } } } - else if (__isinf (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = negate ? -HUGE_VAL : HUGE_VAL; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysign (HUGE_VAL, __cos (__imag__ x)); __imag__ retval = __copysign (HUGE_VAL, __sin (__imag__ x)); @@ -81,16 +90,8 @@ __csinh (__complex__ double x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nan (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nan (""); - __imag__ retval = __nan (""); - } + __real__ retval = __nan (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nan (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_csinhf.c b/sysdeps/libm-ieee754/s_csinhf.c index 204bbfebb9..93f83cf7b6 100644 --- a/sysdeps/libm-ieee754/s_csinhf.c +++ b/sysdeps/libm-ieee754/s_csinhf.c @@ -21,34 +21,38 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ float __csinhf (__complex__ float x) { __complex__ float retval; int negate = signbit (__real__ x); + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); __real__ x = fabsf (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - float exp_val = __expf (__real__ x); - float rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + float sinh_val = __ieee754_sinhf (__real__ x); - __real__ retval = (0.5 * (exp_val - rec_exp_val) - * __cosf (__imag__ x)); - __imag__ retval = (0.5 * (exp_val - rec_exp_val) - * __sinf (__imag__ x)); + __real__ retval = sinh_val * __cosf (__imag__ x); + __imag__ retval = sinh_val * __sinf (__imag__ x); if (negate) __real__ retval = -__real__ retval; } else { - if (__real__ x == 0) + if (rcls == FP_ZERO) { + /* Real part is 0.0. */ __real__ retval = __copysignf (0.0, negate ? -1.0 : 1.0); __imag__ retval = __nanf ("") + __nanf (""); } @@ -59,15 +63,18 @@ __csinhf (__complex__ float x) } } } - else if (__isinff (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = negate ? -HUGE_VALF : HUGE_VALF; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysignf (HUGE_VALF, __cosf (__imag__ x)); __imag__ retval = __copysignf (HUGE_VALF, __sinf (__imag__ x)); @@ -83,16 +90,8 @@ __csinhf (__complex__ float x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nanf (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nanf (""); - __imag__ retval = __nanf (""); - } + __real__ retval = __nanf (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nanf (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_csinhl.c b/sysdeps/libm-ieee754/s_csinhl.c index e403dd4796..8388a40b46 100644 --- a/sysdeps/libm-ieee754/s_csinhl.c +++ b/sysdeps/libm-ieee754/s_csinhl.c @@ -21,34 +21,38 @@ #include <complex.h> #include <math.h> +#include "math_private.h" + __complex__ long double __csinhl (__complex__ long double x) { __complex__ long double retval; int negate = signbit (__real__ x); + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); - __real__ x = fabs (__real__ x); + __real__ x = fabsl (__real__ x); - if (isfinite (__real__ x)) + if (rcls >= FP_ZERO) { - if (isfinite (__imag__ x)) + /* Real part is finite. */ + if (icls >= FP_ZERO) { - long double exp_val = __expl (__real__ x); - long double rec_exp_val = 1.0 / exp_val; + /* Imaginary part is finite. */ + long double sinh_val = __ieee754_sinhl (__real__ x); - __real__ retval = (0.5 * (exp_val - rec_exp_val) - * __cosl (__imag__ x)); - __imag__ retval = (0.5 * (exp_val - rec_exp_val) - * __sinl (__imag__ x)); + __real__ retval = sinh_val * __cosl (__imag__ x); + __imag__ retval = sinh_val * __sinl (__imag__ x); if (negate) __real__ retval = -__real__ retval; } else { - if (__real__ x == 0) + if (rcls == FP_ZERO) { + /* Real part is 0.0. */ __real__ retval = __copysignl (0.0, negate ? -1.0 : 1.0); __imag__ retval = __nanl ("") + __nanl (""); } @@ -59,15 +63,18 @@ __csinhl (__complex__ long double x) } } } - else if (__isinfl (__real__ x)) + else if (rcls == FP_INFINITE) { - if (__imag__ x == 0.0) + /* Real part is infinite. */ + if (icls == FP_ZERO) { + /* Imaginary part is 0.0. */ __real__ retval = negate ? -HUGE_VALL : HUGE_VALL; __imag__ retval = __imag__ x; } - else if (isfinite (__imag__ x)) + else if (icls > FP_ZERO) { + /* Imaginary part is finite. */ __real__ retval = __copysignl (HUGE_VALL, __cosl (__imag__ x)); __imag__ retval = __copysignl (HUGE_VALL, __sinl (__imag__ x)); @@ -83,16 +90,8 @@ __csinhl (__complex__ long double x) } else { - if (__imag__ x == 0.0) - { - __real__ retval = __nanl (""); - __imag__ retval = __imag__ x; - } - else - { - __real__ retval = __nanl (""); - __imag__ retval = __nanl (""); - } + __real__ retval = __nanl (""); + __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nanl (""); } return retval; diff --git a/sysdeps/libm-ieee754/s_csinl.c b/sysdeps/libm-ieee754/s_csinl.c new file mode 100644 index 0000000000..513c144198 --- /dev/null +++ b/sysdeps/libm-ieee754/s_csinl.c @@ -0,0 +1,63 @@ +/* Complex sine function for long double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + + +__complex__ long double +__csinl (__complex__ long double x) +{ + __complex__ long double res; + + if (!isfinite (__real__ x) || isnan (__imag__ x)) + { + if (__real__ x == 0.0 || __imag__ x == 0.0) + { + __real__ res = __nanl (""); + __imag__ res = 0.0; + } + else if (__isinfl (__imag__ x)) + { + __real__ res = __nanl (""); + __imag__ res = __imag__ x; + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + __complex__ long double y; + + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __csinhl (y); + + __real__ res = __imag__ y; + __imag__ res = -__real__ y; + } + + return res; +} +weak_alias (__csinl, csinl) diff --git a/sysdeps/libm-ieee754/s_csqrt.c b/sysdeps/libm-ieee754/s_csqrt.c new file mode 100644 index 0000000000..c5c609bd8c --- /dev/null +++ b/sysdeps/libm-ieee754/s_csqrt.c @@ -0,0 +1,111 @@ +/* Complex square root of double value. + Copyright (C) 1997 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. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ double +__csqrt (__complex__ double x) +{ + __complex__ double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VAL; + __imag__ res = __imag__ x; + } + else if (rcls == FP_INFINITE) + { + if (__real__ x < 0.0) + { + __real__ res = icls == FP_NAN ? __nan ("") : 0; + __imag__ res = __copysign (HUGE_VAL, __imag__ x); + } + else + { + __real__ res = __real__ x; + __imag__ res = (icls == FP_NAN + ? __nan ("") : __copysign (0.0, __imag__ x)); + } + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + if (icls == FP_ZERO) + { + if (__real__ x < 0.0) + { + __real__ res = 0.0; + __imag__ res = __copysign (__ieee754_sqrt (-__real__ x), + __imag__ x); + } + else + { + __real__ res = fabs (__ieee754_sqrt (__real__ x)); + __imag__ res = __copysign (0.0, __imag__ x); + } + } + else if (rcls == FP_ZERO) + { + double r = __ieee754_sqrt (0.5 * fabs (__imag__ x)); + + __real__ res = __copysign (r, __imag__ x); + __imag__ res = r; + } + else + { + __complex__ double q; + double t, r; + + if (fabs (__imag__ x) < 2.0e-4 * fabs (__real__ x)) + t = 0.25 * __imag__ x * (__imag__ x / __real__ x); + else + t = 0.5 * (__ieee754_hypot (__real__ x, __imag__ x) - __real__ x); + + r = __ieee754_sqrt (t); + + __real__ q = __imag__ x / (2.0 * r); + __imag__ q = r; + + /* Heron iteration in complex arithmetic. */ + res = 0.5 * (q + q / x); + } + } + + return res; +} +weak_alias (__csqrt, csqrt) +#ifdef NO_LONG_DOUBLE +strong_alias (__csqrt, __csqrtl) +weak_alias (__csqrt, csqrtl) +#endif diff --git a/sysdeps/libm-ieee754/s_csqrtf.c b/sysdeps/libm-ieee754/s_csqrtf.c new file mode 100644 index 0000000000..2289045cfd --- /dev/null +++ b/sysdeps/libm-ieee754/s_csqrtf.c @@ -0,0 +1,107 @@ +/* Complex square root of float value. + Copyright (C) 1997 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. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ float +__csqrtf (__complex__ float x) +{ + __complex__ float res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VALF; + __imag__ res = __imag__ x; + } + else if (rcls == FP_INFINITE) + { + if (__real__ x < 0.0) + { + __real__ res = icls == FP_NAN ? __nanf ("") : 0; + __imag__ res = __copysignf (HUGE_VALF, __imag__ x); + } + else + { + __real__ res = __real__ x; + __imag__ res = (icls == FP_NAN + ? __nanf ("") : __copysignf (0.0, __imag__ x)); + } + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + if (icls == FP_ZERO) + { + if (__real__ x < 0.0) + { + __real__ res = 0.0; + __imag__ res = __copysignf (__ieee754_sqrtf (-__real__ x), + __imag__ x); + } + else + { + __real__ res = fabsf (__ieee754_sqrtf (__real__ x)); + __imag__ res = __copysignf (0.0, __imag__ x); + } + } + else if (rcls == FP_ZERO) + { + float r = __ieee754_sqrtf (0.5 * fabsf (__imag__ x)); + + __real__ res = __copysignf (r, __imag__ x); + __imag__ res = r; + } + else + { + __complex__ float q; + float t, r; + + if (fabsf (__imag__ x) < 2.0e-4 * fabsf (__real__ x)) + t = 0.25 * __imag__ x * (__imag__ x / __real__ x); + else + t = 0.5 * (__ieee754_hypotf (__real__ x, __imag__ x) - __real__ x); + + r = __ieee754_sqrtf (t); + + __real__ q = __imag__ x / (2.0 * r); + __imag__ q = r; + + /* Heron iteration in complex arithmetic. */ + res = 0.5 * (q + q / x); + } + } + + return res; +} +weak_alias (__csqrtf, csqrtf) diff --git a/sysdeps/libm-ieee754/s_csqrtl.c b/sysdeps/libm-ieee754/s_csqrtl.c new file mode 100644 index 0000000000..3de7310c73 --- /dev/null +++ b/sysdeps/libm-ieee754/s_csqrtl.c @@ -0,0 +1,107 @@ +/* Complex square root of long double value. + Copyright (C) 1997 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. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ long double +__csqrtl (__complex__ long double x) +{ + __complex__ long double res; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE) + { + if (icls == FP_INFINITE) + { + __real__ res = HUGE_VALL; + __imag__ res = __imag__ x; + } + else if (rcls == FP_INFINITE) + { + if (__real__ x < 0.0) + { + __real__ res = icls == FP_NAN ? __nanl ("") : 0; + __imag__ res = __copysignl (HUGE_VALL, __imag__ x); + } + else + { + __real__ res = __real__ x; + __imag__ res = (icls == FP_NAN + ? __nanl ("") : __copysignl (0.0, __imag__ x)); + } + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + if (icls == FP_ZERO) + { + if (__real__ x < 0.0) + { + __real__ res = 0.0; + __imag__ res = __copysignl (__ieee754_sqrtl (-__real__ x), + __imag__ x); + } + else + { + __real__ res = fabsl (__ieee754_sqrtl (__real__ x)); + __imag__ res = __copysignl (0.0, __imag__ x); + } + } + else if (rcls == FP_ZERO) + { + long double r = __ieee754_sqrtl (0.5 * fabsl (__imag__ x)); + + __real__ res = __copysignl (r, __imag__ x); + __imag__ res = r; + } + else + { + __complex__ long double q; + long double t, r; + + if (fabsl (__imag__ x) < 2.0e-4 * fabsl (__real__ x)) + t = 0.25 * __imag__ x * (__imag__ x / __real__ x); + else + t = 0.5 * (__ieee754_hypotl (__real__ x, __imag__ x) - __real__ x); + + r = __ieee754_sqrtl (t); + + __real__ q = __imag__ x / (2.0 * r); + __imag__ q = r; + + /* Heron iteration in complex arithmetic. */ + res = 0.5 * (q + q / x); + } + } + + return res; +} +weak_alias (__csqrtl, csqrtl) diff --git a/sysdeps/libm-ieee754/s_ctan.c b/sysdeps/libm-ieee754/s_ctan.c new file mode 100644 index 0000000000..f448395c7e --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctan.c @@ -0,0 +1,64 @@ +/* Complex tangent function for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ double +__ctan (__complex__ double x) +{ + __complex__ double res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinf (__imag__ x)) + { + __real__ res = __copysign (0.0, __real__ x); + __imag__ res = __copysign (1.0, __imag__ x); + } + else if (__real__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + double den = (__cos (2.0 * __real__ x) + + __ieee754_cosh (2.0 * __imag__ x)); + + __real__ res = __sin (2.0 * __real__ x) / den; + __imag__ res = __ieee754_sinh (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctan, ctan) +#ifdef NO_LONG_DOUBLE +strong_alias (__ctan, __ctanl) +weak_alias (__ctan, ctanl) +#endif diff --git a/sysdeps/libm-ieee754/s_ctanf.c b/sysdeps/libm-ieee754/s_ctanf.c new file mode 100644 index 0000000000..99011fa41d --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctanf.c @@ -0,0 +1,60 @@ +/* Complex tangent function for float. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ float +__ctanf (__complex__ float x) +{ + __complex__ float res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinff (__imag__ x)) + { + __real__ res = __copysignf (0.0, __real__ x); + __imag__ res = __copysignf (1.0, __imag__ x); + } + else if (__real__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + float den = (__cosf (2.0 * __real__ x) + + __ieee754_coshf (2.0 * __imag__ x)); + + __real__ res = __sinf (2.0 * __real__ x) / den; + __imag__ res = __ieee754_sinhf (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctanf, ctanf) diff --git a/sysdeps/libm-ieee754/s_ctanh.c b/sysdeps/libm-ieee754/s_ctanh.c new file mode 100644 index 0000000000..7c9b3197ef --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctanh.c @@ -0,0 +1,64 @@ +/* Complex hyperbole tangent for double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ double +__ctanh (__complex__ double x) +{ + __complex__ double res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinf (__real__ x)) + { + __real__ res = __copysign (1.0, __real__ x); + __imag__ res = __copysign (0.0, __imag__ x); + } + else if (__imag__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nan (""); + __imag__ res = __nan (""); + } + } + else + { + double den = (__ieee754_cosh (2.0 * __real__ x) + + __cos (2.0 * __imag__ x)); + + __real__ res = __ieee754_sinh (2.0 * __real__ x) / den; + __imag__ res = __sin (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctanh, ctanh) +#ifdef NO_LONG_DOUBLE +strong_alias (__ctanh, __ctanhl) +weak_alias (__ctanh, ctanhl) +#endif diff --git a/sysdeps/libm-ieee754/s_ctanhf.c b/sysdeps/libm-ieee754/s_ctanhf.c new file mode 100644 index 0000000000..1bdbc0fdc5 --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctanhf.c @@ -0,0 +1,60 @@ +/* Complex hyperbole tangent for float. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ float +__ctanhf (__complex__ float x) +{ + __complex__ float res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinff (__real__ x)) + { + __real__ res = __copysignf (1.0, __real__ x); + __imag__ res = __copysignf (0.0, __imag__ x); + } + else if (__imag__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanf (""); + __imag__ res = __nanf (""); + } + } + else + { + float den = (__ieee754_coshf (2.0 * __real__ x) + + __cosf (2.0 * __imag__ x)); + + __real__ res = __ieee754_sinhf (2.0 * __real__ x) / den; + __imag__ res = __sinf (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctanhf, ctanhf) diff --git a/sysdeps/libm-ieee754/s_ctanhl.c b/sysdeps/libm-ieee754/s_ctanhl.c new file mode 100644 index 0000000000..b34aeb77dd --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctanhl.c @@ -0,0 +1,60 @@ +/* Complex hyperbole tangent for long double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ long double +__ctanhl (__complex__ long double x) +{ + __complex__ long double res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinfl (__real__ x)) + { + __real__ res = __copysignl (1.0, __real__ x); + __imag__ res = __copysignl (0.0, __imag__ x); + } + else if (__imag__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + long double den = (__ieee754_coshl (2.0 * __real__ x) + + __cosl (2.0 * __imag__ x)); + + __real__ res = __ieee754_sinhl (2.0 * __real__ x) / den; + __imag__ res = __sinl (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctanhl, ctanhl) diff --git a/sysdeps/libm-ieee754/s_ctanl.c b/sysdeps/libm-ieee754/s_ctanl.c new file mode 100644 index 0000000000..82f86fc148 --- /dev/null +++ b/sysdeps/libm-ieee754/s_ctanl.c @@ -0,0 +1,60 @@ +/* Complex tangent function for long double. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <complex.h> +#include <math.h> + +#include "math_private.h" + + +__complex__ long double +__ctanl (__complex__ long double x) +{ + __complex__ double res; + + if (!finite (__real__ x) || !finite (__imag__ x)) + { + if (__isinfl (__imag__ x)) + { + __real__ res = __copysignl (0.0, __real__ x); + __imag__ res = __copysignl (1.0, __imag__ x); + } + else if (__real__ x == 0.0) + { + res = x; + } + else + { + __real__ res = __nanl (""); + __imag__ res = __nanl (""); + } + } + else + { + long double den = (__cosl (2.0 * __real__ x) + + __ieee754_coshl (2.0 * __imag__ x)); + + __real__ res = __sinl (2.0 * __real__ x) / den; + __imag__ res = __ieee754_sinhl (2.0 * __imag__ x) / den; + } + + return res; +} +weak_alias (__ctanl, ctanl) diff --git a/sysdeps/libm-ieee754/s_floorl.c b/sysdeps/libm-ieee754/s_floorl.c index 8cd81c6302..0eb0bec9b8 100644 --- a/sysdeps/libm-ieee754/s_floorl.c +++ b/sysdeps/libm-ieee754/s_floorl.c @@ -48,7 +48,7 @@ static long double huge = 1.0e4930; GET_LDOUBLE_WORDS(se,i0,i1,x); sx = (se>>15)&1; j0 = (se&0x7fff)-0x3fff; - if(j0<32) { + if(j0<31) { if(j0<0) { /* raise inexact if x != 0 */ if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ if(sx==0) {se=0;i0=i1=0;} @@ -56,26 +56,26 @@ static long double huge = 1.0e4930; { se=0xbfff;i0;i1=0;} } } else { - i = (0xffffffff)>>j0; + i = (0x7fffffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(sx) { - if (j0>0) i0 += (0x80000000)>>(j0-1); + if (j0>0) i0 += (0x80000000)>>j0; else ++se; i0 &= (~i); i1=0; } } - } else if (j0>63) { + } else if (j0>62) { if(j0==0x4000) return x+x; /* inf or NaN */ else return x; /* x is integral */ } else { - i = ((u_int32_t)(0xffffffff))>>(j0-32); + i = ((u_int32_t)(0xffffffff))>>(j0-31); if((i1&i)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(sx) { - if(j0==32) i0+=1; + if(j0==31) i0+=1; else { - j = i1+(1<<(64-j0)); + j = i1+(1<<(63-j0)); if(j<i1) i0 +=1 ; /* got a carry */ i1=j; } diff --git a/sysdeps/libm-ieee754/s_nearbyint.c b/sysdeps/libm-ieee754/s_nearbyint.c new file mode 100644 index 0000000000..32f5bf9447 --- /dev/null +++ b/sysdeps/libm-ieee754/s_nearbyint.c @@ -0,0 +1,98 @@ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ +/* + * ==================================================== + * 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. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $"; +#endif + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +TWO52[2]={ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + +#ifdef __STDC__ + double __nearbyint(double x) +#else + double __nearbyint(x) + double x; +#endif +{ + fenv_t env; + int32_t i0,j0,sx; + u_int32_t i,i1; + double w,t; + EXTRACT_WORDS(i0,i1,x); + sx = (i0>>31)&1; + j0 = ((i0>>20)&0x7ff)-0x3ff; + if(j0<20) { + if(j0<0) { + if(((i0&0x7fffffff)|i1)==0) return x; + i1 |= (i0&0x0fffff); + i0 &= 0xfffe0000; + i0 |= ((i1|-i1)>>12)&0x80000; + SET_HIGH_WORD(x,i0); + feholdexcept (&env); + w = TWO52[sx]+x; + t = w-TWO52[sx]; + fesetenv (&env); + GET_HIGH_WORD(i0,t); + SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); + return t; + } else { + i = (0x000fffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + i>>=1; + if(((i0&i)|i1)!=0) { + if(j0==19) i1 = 0x40000000; else + i0 = (i0&(~i))|((0x20000)>>j0); + } + } + } else if (j0>51) { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-20); + if((i1&i)==0) return x; /* x is integral */ + i>>=1; + if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); + } + INSERT_WORDS(x,i0,i1); + feholdexcept (&env); + w = TWO52[sx]+x; + t = w-TWO52[sx]; + fesetenv (&env); + return t; +} +weak_alias (__nearbyint, nearbyint) +#ifdef NO_LONG_DOUBLE +strong_alias (__nearbyint, __nearbyintl) +weak_alias (__nearbyint, nearbyintl) +#endif diff --git a/sysdeps/libm-ieee754/s_nearbyintf.c b/sysdeps/libm-ieee754/s_nearbyintf.c new file mode 100644 index 0000000000..dc33fa59f9 --- /dev/null +++ b/sysdeps/libm-ieee754/s_nearbyintf.c @@ -0,0 +1,80 @@ +/* s_rintf.c -- float version of s_rint.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +/* + * ==================================================== + * 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. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: s_rintf.c,v 1.4 1995/05/10 20:48:06 jtc Exp $"; +#endif + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +TWO23[2]={ + 8.3886080000e+06, /* 0x4b000000 */ + -8.3886080000e+06, /* 0xcb000000 */ +}; + +#ifdef __STDC__ + float __rintf(float x) +#else + float __rintf(x) + float x; +#endif +{ + fenv_t env; + int32_t i0,j0,sx; + u_int32_t i,i1; + float w,t; + GET_FLOAT_WORD(i0,x); + sx = (i0>>31)&1; + j0 = ((i0>>23)&0xff)-0x7f; + if(j0<23) { + if(j0<0) { + if((i0&0x7fffffff)==0) return x; + i1 = (i0&0x07fffff); + i0 &= 0xfff00000; + i0 |= ((i1|-i1)>>9)&0x400000; + SET_FLOAT_WORD(x,i0); + feholdexcept (&env); + w = TWO23[sx]+x; + t = w-TWO23[sx]; + fesetenv (&env); + GET_FLOAT_WORD(i0,t); + SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31)); + return t; + } else { + i = (0x007fffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + i>>=1; + if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0); + } + } else { + if(j0==0x80) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + SET_FLOAT_WORD(x,i0); + feholdexcept (&env); + w = TWO23[sx]+x; + t = w-TWO23[sx]; + fesetenv (&env); + return t; +} +weak_alias (__rintf, rintf) diff --git a/sysdeps/libm-ieee754/s_nearbyintl.c b/sysdeps/libm-ieee754/s_nearbyintl.c new file mode 100644 index 0000000000..b6a865443a --- /dev/null +++ b/sysdeps/libm-ieee754/s_nearbyintl.c @@ -0,0 +1,104 @@ +/* s_rintl.c -- long double version of s_rint.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +/* + * ==================================================== + * 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. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* + * rintl(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rintl(x). + */ + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +TWO63[2]={ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + +#ifdef __STDC__ + long double __rintl(long double x) +#else + long double __rintl(x) + long double x; +#endif +{ + fenv_t env; + int32_t se,j0,sx; + u_int32_t i,i0,i1; + long double w,t; + GET_LDOUBLE_WORDS(se,i0,i1,x); + sx = (se>>15)&1; + j0 = (se&0x7fff)-0x3fff; + if(j0<31) { + if(j0<0) { + if(((se&0x7fff)|i0|i1)==0) return x; + i1 |= i0; + i0 &= 0xe0000000; + i0 |= (i1|-i1)&0x80000000; + SET_LDOUBLE_MSW(x,i0); + feholdexcept (&env); + w = TWO63[sx]+x; + t = w-TWO63[sx]; + fesetenv (&env); + GET_LDOUBLE_EXP(i0,t); + SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15)); + return t; + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + i>>=1; + if(((i0&i)|i1)!=0) { + if(j0==31) i1 = 0x40000000; else + i0 = (i0&(~i))|((0x20000000)>>j0); + /* Shouldn't this be + if (j0 >= 30) i1 = 0x80000000 >> (j0 - 30); + i0 = (i0&(~i))|((0x20000000)>>j0); + If yes, this should be correct in s_rint and + s_rintf, too. -- drepper@cygnus.com */ + } + } + } else if (j0>62) { + if(j0==0x4000) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-31); + if((i1&i)==0) return x; /* x is integral */ + i>>=1; + if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31)); + } + SET_LDOUBLE_WORDS(x,se,i0,i1); + feholdexcept (&env); + w = TWO63[sx]+x; + t = w-TWO63[sx]; + fesetenv (&env); + return t; +} +weak_alias (__rintl, rintl) diff --git a/sysdeps/libm-ieee754/s_remquo.c b/sysdeps/libm-ieee754/s_remquo.c index 53f26c6d89..4103155e3f 100644 --- a/sysdeps/libm-ieee754/s_remquo.c +++ b/sysdeps/libm-ieee754/s_remquo.c @@ -49,12 +49,7 @@ __remquo (double x, double y, int *quo) return (x * y) / (x * y); if (hy <= 0x7fbfffff) - { - x = __ieee754_fmod (x, 8 * y); /* now x < 8y */ - - if (fabs (x) >= 4 * fabs (y)) - cquo += 4; - } + x = __ieee754_fmod (x, 8 * y); /* now x < 8y */ if (((hx - hy) | (lx - ly)) == 0) { @@ -66,14 +61,19 @@ __remquo (double x, double y, int *quo) y = fabs (y); cquo = 0; - if (x >= 2 * y) + if (x >= 4 * y) { x -= 4 * y; + cquo += 4; + } + if (x >= 2 * y) + { + x -= 2 * y; cquo += 2; } if (x >= y) { - x -= 2 * y; + x -= y; ++cquo; } @@ -83,24 +83,30 @@ __remquo (double x, double y, int *quo) { x -= y; if (x + x >= y) - x -= y; + { + x -= y; + ++cquo; + } } } else { double y_half = 0.5 * y; - if(x > y_half) + if (x > y_half) { x -= y; if (x >= y_half) - x -= y; + { + x -= y; + ++cquo; + } } } *quo = qs ? -cquo : cquo; - GET_HIGH_WORD (hx, x); - SET_HIGH_WORD (x, hx ^ sx); + if (sx) + x = -x; return x; } weak_alias (__remquo, remquo) diff --git a/sysdeps/libm-ieee754/s_remquof.c b/sysdeps/libm-ieee754/s_remquof.c index 0968fe650b..6fa02e47b3 100644 --- a/sysdeps/libm-ieee754/s_remquof.c +++ b/sysdeps/libm-ieee754/s_remquof.c @@ -48,12 +48,7 @@ __remquof (float x, float y, int *quo) return (x * y) / (x * y); if (hy <= 0x7dffffff) - { - x = __ieee754_fmodf (x, 8 * y); /* now x < 8y */ - - if (fabs (x) >= 4 * fabs (y)) - cquo += 4; - } + x = __ieee754_fmodf (x, 8 * y); /* now x < 8y */ if ((hx - hy) == 0) { @@ -65,14 +60,19 @@ __remquof (float x, float y, int *quo) y = fabsf (y); cquo = 0; - if (x >= 2 * y) + if (x >= 4 * y) { x -= 4 * y; + cquo += 4; + } + if (x >= 2 * y) + { + x -= 2 * y; cquo += 2; } if (x >= y) { - x -= 2 * y; + x -= y; ++cquo; } @@ -82,24 +82,30 @@ __remquof (float x, float y, int *quo) { x -= y; if (x + x >= y) - x -= y; + { + x -= y; + ++cquo; + } } } else { float y_half = 0.5 * y; - if(x > y_half) + if (x > y_half) { x -= y; if (x >= y_half) - x -= y; + { + x -= y; + ++cquo; + } } } *quo = qs ? -cquo : cquo; - GET_FLOAT_WORD (hx, x); - SET_FLOAT_WORD (x, hx ^ sx); + if (sx) + x = -x; return x; } weak_alias (__remquof, remquof) diff --git a/sysdeps/libm-ieee754/s_remquol.c b/sysdeps/libm-ieee754/s_remquol.c index 9515b218c6..9ef424901b 100644 --- a/sysdeps/libm-ieee754/s_remquol.c +++ b/sysdeps/libm-ieee754/s_remquol.c @@ -23,15 +23,15 @@ #include "math_private.h" -static const double zero = 0.0; +static const long double zero = 0.0; long double -__remquol (long double x, long double y, int *quo) +__remquol (long double x, long double p, int *quo) { int32_t ex,ep,hx,hp; u_int32_t sx,lx,lp; - int cquo; + int cquo,qs; GET_LDOUBLE_WORDS (ex, hx, lx, x); GET_LDOUBLE_WORDS (ep, hp, lp, p); @@ -49,12 +49,7 @@ __remquol (long double x, long double y, int *quo) return (x * p) / (x * p); if (ep <= 0x7ffb) - { - x = __ieee754_fmodl (x, 8 * p); /* now x < 8p */ - - if (fabsl (x) >= 4 * fabsl (p)) - cquo += 4; - } + x = __ieee754_fmodl (x, 8 * p); /* now x < 8p */ if (((ex - ep) | (hx - hp) | (lx - lp)) == 0) { @@ -66,14 +61,19 @@ __remquol (long double x, long double y, int *quo) p = fabsl (p); cquo = 0; - if (x >= 2 * p) + if (x >= 4 * p) { x -= 4 * p; + cquo += 4; + } + if (x >= 2 * p) + { + x -= 2 * p; cquo += 2; } if (x >= p) { - x -= 2 * p; + x -= p; ++cquo; } @@ -83,24 +83,30 @@ __remquol (long double x, long double y, int *quo) { x -= p; if (x + x >= p) - x -= p; + { + x -= p; + ++cquo; + } } } else { long double p_half = 0.5 * p; - if(x > p_half) + if (x > p_half) { x -= p; if (x >= p_half) - x -= p; + { + x -= p; + ++cquo; + } } } *quo = qs ? -cquo : cquo; - GET_LDOUBLE_EXP (ex, x); - SET_LDOUBLE_EXP (x, ex ^ sx); + if (sx) + x = -x; return x; } weak_alias (__remquol, remquol) diff --git a/sysdeps/libm-ieee754/s_round.c b/sysdeps/libm-ieee754/s_round.c new file mode 100644 index 0000000000..fdb17f8de8 --- /dev/null +++ b/sysdeps/libm-ieee754/s_round.c @@ -0,0 +1,97 @@ +/* Round double to integer away from zero. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const double huge = 1.0e300; + + +double +__round (double x) +{ + int32_t i0, j0; + u_int32_t i1; + + EXTRACT_WORDS (i0, i1, x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) + { + if (j0 < 0) + { + if (huge + x > 0.0) + { + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3ff00000; + i1 = 0; + } + } + else + { + u_int32_t i = 0x000fffff >> j0; + if (((i0 & i) | i1) == 0) + /* X is integral. */ + return x; + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + i0 += 0x00080000 >> j0; + i0 &= ~i; + i1 = 0; + } + } + } + else if (j0 > 51) + { + if (j0 == 0x400) + /* Inf or NaN. */ + return x + x; + else + return x; + } + else + { + u_int32_t i = 0xffffffff >> (j0 - 20); + if ((i1 & i) == 0) + /* X is integral. */ + return x; + + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (51 - j0)); + if (j < i1) + i0 += 1; + i1 = j; + } + i1 &= ~i; + } + + INSERT_WORDS (x, i0, i1); + return x; +} +weak_alias (__round, round) +#ifdef NO_LONG_DOUBLE +strong_alias (__round, __roundl) +weak_alias (__round, roundl) +#endif diff --git a/sysdeps/libm-ieee754/s_roundf.c b/sysdeps/libm-ieee754/s_roundf.c new file mode 100644 index 0000000000..5dc0e368ff --- /dev/null +++ b/sysdeps/libm-ieee754/s_roundf.c @@ -0,0 +1,73 @@ +/* Round float to integer away from zero. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const float huge = 1.0e30; + + +float +__roundf (float x) +{ + int32_t i0, j0; + + GET_FLOAT_WORD (i0, x); + j0 = ((i0 >> 23) & 0xff) - 0x7f; + if (j0 < 23) + { + if (j0 < 0) + { + if (huge + x > 0.0F) + { + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3f800000; + } + } + else + { + u_int32_t i = 0x007fffff >> j0; + if ((i0 & i) == 0) + /* X is integral. */ + return x; + if (huge + x > 0.0F) + { + /* Raise inexact if x != 0. */ + i0 += 0x00400000 >> j0; + i0 &= ~i; + } + } + } + else + { + if (j0 == 0x80) + /* Inf or NaN. */ + return x + x; + else + return x; + } + + SET_FLOAT_WORD (x, i0); + return x; +} +weak_alias (__roundf, roundf) diff --git a/sysdeps/libm-ieee754/s_roundl.c b/sysdeps/libm-ieee754/s_roundl.c new file mode 100644 index 0000000000..db87154089 --- /dev/null +++ b/sysdeps/libm-ieee754/s_roundl.c @@ -0,0 +1,100 @@ +/* Round long double to integer away from zero. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const long double huge = 1.0e4930; + + +long double +__roundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + { + if (huge + x > 0.0) + { + se &= 0x8000; + if (j0 == -1) + se |= 0x3fff; + i0 = i1 = 0; + } + } + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + /* X is integral. */ + return x; + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + u_int32_t j = i0 + 0x40000000 >> j0; + if (j < i0) + se += 1; + i0 = (j & ~i) | 0x80000000; + i1 = 0; + } + } + } + else if (j0 > 62) + { + if (j0 == 0x4000) + /* Inf or NaN. */ + return x + x; + else + return x; + } + else + { + u_int32_t i = 0xffffffff >> (j0 - 31); + if ((i1 & i) == 0) + /* X is integral. */ + return x; + + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (62 - j0)); + if (j < i1) + { + u_int32_t k = i0 + 1; + if (k < i0) + se += 1; + i0 = k; + } + i1 = j; + } + i1 &= ~i; + } + + SET_LDOUBLE_WORDS (x, se, i0, i1); + return x; +} +weak_alias (__roundl, roundl) diff --git a/sysdeps/libm-ieee754/s_roundtol.c b/sysdeps/libm-ieee754/s_roundtol.c new file mode 100644 index 0000000000..6649369b06 --- /dev/null +++ b/sysdeps/libm-ieee754/s_roundtol.c @@ -0,0 +1,177 @@ +/* Round long double value to long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ + +long int +__roundtol (long double x) +{ + int32_t j0; + u_int32_t i1, i0; + long int result; + + EXTRACT_WORDS (i0, i1, x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) + { + if (j0 < 0) + result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); + else + { + u_int32_t i = 0xfffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long int) ((i0 & 0xfffff) | 0x100000) >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x80000 >> j0); + if (j < i0) + result = (long int) 0x80000 >> (20 - j0); + else + result = (j | 0x100000) >> (20 - j0); + } + } + } + else if (j0 >= 8 * sizeof (long int) || j0 > 51) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long int) x; + } + else + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + if (j < i1) + { + j = i0 + 1; + if ((j & 0xfffff) == 0) + { + if (sizeof (long int) <= 4) + /* Overflow. */ + result = (long int) x; + else + result = 1l << (j0 + 1); + } + else + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + } + else + { + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + + return i0 & 0x80000000 ? -result : result; +} +#else +long int +__roundtol (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long int result; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + result = j0 < -1 ? 0 : 1; + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long int) i0 >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + result = 0x80000000l >> (30 - j0); + else + result = j >> (31 - j0); + } + } + } + else if ((unsigned int) j0 >= 8 * sizeof (long int) || j0 > 62) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + { + j = i0 + 1; + if (j == 0) + { + if (sizeof (long int) <= 4) + /* Overflow. */ + result = (long int) x; + else + result = 1l << (j0 + 1); + } + else + result = (long int) i0 << (j0 - 31); + } + else + { + result = (long int) i0 << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long int) i0 << (j0 - 31); + if (sizeof (long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return se & 0x8000 ? -result : result; +} +#endif +weak_alias (__roundtol, roundtol) diff --git a/sysdeps/libm-ieee754/s_roundtoll.c b/sysdeps/libm-ieee754/s_roundtoll.c new file mode 100644 index 0000000000..8d99130697 --- /dev/null +++ b/sysdeps/libm-ieee754/s_roundtoll.c @@ -0,0 +1,179 @@ +/* Round long double value to long long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +#ifdef NO_LONG_DOUBLE +/* The `long double' is in fact the IEEE `double' type. */ + +long long int +__roundtoll (long double x) +{ + int32_t j0; + u_int32_t i1, i0; + long long int result; + + EXTRACT_WORDS (i0, i1, x); + j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; + if (j0 < 20) + { + if (j0 < 0) + result = j0 < -1 ? 0 : ((i0 & 0x80000000) ? -1 : 1); + else + { + u_int32_t i = 0xfffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long long int) ((i0 & 0xfffff) | 0x100000) >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x80000 >> j0); + if (j < i0) + result = (long long int) 0x80000 >> (20 - j0); + else + result = (j | 0x100000) >> (20 - j0); + } + } + } + else if (j0 >= 8 * sizeof (long long int) || j0 > 51) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long long int) x; + } + else + { + i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + if (j < i1) + { + j = i0 + 1; + if ((j & 0xfffff) == 0) + { + if (sizeof (long long int) <= 4) + /* Overflow. */ + result = (long long int) x; + else + result = 1ll << (j0 + 1); + } + else + result = ((long long int) ((i0 & 0xfffff) | 0x100000) + << (j0 - 31)); + } + else + { + result = ((long long int) ((i0 & 0xfffff) | 0x100000) + << (j0 - 31)); + if (sizeof (long long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long long int) ((i0 & 0xfffff) | 0x100000) << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + + return i0 & 0x80000000 ? -result : result; +} +#else +long long int +__roundtoll (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long long int result; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + result = j0 < -1 ? 0 : 1; + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + result = (long long int) i0 >> j0; + else + { + /* X is not integral. */ + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + result = 0x80000000l >> (30 - j0); + else + result = j >> (31 - j0); + } + } + } + else if ((unsigned int) j0 >= 8 * sizeof (long long int) || j0 > 62) + { + /* The number is too large. It is left implementation defined + what happens. */ + result = (long long int) x; + } + else + { + u_int32_t i = ((u_int32_t) (0xffffffff)) >> (j0 - 31); + if ((i1 & i) != 0) + { + /* x is not integral. */ + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + { + j = i0 + 1; + if (j == 0) + { + if (sizeof (long long int) <= 4) + /* Overflow. */ + result = (long long int) x; + else + result = 1ll << (j0 + 1); + } + else + result = (long long int) i0 << (j0 - 31); + } + else + { + result = (long long int) i0 << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= j >> (63 - j0); + } + } + else + { + result = (long long int) i0 << (j0 - 31); + if (sizeof (long long int) > 4 && j0 > 31) + result |= i1 >> (63 - j0); + } + } + + return se & 0x8000 ? -result : result; +} +#endif +weak_alias (__roundtoll, roundtoll) |