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  __divtc3 (long double a, long double b, long double c, long double d)
      25  {
      26    long double denom, ratio, x, y;
      27  
      28    /* ??? We can get better behavior from logarithmic scaling instead of
      29       the division.  But that would mean starting to link libgcc against
      30       libm.  We could implement something akin to ldexp/frexp as gcc builtins
      31       fairly easily...  */
      32    if (fabsl (c) < fabsl (d))
      33      {
      34        ratio = c / d;
      35        denom = (c * ratio) + d;
      36        x = ((a * ratio) + b) / denom;
      37        y = ((b * ratio) - a) / denom;
      38      }
      39    else
      40      {
      41        ratio = d / c;
      42        denom = (d * ratio) + c;
      43        x = ((b * ratio) + a) / denom;
      44        y = (b - (a * ratio)) / denom;
      45      }
      46  
      47    /* Recover infinities and zeros that computed as NaN+iNaN; the only cases
      48       are nonzero/zero, infinite/finite, and finite/infinite.  */
      49    if (isnan (x) && isnan (y))
      50      {
      51        if (denom == 0.0 && (!isnan (a) || !isnan (b)))
      52  	{
      53  	  x = copysignl (INFINITY, c) * a;
      54  	  y = copysignl (INFINITY, c) * b;
      55  	}
      56        else if ((isinf (a) || isinf (b))
      57  	       && isfinite (c) && isfinite (d))
      58  	{
      59  	  a = copysignl (isinf (a) ? 1 : 0, a);
      60  	  b = copysignl (isinf (b) ? 1 : 0, b);
      61  	  x = INFINITY * (a * c + b * d);
      62  	  y = INFINITY * (b * c - a * d);
      63  	}
      64        else if ((isinf (c) || isinf (d))
      65  	       && isfinite (a) && isfinite (b))
      66  	{
      67  	  c = copysignl (isinf (c) ? 1 : 0, c);
      68  	  d = copysignl (isinf (d) ? 1 : 0, d);
      69  	  x = 0.0 * (a * c + b * d);
      70  	  y = 0.0 * (b * c - a * d);
      71  	}
      72      }
      73  
      74    return x + I * y;
      75  }