home *** CD-ROM | disk | FTP | other *** search
/ Serving the Web / ServingTheWeb1995.disc1of1.iso / linux / slacksrce / d / libc / libc-4.6 / libc-4 / libc-linux / sysdeps / generic / ldexp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-11-02  |  4.1 KB  |  140 lines
  1. /* Copyright (C) 1992 Free Software Foundation, Inc.
  2. This file is part of the GNU C Library.
  3.  
  4. The GNU C Library is free software; you can redistribute it and/or
  5. modify it under the terms of the GNU Library General Public License as
  6. published by the Free Software Foundation; either version 2 of the
  7. License, or (at your option) any later version.
  8.  
  9. The GNU C Library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12. Library General Public License for more details.
  13.  
  14. You should have received a copy of the GNU Library General Public
  15. License along with the GNU C Library; see the file COPYING.LIB.  If
  16. not, write to the Free Software Foundation, Inc., 675 Mass Ave,
  17. Cambridge, MA 02139, USA.  */
  18.  
  19. /*
  20.  * Copyright (c) 1985 Regents of the University of California.
  21.  * All rights reserved.
  22.  *
  23.  * Redistribution and use in source and binary forms are permitted provided
  24.  * that: (1) source distributions retain this entire copyright notice and
  25.  * comment, and (2) distributions including binaries display the following
  26.  * acknowledgement:  ``This product includes software developed by the
  27.  * University of California, Berkeley and its contributors'' in the
  28.  * documentation or other materials provided with the distribution and in
  29.  * all advertising materials mentioning features or use of this software.
  30.  * Neither the name of the University nor the names of its contributors may
  31.  * be used to endorse or promote products derived from this software without
  32.  * specific prior written permission.
  33.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
  34.  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
  35.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  36.  */
  37.  
  38. #include <ansidecl.h>
  39. #include <math.h>
  40. #include <float.h>
  41. #include <errno.h>
  42. #include "ieee754.h"
  43.  
  44. double
  45. DEFUN(ldexp, (x, exp),
  46.       double x AND int exp)
  47. {
  48.   union ieee754_double *u = (union ieee754_double *) &x;
  49.   unsigned int exponent;
  50.  
  51.   exponent = u->ieee.exponent;
  52.  
  53.   /* The order of the tests is carefully chosen to handle
  54.      the usual case first, with no branches taken.  */
  55.  
  56.   if (exponent != 0)
  57.     {
  58.       /* X is nonzero and not denormalized.  */
  59.  
  60.       if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
  61.       {
  62.       /* X is finite.  When EXP < 0, overflow is actually underflow.  */
  63.  
  64.       exponent += exp;
  65.  
  66.       if (exponent != 0)
  67.         {
  68.           if (exponent <= DBL_MAX_EXP - DBL_MIN_EXP + 1)
  69.         {
  70.           /* In range.  */
  71.           u->ieee.exponent = exponent;
  72.           return x;
  73.         }
  74.  
  75.           if (exp >= 0)
  76.           overflow:
  77.         {
  78.           CONST int negative = u->ieee.negative;
  79.           x = HUGE_VAL;
  80.           u->ieee.negative = negative;
  81.           errno = ERANGE;
  82.           return x;
  83.         }
  84.  
  85.           if (exponent <= - (unsigned int) (DBL_MANT_DIG + 1))
  86.         {
  87.           /* Underflow.  */
  88.           CONST int negative = u->ieee.negative;
  89.           x = 0.0;
  90.           u->ieee.negative = negative;
  91.           errno = ERANGE;
  92.           return x;
  93.         }
  94.         }
  95.  
  96.       /* Gradual underflow.  */
  97.       u->ieee.exponent = 1;
  98.       x *= ldexp (1.0, (int) exponent - 1);
  99.       if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  100.         /* Underflow.  */
  101.         errno = ERANGE;
  102.       return x;
  103.       }
  104.  
  105.       /* X is +-infinity or NaN.  */
  106.       if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  107.       {
  108.       /* X is +-infinity.  */
  109.       if (exp >= 0)
  110.         goto overflow;
  111.       else
  112.         {
  113.           /* (infinity * number < 1).  With infinite precision,
  114.          (infinity / finite) would be infinity, but otherwise it's
  115.          safest to regard (infinity / 2) as indeterminate.  The
  116.          infinity might be (2 * finite).  */
  117.           CONST int negative = u->ieee.negative;
  118.           x = NAN;
  119.           u->ieee.negative = negative;
  120.           errno = EDOM;
  121.           return x;
  122.         }
  123.     }
  124.  
  125.       /* X is NaN.  */
  126.       errno = EDOM;
  127.       return x;
  128.     }
  129.  
  130.   /* X is zero or denormalized.  */
  131.   if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  132.     /* X is +-0.0. */
  133.     return x;
  134.  
  135.   /* X is denormalized.
  136.      Multiplying by 2 ** DBL_MANT_DIG normalizes it;
  137.      we then subtract the DBL_MANT_DIG we added to the exponent.  */
  138.   return ldexp (x * ldexp (1.0, DBL_MANT_DIG), exp - DBL_MANT_DIG);
  139. }
  140.