minix/lib/gnu/ieee_float/ldexp.c

100 lines
1.7 KiB
C
Raw Normal View History

/*
libc/ieee_float/ldexp.c
Created: Oct 14, 1993 by Philip Homburg <philip@cs.vu.nl>
Implementation of ldexp that directly manipulates the exponent bits in an
ieee float
*/
#include <sys/types.h>
#include <errno.h>
#include <math.h>
#include "ieee_float.h"
double ldexp(value, exp)
double value;
int exp;
{
struct f64 *f64p;
int oldexp, exp_bias;
double factor;
f64p= (struct f64 *)&value;
exp_bias= 0;
oldexp= F64_GET_EXP(f64p);
if (oldexp == F64_EXP_MAX)
{ /* Either infinity or Nan */
return value;
}
if (oldexp == 0)
{
/* Either 0 or denormal */
if (F64_GET_MANT_LOW(f64p) == 0 &&
F64_GET_MANT_HIGH(f64p) == 0)
{
return value;
}
}
/* If exp is too large (> 2*F64_EXP_MAX) or too small
* (< -2*F64_EXP_MAX) return HUGE_VAL or 0. This prevents overflows
* in exp if exp is really weird
*/
if (exp >= 2*F64_EXP_MAX)
{
errno= ERANGE;
return HUGE_VAL;
}
if (exp <= -2*F64_EXP_MAX)
{
errno= ERANGE;
return 0;
}
/* Normalize a denormal */
if (oldexp == 0)
{
/* Multiply by 2^64 */
factor= 65536.0; /* 2^16 */
factor *= factor; /* 2^32 */
factor *= factor; /* 2^64 */
value *= factor;
exp= -64;
oldexp= F64_GET_EXP(f64p);
}
exp= oldexp + exp;
if (exp >= F64_EXP_MAX)
{ /* Overflow */
errno= ERANGE;
return HUGE_VAL;
}
if (exp > 0)
{
/* Normal */
F64_SET_EXP(f64p, exp);
return value;
}
/* Denormal, or underflow. */
exp += 64;
F64_SET_EXP(f64p, exp);
/* Divide by 2^64 */
factor= 65536.0; /* 2^16 */
factor *= factor; /* 2^32 */
factor *= factor; /* 2^64 */
value /= factor;
if (value == 0.0)
{
/* Underflow */
errno= ERANGE;
}
return value;
}
/*
* $PchId: ldexp.c,v 1.3 1996/02/22 21:01:39 philip Exp $
*/