1/* mpfr_scale2 -- multiply a double float by 2^exp
2
3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include <float.h> /* for DBL_EPSILON */
24#include "mpfr-impl.h"
25
26/* Note: we could use the ldexp function, but since we want not to depend on
27   math.h, we write our own implementation. */
28
29/* multiplies 1/2 <= d <= 1 by 2^exp */
30double
31mpfr_scale2 (double d, int exp)
32{
33#if _GMP_IEEE_FLOATS
34  {
35    union ieee_double_extract x;
36
37    if (MPFR_UNLIKELY (d == 1.0))
38      {
39        d = 0.5;
40        exp ++;
41      }
42
43    /* now 1/2 <= d < 1 */
44
45    /* infinities and zeroes have already been checked */
46    MPFR_ASSERTD (-1073 <= exp && exp <= 1025);
47
48    x.d = d;
49    if (MPFR_UNLIKELY (exp < -1021)) /* subnormal case */
50      {
51        x.s.exp += exp + 52;
52        x.d *= DBL_EPSILON;
53      }
54    else /* normalized case */
55      {
56        x.s.exp += exp;
57      }
58    return x.d;
59  }
60#else /* _GMP_IEEE_FLOATS */
61  {
62    double factor;
63
64    /* An overflow may occurs (example: 0.5*2^1024) */
65    if (d < 1.0)
66      {
67        d += d;
68        exp--;
69      }
70    /* Now 1.0 <= d < 2.0 */
71
72    if (exp < 0)
73      {
74        factor = 0.5;
75        exp = -exp;
76      }
77    else
78      {
79        factor = 2.0;
80      }
81    while (exp != 0)
82      {
83        if ((exp & 1) != 0)
84          d *= factor;
85        exp >>= 1;
86        factor *= factor;
87      }
88    return d;
89  }
90#endif
91}
92