|
Packit |
67cb25 |
/* sys/ldfrexp.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2002, Gert Van den Eynde
|
|
Packit |
67cb25 |
* Copyright (C) 2007, Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ldexp (const double x, const int e)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int ex;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double y = gsl_frexp (x, &ex);
|
|
Packit |
67cb25 |
double e2 = e + ex, p2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (e2 >= DBL_MAX_EXP)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
y *= pow (2.0, e2 - DBL_MAX_EXP + 1);
|
|
Packit |
67cb25 |
e2 = DBL_MAX_EXP - 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (e2 <= DBL_MIN_EXP)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
y *= pow (2.0, e2 - DBL_MIN_EXP - 1);
|
|
Packit |
67cb25 |
e2 = DBL_MIN_EXP + 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
p2 = pow (2.0, e2);
|
|
Packit |
67cb25 |
return y * p2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_frexp (const double x, int *e)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (x == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*e = 0;
|
|
Packit |
67cb25 |
return 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (!gsl_finite (x))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*e = 0;
|
|
Packit |
67cb25 |
return x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (fabs (x) >= 0.5 && fabs (x) < 1) /* Handle the common case */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*e = 0;
|
|
Packit |
67cb25 |
return x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ex = ceil (log (fabs (x)) / M_LN2);
|
|
Packit |
67cb25 |
int ei = (int) ex;
|
|
Packit |
67cb25 |
double f;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Prevent underflow and overflow of 2**(-ei) */
|
|
Packit |
67cb25 |
if (ei < DBL_MIN_EXP)
|
|
Packit |
67cb25 |
ei = DBL_MIN_EXP;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ei > -DBL_MIN_EXP)
|
|
Packit |
67cb25 |
ei = -DBL_MIN_EXP;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
f = x * pow (2.0, -ei);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!gsl_finite (f))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* This should not happen */
|
|
Packit |
67cb25 |
*e = 0;
|
|
Packit |
67cb25 |
return f;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (fabs (f) >= 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ei++;
|
|
Packit |
67cb25 |
f /= 2.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while (fabs (f) > 0 && fabs (f) < 0.5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ei--;
|
|
Packit |
67cb25 |
f *= 2.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*e = ei;
|
|
Packit |
67cb25 |
return f;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|