newlib/libm/common/frexpl.c

   1 /*
   2 (C) Copyright IBM Corp. 2009
   3
   4 All rights reserved.
   5
   6 Redistribution and use in source and binary forms, with or without
   7 modification, are permitted provided that the following conditions are met:
   8
   9 * Redistributions of source code must retain the above copyright notice,
  10 this list of conditions and the following disclaimer.
  11 * Redistributions in binary form must reproduce the above copyright
  12 notice, this list of conditions and the following disclaimer in the
  13 documentation and/or other materials provided with the distribution.
  14 * Neither the name of IBM nor the names of its contributors may be
  15 used to endorse or promote products derived from this software without
  16 specific prior written permission.
  17
  18 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  22 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 POSSIBILITY OF SUCH DAMAGE.
  29 */
  30
  31 #include <math.h>
  32 #include <float.h>
  33 #include "local.h"
  34
  35 /* On platforms where long double is as wide as double.  */
  36 #if defined(_LDBL_EQ_DBL)
  37 long double
  38 frexpl (long double x, int *eptr)
  39 {
  40   return frexp(x, eptr);
  41 }
  42 #else  /* !_DBL_EQ_DBL */
  43 # if (LDBL_MANT_DIG == 53) /* 64-bit long double */
  44 static const double scale = 0x1p54;
  45
  46 union ldbl {
  47   long double x;
  48   struct {
  49 #  ifdef __IEEE_LITTLE_ENDIAN /* for Intel CPU */
  50     __uint32_t fracl;
  51     __uint32_t frach:20;
  52     __uint32_t exp:11;
  53     __uint32_t sign:1;
  54 #  endif
  55 #  ifdef __IEEE_BIG_ENDIAN
  56     __uint32_t sign:1;
  57     __uint32_t exp:11;
  58     __uint32_t frach:20;
  59 #   ifndef ___IEEE_BYTES_LITTLE_ENDIAN
  60 #   else /* ARMEL without __VFP_FP__ */
  61     __uint32_t frach:20;
  62     __uint32_t exp:11;
  63     __uint32_t sign:1;
  64 #   endif
  65     __uint32_t fracl;
  66 #  endif
  67   } u32;
  68 };
  69 # elif (LDBL_MANT_DIG == 64) /* 80-bit long double */
  70 static const double scale = 0x1p65;
  71
  72 union ldbl {
  73   long double x;
  74   struct {
  75 #  ifdef __IEEE_LITTLE_ENDIAN /* for Intel CPU */
  76     __uint32_t fracl;
  77     __uint32_t frach;
  78     __uint32_t exp:15;
  79     __uint32_t sign:1;
  80     __uint32_t pad:16;
  81 #  endif
  82 #  ifdef __IEEE_BIG_ENDIAN
  83 #   ifndef ___IEEE_BYTES_LITTLE_ENDIAN /* for m86k */
  84     __uint32_t sign:1;
  85     __uint32_t exp:15;
  86     __uint32_t pad:16;
  87 #   else /* ARM FPA10 math copprocessor */
  88     __uint32_t exp:15;
  89     __uint32_t pad:16;
  90     __uint32_t sign:1;
  91 #   endif
  92     __uint32_t frach;
  93     __uint32_t fracl;
  94 #  endif
  95   } u32;
  96 };
  97 # elif (LDBL_MANT_DIG == 113) /* 128-bit long double */
  98 static const double scale = 0x1p114;
  99
 100 union ldbl {
 101   long double x;
 102   struct {
 103 #  ifdef __IEEE_LITTLE_ENDIAN
 104     __uint32_t fracl;
 105     __uint32_t fraclm;
 106     __uint32_t frachm;
 107     __uint32_t frach:16;
 108     __uint32_t exp:15;
 109     __uint32_t sign:1;
 110 #  endif
 111 #  ifdef __IEEE_BIG_ENDIAN
 112 #   ifndef ___IEEE_BYTES_LITTLE_ENDIAN
 113     __uint32_t sign:1;
 114     __uint32_t exp:15;
 115     __uint32_t frach:16;
 116 #   else /* ARMEL without __VFP_FP__ */
 117     __uint32_t frach:16;
 118     __uint32_t exp:15;
 119     __uint32_t sign:1;
 120 #   endif
 121     __uint32_t frachm;
 122     __uint32_t fraclm;
 123     __uint32_t fracl;
 124 #  endif
 125   } u32;
 126 };
 127 # else
 128 #  error Unsupported long double format.
 129 # endif
 130
 131 static const int scale_exp = LDBL_MANT_DIG + 1;
 132
 133 long double
 134 frexpl (long double x, int *eptr)
 135 {
 136   union ldbl u;
 137   u.x = x;
 138   int e = u.u32.exp;
 139   *eptr = 0;
 140   if (e == (LDBL_MAX_EXP*2 - 1) || x == 0)
 141     return x; /* inf,nan,0 */
 142   if (e == 0) /* subnormal */
 143     {
 144       u.x *= scale;
 145       e = u.u32.exp;
 146       *eptr -= scale_exp;
 147     }
 148   *eptr += e - (LDBL_MAX_EXP - 2);
 149   u.u32.exp = LDBL_MAX_EXP - 2; /* -1 */
 150   return u.x;
 151 }
 152 #endif /* !_LDBL_EQ_DBL */