Back to index

glibc  2.9
ldbl2mpn.c
Go to the documentation of this file.
00001 /* Copyright (C) 1995, 1996, 1997, 2000, 2007 Free Software Foundation, Inc.
00002    This file is part of the GNU C Library.
00003 
00004    The GNU C Library is free software; you can redistribute it and/or
00005    modify it under the terms of the GNU Lesser General Public
00006    License as published by the Free Software Foundation; either
00007    version 2.1 of the License, or (at your option) any later version.
00008 
00009    The GNU C Library is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012    Lesser General Public License for more details.
00013 
00014    You should have received a copy of the GNU Lesser General Public
00015    License along with the GNU C Library; if not, write to the Free
00016    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00017    02111-1307 USA.  */
00018 
00019 #include "gmp.h"
00020 #include "gmp-impl.h"
00021 #include "longlong.h"
00022 #include <ieee754.h>
00023 #include <float.h>
00024 #include <stdlib.h>
00025 
00026 /* Convert a `long double' in IEEE854 standard double-precision format to a
00027    multi-precision integer representing the significand scaled up by its
00028    number of bits (64 for long double) and an integral power of two
00029    (MPN frexpl). */
00030 
00031 mp_size_t
00032 __mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size,
00033                         int *expt, int *is_neg,
00034                         long double value)
00035 {
00036   union ieee854_long_double u;
00037   u.d = value;
00038 
00039   *is_neg = u.ieee.negative;
00040   *expt = (int) u.ieee.exponent - IEEE854_LONG_DOUBLE_BIAS;
00041 
00042 #if BITS_PER_MP_LIMB == 32
00043   res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction.  */
00044   res_ptr[1] = u.ieee.mantissa0; /* High-order 32 bits.  */
00045   #define N 2
00046 #elif BITS_PER_MP_LIMB == 64
00047   /* Hopefully the compiler will combine the two bitfield extracts
00048      and this composition into just the original quadword extract.  */
00049   res_ptr[0] = ((mp_limb_t) u.ieee.mantissa0 << 32) | u.ieee.mantissa1;
00050   #define N 1
00051 #else
00052   #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
00053 #endif
00054 
00055   if (u.ieee.exponent == 0)
00056     {
00057       /* A biased exponent of zero is a special case.
00058         Either it is a zero or it is a denormal number.  */
00059       if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2.  */
00060        /* It's zero.  */
00061        *expt = 0;
00062       else
00063        {
00064          /* It is a denormal number, meaning it has no implicit leading
00065             one bit, and its exponent is in fact the format minimum.  */
00066          int cnt;
00067 
00068          /* One problem with Intel's 80-bit format is that the explicit
00069             leading one in the normalized representation has to be zero
00070             for denormalized number.  If it is one, the number is according
00071             to Intel's specification an invalid number.  We make the
00072             representation unique by explicitly clearing this bit.  */
00073          res_ptr[N - 1] &= ~(1L << ((LDBL_MANT_DIG - 1) % BITS_PER_MP_LIMB));
00074 
00075          if (res_ptr[N - 1] != 0)
00076            {
00077              count_leading_zeros (cnt, res_ptr[N - 1]);
00078              if (cnt != 0)
00079               {
00080 #if N == 2
00081                 res_ptr[N - 1] = res_ptr[N - 1] << cnt
00082                                | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt));
00083                 res_ptr[0] <<= cnt;
00084 #else
00085                 res_ptr[N - 1] <<= cnt;
00086 #endif
00087               }
00088              *expt = LDBL_MIN_EXP - 1 - cnt;
00089            }
00090          else if (res_ptr[0] != 0)
00091            {
00092              count_leading_zeros (cnt, res_ptr[0]);
00093              res_ptr[N - 1] = res_ptr[0] << cnt;
00094              res_ptr[0] = 0;
00095              *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt;
00096            }
00097          else
00098            {
00099              /* This is the special case of the pseudo denormal number
00100                with only the implicit leading bit set.  The value is
00101                in fact a normal number and so we have to treat this
00102                case differently.  */
00103 #if N == 2
00104              res_ptr[N - 1] = 0x80000000ul;
00105 #else
00106              res_ptr[0] = 0x8000000000000000ul;
00107 #endif
00108              *expt = LDBL_MIN_EXP - 1;
00109            }
00110        }
00111     }
00112   else if (u.ieee.exponent < 0x7fff
00113 #if N == 2
00114           && res_ptr[0] == 0
00115 #endif
00116           && res_ptr[N - 1] == 0)
00117     /* Pseudo zero.  */
00118     *expt = 0;
00119 
00120   return N;
00121 }