lib/ldexp.c

   1 /* Multiply a 'float' by a power of 2.
   2    Copyright 2002-2003, 2007-2024 Free Software Foundation, Inc.
   3
   4    This file is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU Lesser General Public License as
   6    published by the Free Software Foundation; either version 2.1 of the
   7    License, or (at your option) any later version.
   8
   9    This file 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
  12    GNU Lesser General Public License for more details.
  13
  14    You should have received a copy of the GNU Lesser General Public License
  15    along with this program.  If not, see <https://www.gnu.org/licenses/>.  */
  16
  17 /* Written by Paolo Bonzini and Bruno Haible.  */
  18
  19 #include <config.h>
  20
  21 /* Specification.  */
  22 #include <math.h>
  23
  24 # include <float.h>
  25 #ifdef USE_LONG_DOUBLE
  26 # include "fpucw.h"
  27 #endif
  28
  29 /* Avoid some warnings from "gcc -Wshadow".
  30    This file doesn't use the exp() function.  */
  31 #undef exp
  32 #define exp exponent
  33
  34 #ifdef USE_LONG_DOUBLE
  35 # define FUNC ldexpl
  36 # define DOUBLE long double
  37 # define ISNAN isnanl
  38 # define DECL_ROUNDING DECL_LONG_DOUBLE_ROUNDING
  39 # define BEGIN_ROUNDING() BEGIN_LONG_DOUBLE_ROUNDING ()
  40 # define END_ROUNDING() END_LONG_DOUBLE_ROUNDING ()
  41 # define L_(literal) literal##L
  42 #else
  43 # define FUNC ldexp
  44 # define DOUBLE double
  45 # define ISNAN isnand
  46 # define DECL_ROUNDING
  47 # define BEGIN_ROUNDING()
  48 # define END_ROUNDING()
  49 # define L_(literal) literal
  50 #endif
  51
  52 DOUBLE
  53 FUNC (DOUBLE x, int exp)
  54 {
  55   DECL_ROUNDING
  56
  57   BEGIN_ROUNDING ();
  58
  59   /* Check for zero, nan and infinity. */
  60   if (!(ISNAN (x) || x + x == x))
  61     {
  62       DOUBLE factor = exp < 0 ? L_(0.5) : L_(2.0);
  63
  64       while (true)
  65         {
  66           if (exp % 2 != 0)
  67             x *= factor;
  68           exp /= 2;
  69           if (exp == 0)
  70             break;
  71           factor = factor * factor;
  72         }
  73     }
  74
  75   END_ROUNDING ();
  76
  77   return x;
  78 }