1  /* Emulation for ldexpl.
       2     Contributed by Paolo Bonzini
       3  
       4     Copyright 2002-2003, 2007-2021 Free Software Foundation, Inc.
       5  
       6     This file is part of gnulib.
       7  
       8     This file is free software: you can redistribute it and/or modify
       9     it under the terms of the GNU Lesser General Public License as
      10     published by the Free Software Foundation; either version 3 of the
      11     License, or (at your option) any later version.
      12  
      13     This file is distributed in the hope that it will be useful,
      14     but WITHOUT ANY WARRANTY; without even the implied warranty of
      15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      16     GNU Lesser General Public License for more details.
      17  
      18     You should have received a copy of the GNU Lesser General Public License
      19     along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
      20  
      21  #include <config.h>
      22  
      23  /* Specification.  */
      24  #include <math.h>
      25  
      26  #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
      27  
      28  long double
      29  ldexpl (long double x, int exp)
      30  {
      31    return ldexp (x, exp);
      32  }
      33  
      34  #else
      35  
      36  # include <float.h>
      37  # include "fpucw.h"
      38  
      39  long double
      40  ldexpl (long double x, int exp)
      41  {
      42    long double factor;
      43    int bit;
      44    DECL_LONG_DOUBLE_ROUNDING
      45  
      46    BEGIN_LONG_DOUBLE_ROUNDING ();
      47  
      48    /* Check for zero, nan and infinity. */
      49    if (!(isnanl (x) || x + x == x))
      50      {
      51        if (exp < 0)
      52          {
      53            exp = -exp;
      54            factor = 0.5L;
      55          }
      56        else
      57          factor = 2.0L;
      58  
      59        if (exp > 0)
      60          for (bit = 1;;)
      61            {
      62              /* Invariant: Here bit = 2^i, factor = 2^-2^i or = 2^2^i,
      63                 and bit <= exp.  */
      64              if (exp & bit)
      65                x *= factor;
      66              bit <<= 1;
      67              if (bit > exp)
      68                break;
      69              factor = factor * factor;
      70            }
      71      }
      72  
      73    END_LONG_DOUBLE_ROUNDING ();
      74  
      75    return x;
      76  }
      77  
      78  #endif
      79  
      80  #if 0
      81  int
      82  main (void)
      83  {
      84    long double x;
      85    int y;
      86    for (y = 0; y < 29; y++)
      87      printf ("%5d %.16Lg %.16Lg\n", y, ldexpl (0.8L, y), ldexpl (0.8L, -y) * ldexpl (0.8L, y));
      88  }
      89  #endif