htswanalysis/MACS/lib/gsl/gsl-1.11/sys/ldfrexp.c

   1 /* sys/ldfrexp.c
   2  *
   3  * Copyright (C) 2002, Gert Van den Eynde
   4  * Copyright (C) 2007, Brian Gough
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <math.h>
  23 #include <gsl/gsl_math.h>
  24
  25 double
  26 gsl_ldexp (const double x, const int e)
  27 {
  28   int ex;
  29
  30   if (x == 0.0)
  31     {
  32       return x;
  33     }
  34
  35   {
  36     double y = gsl_frexp (x, &ex);
  37     double e2 = e + ex, p2;
  38
  39     if (e2 >= DBL_MAX_EXP)
  40       {
  41         y *= pow (2.0, e2 - DBL_MAX_EXP + 1);
  42         e2 = DBL_MAX_EXP - 1;
  43       }
  44     else if (e2 <= DBL_MIN_EXP)
  45       {
  46         y *= pow (2.0, e2 - DBL_MIN_EXP - 1);
  47         e2 = DBL_MIN_EXP + 1;
  48       }
  49
  50     p2 = pow (2.0, e2);
  51     return y * p2;
  52   }
  53 }
  54
  55 double
  56 gsl_frexp (const double x, int *e)
  57 {
  58   if (x == 0.0)
  59     {
  60       *e = 0;
  61       return 0.0;
  62     }
  63   else if (!finite (x))
  64     {
  65       *e = 0;
  66       return x;
  67     }
  68   else if (fabs (x) >= 0.5 && fabs (x) < 1)     /* Handle the common case */
  69     {
  70       *e = 0;
  71       return x;
  72     }
  73   else
  74     {
  75       double ex = ceil (log (fabs (x)) / M_LN2);
  76       int ei = (int) ex;
  77       double f;
  78
  79       /* Prevent underflow and overflow of 2**(-ei) */
  80       if (ei < DBL_MIN_EXP)
  81         ei = DBL_MIN_EXP;
  82
  83       if (ei > -DBL_MIN_EXP)
  84         ei = -DBL_MIN_EXP;
  85
  86       f = x * pow (2.0, -ei);
  87
  88       if (!finite (f))
  89         {
  90           /* This should not happen */
  91           *e = 0;
  92           return f;
  93         }
  94
  95       while (fabs (f) >= 1.0)
  96         {
  97           ei++;
  98           f /= 2.0;
  99         }
 100
 101       while (fabs (f) > 0 && fabs (f) < 0.5)
 102         {
 103           ei--;
 104           f *= 2.0;
 105         }
 106
 107       *e = ei;
 108       return f;
 109     }
 110 }