multc3.c 2.3 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980
  1. /* Copyright (C) 2005-2019 Free Software Foundation, Inc.
  2. This file is part of the GNU C Library.
  3. Contributed by Richard Henderson <rth@redhat.com>, 2005.
  4. The GNU C Library is free software; you can redistribute it and/or
  5. modify it under the terms of the GNU Lesser General Public
  6. License as published by the Free Software Foundation; either
  7. version 2.1 of the License, or (at your option) any later version.
  8. The GNU C Library is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  11. Lesser General Public License for more details.
  12. You should have received a copy of the GNU Lesser General Public
  13. License along with the GNU C Library; if not, see
  14. <http://www.gnu.org/licenses/>. */
  15. #include <stdbool.h>
  16. #include <math.h>
  17. #include <complex.h>
  18. attribute_hidden
  19. long double _Complex
  20. __multc3 (long double a, long double b, long double c, long double d)
  21. {
  22. long double ac, bd, ad, bc, x, y;
  23. ac = a * c;
  24. bd = b * d;
  25. ad = a * d;
  26. bc = b * c;
  27. x = ac - bd;
  28. y = ad + bc;
  29. if (isnan (x) && isnan (y))
  30. {
  31. /* Recover infinities that computed as NaN + iNaN. */
  32. bool recalc = 0;
  33. if (isinf (a) || isinf (b))
  34. {
  35. /* z is infinite. "Box" the infinity and change NaNs in
  36. the other factor to 0. */
  37. a = copysignl (isinf (a) ? 1 : 0, a);
  38. b = copysignl (isinf (b) ? 1 : 0, b);
  39. if (isnan (c)) c = copysignl (0, c);
  40. if (isnan (d)) d = copysignl (0, d);
  41. recalc = 1;
  42. }
  43. if (isinf (c) || isinf (d))
  44. {
  45. /* w is infinite. "Box" the infinity and change NaNs in
  46. the other factor to 0. */
  47. c = copysignl (isinf (c) ? 1 : 0, c);
  48. d = copysignl (isinf (d) ? 1 : 0, d);
  49. if (isnan (a)) a = copysignl (0, a);
  50. if (isnan (b)) b = copysignl (0, b);
  51. recalc = 1;
  52. }
  53. if (!recalc
  54. && (isinf (ac) || isinf (bd)
  55. || isinf (ad) || isinf (bc)))
  56. {
  57. /* Recover infinities from overflow by changing NaNs to 0. */
  58. if (isnan (a)) a = copysignl (0, a);
  59. if (isnan (b)) b = copysignl (0, b);
  60. if (isnan (c)) c = copysignl (0, c);
  61. if (isnan (d)) d = copysignl (0, d);
  62. recalc = 1;
  63. }
  64. if (recalc)
  65. {
  66. x = INFINITY * (a * c - b * d);
  67. y = INFINITY * (a * d + b * c);
  68. }
  69. }
  70. return x + I * y;
  71. }