#include "ieeeftn.h"

#if STDC
DOUBLEPRECISION
dsetxp(DOUBLEPRECISION *x, INTEGER *n)
#else /* NOT STDC */
DOUBLEPRECISION
dsetxp(x,n)			/* return (fraction of x) * base**n */
DOUBLEPRECISION *x;		/* i.e. replace the exponent field of x */
INTEGER *n;
#endif /* STDC */
{
    DOUBLEPRECISION_PARTS w;

    if (disden(x))			/* then x is denormalized */
        w.r = *x * BASE_TO_THE_T_DP;	/* so must first normalize fraction */
    else
	w.r = *x;			/* x is normal, NaN, or Inf */
    if (w.r == 0.0)			/* zeroes must be left intact */
        /* NO-OP */;
    else if (disnan(x))		/* NaNs must be left intact */
        /* NO-OP */;
    else if (disinf(x))		/* Infs must be left intact */
        /* NO-OP */;
    else if (*n <= (EXPONENT_DENORM_DP - T_DP)) /* underflow to zero */
	w.r = (DOUBLEPRECISION)0.0;
    else if (*n <= EXPONENT_DENORM_DP)	/* generate denormalized value */
    {
	w.i[DP_HIGH] &= ~EXPONENT_MASK_DP; /* clear exponent field */
	w.i[DP_HIGH] |= SET_EXPONENT_DP(*n + T_DP); /* set new */
					/* exponent of scaled normal value */
	w.r /= BASE_TO_THE_T_DP;	/* and denormalize by division */
    }
    else if (*n >= EXPONENT_INFNAN_DP)	/* generate infinity */
    {
	w.i[DP_HIGH] = ISNEG_DP(w.i[DP_HIGH]) ? NegInf_DP : Inf_DP;
	w.i[DP_LOW] = ISNEG_DP(w.i[DP_HIGH]) ? NegInf_Low_DP : Inf_Low_DP;
    }
    else				/* result is normal number */
    {					/* and exponent in range */
	w.i[DP_HIGH] &= ~EXPONENT_MASK_DP;	/* clear exponent field */
	w.i[DP_HIGH] |= SET_EXPONENT_DP(*n);	/* set new exponent */
    }

    return (w.r);
}