/*libc/ieee_float/ldexp.cCreated:	Oct 14, 1993 by Philip Homburg <philip@cs.vu.nl>Implementation of ldexp that directly manipulates the exponent bits in anieee 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 $ */