1  /* Copyright (C) 2005-2023 Free Software Foundation, Inc.
       2     This file is part of the GNU C Library.
       3  
       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  
       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  
      14     You should have received a copy of the GNU Lesser General Public
      15     License along with the GNU C Library; if not, see
      16     <https://www.gnu.org/licenses/>.  */
      17  
      18  #include <stdbool.h>
      19  #include <math.h>
      20  #include <complex.h>
      21  
      22  attribute_hidden
      23  long double _Complex
      24  __multc3 (long double a, long double b, long double c, long double d)
      25  {
      26    long double ac, bd, ad, bc, x, y;
      27  
      28    ac = a * c;
      29    bd = b * d;
      30    ad = a * d;
      31    bc = b * c;
      32  
      33    x = ac - bd;
      34    y = ad + bc;
      35  
      36    if (isnan (x) && isnan (y))
      37      {
      38        /* Recover infinities that computed as NaN + iNaN.  */
      39        bool recalc = 0;
      40        if (isinf (a) || isinf (b))
      41  	{
      42  	  /* z is infinite.  "Box" the infinity and change NaNs in
      43  	     the other factor to 0.  */
      44  	  a = copysignl (isinf (a) ? 1 : 0, a);
      45  	  b = copysignl (isinf (b) ? 1 : 0, b);
      46  	  if (isnan (c)) c = copysignl (0, c);
      47  	  if (isnan (d)) d = copysignl (0, d);
      48  	  recalc = 1;
      49  	}
      50       if (isinf (c) || isinf (d))
      51  	{
      52  	  /* w is infinite.  "Box" the infinity and change NaNs in
      53  	     the other factor to 0.  */
      54  	  c = copysignl (isinf (c) ? 1 : 0, c);
      55  	  d = copysignl (isinf (d) ? 1 : 0, d);
      56  	  if (isnan (a)) a = copysignl (0, a);
      57  	  if (isnan (b)) b = copysignl (0, b);
      58  	  recalc = 1;
      59  	}
      60       if (!recalc
      61  	  && (isinf (ac) || isinf (bd)
      62  	      || isinf (ad) || isinf (bc)))
      63  	{
      64  	  /* Recover infinities from overflow by changing NaNs to 0.  */
      65  	  if (isnan (a)) a = copysignl (0, a);
      66  	  if (isnan (b)) b = copysignl (0, b);
      67  	  if (isnan (c)) c = copysignl (0, c);
      68  	  if (isnan (d)) d = copysignl (0, d);
      69  	  recalc = 1;
      70  	}
      71        if (recalc)
      72  	{
      73  	  x = INFINITY * (a * c - b * d);
      74  	  y = INFINITY * (a * d + b * c);
      75  	}
      76      }
      77  
      78    return x + I * y;
      79  }