/math/s_catanh_template.c

https://github.com/hjl-tools/glibc · C · 136 lines · 105 code · 13 blank · 18 comment · 34 complexity · ee51a28c1b8e23a51599f186b57d5d18 MD5 · raw file

  1. /* Return arc hyperbolic tangent for a complex float type.
  2. Copyright (C) 1997-2020 Free Software Foundation, Inc.
  3. This file is part of the GNU C Library.
  4. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
  5. The GNU C Library is free software; you can redistribute it and/or
  6. modify it under the terms of the GNU Lesser General Public
  7. License as published by the Free Software Foundation; either
  8. version 2.1 of the License, or (at your option) any later version.
  9. The GNU C Library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  12. Lesser General Public License for more details.
  13. You should have received a copy of the GNU Lesser General Public
  14. License along with the GNU C Library; if not, see
  15. <https://www.gnu.org/licenses/>. */
  16. #include <complex.h>
  17. #include <math.h>
  18. #include <math_private.h>
  19. #include <math-underflow.h>
  20. #include <float.h>
  21. CFLOAT
  22. M_DECL_FUNC (__catanh) (CFLOAT x)
  23. {
  24. CFLOAT res;
  25. int rcls = fpclassify (__real__ x);
  26. int icls = fpclassify (__imag__ x);
  27. if (__glibc_unlikely (rcls <= FP_INFINITE || icls <= FP_INFINITE))
  28. {
  29. if (icls == FP_INFINITE)
  30. {
  31. __real__ res = M_COPYSIGN (0, __real__ x);
  32. __imag__ res = M_COPYSIGN (M_MLIT (M_PI_2), __imag__ x);
  33. }
  34. else if (rcls == FP_INFINITE || rcls == FP_ZERO)
  35. {
  36. __real__ res = M_COPYSIGN (0, __real__ x);
  37. if (icls >= FP_ZERO)
  38. __imag__ res = M_COPYSIGN (M_MLIT (M_PI_2), __imag__ x);
  39. else
  40. __imag__ res = M_NAN;
  41. }
  42. else
  43. {
  44. __real__ res = M_NAN;
  45. __imag__ res = M_NAN;
  46. }
  47. }
  48. else if (__glibc_unlikely (rcls == FP_ZERO && icls == FP_ZERO))
  49. {
  50. res = x;
  51. }
  52. else
  53. {
  54. if (M_FABS (__real__ x) >= 16 / M_EPSILON
  55. || M_FABS (__imag__ x) >= 16 / M_EPSILON)
  56. {
  57. __imag__ res = M_COPYSIGN (M_MLIT (M_PI_2), __imag__ x);
  58. if (M_FABS (__imag__ x) <= 1)
  59. __real__ res = 1 / __real__ x;
  60. else if (M_FABS (__real__ x) <= 1)
  61. __real__ res = __real__ x / __imag__ x / __imag__ x;
  62. else
  63. {
  64. FLOAT h = M_HYPOT (__real__ x / 2, __imag__ x / 2);
  65. __real__ res = __real__ x / h / h / 4;
  66. }
  67. }
  68. else
  69. {
  70. if (M_FABS (__real__ x) == 1
  71. && M_FABS (__imag__ x) < M_EPSILON * M_EPSILON)
  72. __real__ res = (M_COPYSIGN (M_LIT (0.5), __real__ x)
  73. * ((FLOAT) M_MLIT (M_LN2)
  74. - M_LOG (M_FABS (__imag__ x))));
  75. else
  76. {
  77. FLOAT i2 = 0;
  78. if (M_FABS (__imag__ x) >= M_EPSILON * M_EPSILON)
  79. i2 = __imag__ x * __imag__ x;
  80. FLOAT num = 1 + __real__ x;
  81. num = i2 + num * num;
  82. FLOAT den = 1 - __real__ x;
  83. den = i2 + den * den;
  84. FLOAT f = num / den;
  85. if (f < M_LIT (0.5))
  86. __real__ res = M_LIT (0.25) * M_LOG (f);
  87. else
  88. {
  89. num = 4 * __real__ x;
  90. __real__ res = M_LIT (0.25) * M_LOG1P (num / den);
  91. }
  92. }
  93. FLOAT absx, absy, den;
  94. absx = M_FABS (__real__ x);
  95. absy = M_FABS (__imag__ x);
  96. if (absx < absy)
  97. {
  98. FLOAT t = absx;
  99. absx = absy;
  100. absy = t;
  101. }
  102. if (absy < M_EPSILON / 2)
  103. {
  104. den = (1 - absx) * (1 + absx);
  105. if (den == 0)
  106. den = 0;
  107. }
  108. else if (absx >= 1)
  109. den = (1 - absx) * (1 + absx) - absy * absy;
  110. else if (absx >= M_LIT (0.75) || absy >= M_LIT (0.5))
  111. den = -M_SUF (__x2y2m1) (absx, absy);
  112. else
  113. den = (1 - absx) * (1 + absx) - absy * absy;
  114. __imag__ res = M_LIT (0.5) * M_ATAN2 (2 * __imag__ x, den);
  115. }
  116. math_check_force_underflow_complex (res);
  117. }
  118. return res;
  119. }
  120. declare_mgen_alias (__catanh, catanh)