glibc/sysdeps/ieee754/ldexp.c
1995-02-18 01:27:10 +00:00

146 lines
4.2 KiB
C

/* Copyright (C) 1992, 1995 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 "ieee754.h"
double
DEFUN(ldexp, (x, exp),
double x AND int exp)
{
union ieee754_double u;
unsigned int exponent;
u.d = x;
#define x u.d
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;
u.d = HUGE_VAL;
u.ieee.negative = negative;
errno = ERANGE;
return u.d;
}
if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
{
/* Underflow. */
CONST int negative = u.ieee.negative;
u.d = 0.0;
u.ieee.negative = negative;
errno = ERANGE;
return u.d;
}
}
/* Gradual underflow. */
u.ieee.exponent = 1;
u.d *= ldexp (1.0, (int) exponent - 1);
if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
/* Underflow. */
errno = ERANGE;
return u.d;
}
/* 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;
u.d = NAN;
u.ieee.negative = negative;
errno = EDOM;
return u.d;
}
}
/* X is NaN. */
errno = EDOM;
return u.d;
}
/* X is zero or denormalized. */
if (u.ieee.mantissa0 == 0 && u.ieee.mantissa1 == 0)
/* X is +-0.0. */
return x;
/* 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);
}
/* Compatibility names for the same function. */
weak_alias (ldexp, __scalb)
weak_alias (ldexp, scalb)