/* Copyright (C) 1992 Free Software Foundation, Inc.
This file is part of the GNU C Library.

The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with the GNU C Library; see the file COPYING.LIB.  If
not, write to the Free Software Foundation, Inc., 675 Mass Ave,
Cambridge, MA 02139, USA.  */

/*
 * Copyright (c) 1985 Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms are permitted provided
 * that: (1) source distributions retain this entire copyright notice and
 * comment, and (2) distributions including binaries display the following
 * acknowledgement:  ``This product includes software developed by the
 * University of California, Berkeley and its contributors'' in the
 * documentation or other materials provided with the distribution and in
 * all advertising materials mentioning features or use of this software.
 * Neither the name of the University nor the names of its contributors may
 * be used to endorse or promote products derived from this software without
 * specific prior written permission.
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 */

#include <ansidecl.h>
#include <math.h>
#include <float.h>
#include <errno.h>
#include "ieee754.h"

const double ldexp (double x, int exp)
  {
    union ieee754_double *u = (union ieee754_double *) &x;
    unsigned int exponent;

    exponent = u->ieee.exponent;

  /* The order of the tests is carefully chosen to handle
     the usual case first, with no branches taken.  */

    if (exponent != 0)
      {
      /* X is nonzero and not denormalized.  */

	if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
	  {
	  /* X is finite.  When EXP < 0, overflow is actually underflow.  */

	    exponent += exp;

	    if (exponent != 0)
	      {
		if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
		  {
		    /* In range.  */
		    u->ieee.exponent = exponent;
		    return x;
		  }

		if (exp >= 0)
		overflow:
		  {
		    CONST int negative = u->ieee.negative;
		    x = HUGE_VAL;
		    u->ieee.negative = negative;
		    errno = ERANGE;
		    return x;
		  }

		if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
		  {
		    /* Underflow.  */
		    CONST int negative = u->ieee.negative;
		    x = 0.0;
		    u->ieee.negative = negative;
		    errno = ERANGE;
		    return x;
		  }
	      }

	    /* Gradual underflow.  */
	    u->ieee.exponent = 1;
	    x *= ldexp (1.0, (int) exponent - 1);
	    if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
	      /* Underflow.  */
	      errno = ERANGE;
	    return x;
	  }

	/* X is +-infinity or NaN.  */
	if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
	  {
	    /* X is +-infinity.  */
	    if (exp >= 0)
	      goto overflow;
	    else
	      {
		/* (infinity * number < 1).  With infinite precision,
		   (infinity / finite) would be infinity, but otherwise it's
		   safest to regard (infinity / 2) as indeterminate.  The
		   infinity might be (2 * finite).  */
		CONST int negative = u->ieee.negative;
		x = NAN;
		u->ieee.negative = negative;
		errno = EDOM;
		return x;
	      }
	  }
	
	errno = EDOM;
	return x;        /* X is NaN.  */
      }
    
    /* X is zero or denormalized.  */
    if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
      {
	return x;       /* X is +-0.0. */
      }

    /* X is denormalized.
       Multiplying by 2 ** DBL_MANT_DIG normalizes it;
       we then subtract the DBL_MANT_DIG we added to the exponent.  */
    return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
  }
