home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / ldexpl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-26  |  4.2 KB  |  143 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 "mathl.h"
  40. #include <float.h>
  41. #include <errno.h>
  42. #include "ieee854.h"
  43. long double __infnanl (int);
  44.  
  45. long double
  46. DEFUN(ldexpl, (x, exp),
  47.       long double x AND int exp)
  48. {
  49.   union ieee854_double *u = (union ieee854_double *) &x;
  50.   unsigned int exponent;
  51.  
  52.   exponent = u->ieee.exponent;
  53.  
  54.   /* The order of the tests is carefully chosen to handle
  55.      the usual case first, with no branches taken.  */
  56.  
  57.   if (exponent != 0)
  58.     {
  59.       /* X is nonzero and not denormalized.  */
  60.  
  61.       if (exponent <= LDBL_MAX_EXP - LDBL_MIN_EXP + 1)
  62.       {
  63.       /* X is finite.  When EXP < 0, overflow is actually underflow.  */
  64.  
  65.       exponent += exp;
  66.  
  67.       if (exponent != 0)
  68.         {
  69.           if (exponent <= LDBL_MAX_EXP - LDBL_MIN_EXP + 1)
  70.         {
  71.           /* In range.  */
  72.           u->ieee.exponent = exponent;
  73.           return x;
  74.         }
  75.  
  76.           if (exp >= 0)
  77.           overflow:
  78.         {
  79.           CONST int negative = u->ieee.negative;
  80. /*          x = HUGE_VAL; */
  81.           x = 1.189731495357231765021263853E4932L;
  82.           u->ieee.negative = negative;
  83.           errno = ERANGE;
  84.           return x;
  85.         }
  86.  
  87.           if (exponent <= - (unsigned int) (LDBL_MANT_DIG + 1))
  88.         {
  89.           /* Underflow.  */
  90.           CONST int negative = u->ieee.negative;
  91.           x = 0.0L;
  92.           u->ieee.negative = negative;
  93.           errno = ERANGE;
  94.           return x;
  95.         }
  96.         }
  97.  
  98.       /* Gradual underflow.  */
  99.       u->ieee.exponent = 1;
  100.       x *= ldexpl (1.0L, (int) exponent - 1);
  101.       if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  102.         /* Underflow.  */
  103.         errno = ERANGE;
  104.       return x;
  105.       }
  106.  
  107.       /* X is +-infinity or NaN.  */
  108.       if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  109.       {
  110.       /* X is +-infinity.  */
  111.       if (exp >= 0)
  112.         goto overflow;
  113.       else
  114.         {
  115.           /* (infinity * number < 1).  With infinite precision,
  116.          (infinity / finite) would be infinity, but otherwise it's
  117.          safest to regard (infinity / 2) as indeterminate.  The
  118.          infinity might be (2 * finite).  */
  119.           CONST int negative = u->ieee.negative;
  120.           /* Create a NaN.  */
  121.           x = __infnanl(EDOM);
  122.           u->ieee.negative = negative;
  123.           errno = EDOM;
  124.           return x;
  125.         }
  126.     }
  127.  
  128.       /* X is NaN.  */
  129.       errno = EDOM;
  130.       return x;
  131.     }
  132.  
  133.   /* X is zero or denormalized.  */
  134.   if (u->ieee.mantissa0 == 0 && u->ieee.mantissa1 == 0)
  135.     /* X is +-0.0. */
  136.     return x;
  137.  
  138.   /* X is denormalized.
  139.      Multiplying by 2 ** LDBL_MANT_DIG normalizes it;
  140.      we then subtract the LDBL_MANT_DIG we added to the exponent.  */
  141.   return ldexpl (x * ldexpl (1.0L, LDBL_MANT_DIG), exp - LDBL_MANT_DIG);
  142. }
  143.