Home | History | Annotate | Download | only in lib
      1 /* Test for NaN that does not need libm.
      2    Copyright (C) 2007-2012 Free Software Foundation, Inc.
      3 
      4    This program is free software: you can redistribute it and/or modify
      5    it under the terms of the GNU General Public License as published by
      6    the Free Software Foundation; either version 3 of the License, or
      7    (at your option) any later version.
      8 
      9    This program 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
     12    GNU General Public License for more details.
     13 
     14    You should have received a copy of the GNU General Public License
     15    along with this program.  If not, see <http://www.gnu.org/licenses/>.  */
     16 
     17 /* Written by Bruno Haible <bruno (at) clisp.org>, 2007.  */
     18 
     19 #include <config.h>
     20 
     21 /* Specification.  */
     22 #ifdef USE_LONG_DOUBLE
     23 /* Specification found in math.h or isnanl-nolibm.h.  */
     24 extern int rpl_isnanl (long double x) _GL_ATTRIBUTE_CONST;
     25 #elif ! defined USE_FLOAT
     26 /* Specification found in math.h or isnand-nolibm.h.  */
     27 extern int rpl_isnand (double x);
     28 #else /* defined USE_FLOAT */
     29 /* Specification found in math.h or isnanf-nolibm.h.  */
     30 extern int rpl_isnanf (float x);
     31 #endif
     32 
     33 #include <float.h>
     34 #include <string.h>
     35 
     36 #include "float+.h"
     37 
     38 #ifdef USE_LONG_DOUBLE
     39 # define FUNC rpl_isnanl
     40 # define DOUBLE long double
     41 # define MAX_EXP LDBL_MAX_EXP
     42 # define MIN_EXP LDBL_MIN_EXP
     43 # if defined LDBL_EXPBIT0_WORD && defined LDBL_EXPBIT0_BIT
     44 #  define KNOWN_EXPBIT0_LOCATION
     45 #  define EXPBIT0_WORD LDBL_EXPBIT0_WORD
     46 #  define EXPBIT0_BIT LDBL_EXPBIT0_BIT
     47 # endif
     48 # define SIZE SIZEOF_LDBL
     49 # define L_(literal) literal##L
     50 #elif ! defined USE_FLOAT
     51 # define FUNC rpl_isnand
     52 # define DOUBLE double
     53 # define MAX_EXP DBL_MAX_EXP
     54 # define MIN_EXP DBL_MIN_EXP
     55 # if defined DBL_EXPBIT0_WORD && defined DBL_EXPBIT0_BIT
     56 #  define KNOWN_EXPBIT0_LOCATION
     57 #  define EXPBIT0_WORD DBL_EXPBIT0_WORD
     58 #  define EXPBIT0_BIT DBL_EXPBIT0_BIT
     59 # endif
     60 # define SIZE SIZEOF_DBL
     61 # define L_(literal) literal
     62 #else /* defined USE_FLOAT */
     63 # define FUNC rpl_isnanf
     64 # define DOUBLE float
     65 # define MAX_EXP FLT_MAX_EXP
     66 # define MIN_EXP FLT_MIN_EXP
     67 # if defined FLT_EXPBIT0_WORD && defined FLT_EXPBIT0_BIT
     68 #  define KNOWN_EXPBIT0_LOCATION
     69 #  define EXPBIT0_WORD FLT_EXPBIT0_WORD
     70 #  define EXPBIT0_BIT FLT_EXPBIT0_BIT
     71 # endif
     72 # define SIZE SIZEOF_FLT
     73 # define L_(literal) literal##f
     74 #endif
     75 
     76 #define EXP_MASK ((MAX_EXP - MIN_EXP) | 7)
     77 
     78 #define NWORDS \
     79   ((sizeof (DOUBLE) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
     80 typedef union { DOUBLE value; unsigned int word[NWORDS]; } memory_double;
     81 
     82 int
     83 FUNC (DOUBLE x)
     84 {
     85 #ifdef KNOWN_EXPBIT0_LOCATION
     86 # if defined USE_LONG_DOUBLE && ((defined __ia64 && LDBL_MANT_DIG == 64) || (defined __x86_64__ || defined __amd64__) || (defined __i386 || defined __i386__ || defined _I386 || defined _M_IX86 || defined _X86_)) && !HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
     87   /* Special CPU dependent code is needed to treat bit patterns outside the
     88      IEEE 754 specification (such as Pseudo-NaNs, Pseudo-Infinities,
     89      Pseudo-Zeroes, Unnormalized Numbers, and Pseudo-Denormals) as NaNs.
     90      These bit patterns are:
     91        - exponent = 0x0001..0x7FFF, mantissa bit 63 = 0,
     92        - exponent = 0x0000, mantissa bit 63 = 1.
     93      The NaN bit pattern is:
     94        - exponent = 0x7FFF, mantissa >= 0x8000000000000001.  */
     95   memory_double m;
     96   unsigned int exponent;
     97 
     98   m.value = x;
     99   exponent = (m.word[EXPBIT0_WORD] >> EXPBIT0_BIT) & EXP_MASK;
    100 #  ifdef WORDS_BIGENDIAN
    101   /* Big endian: EXPBIT0_WORD = 0, EXPBIT0_BIT = 16.  */
    102   if (exponent == 0)
    103     return 1 & (m.word[0] >> 15);
    104   else if (exponent == EXP_MASK)
    105     return (((m.word[0] ^ 0x8000U) << 16) | m.word[1] | (m.word[2] >> 16)) != 0;
    106   else
    107     return 1 & ~(m.word[0] >> 15);
    108 #  else
    109   /* Little endian: EXPBIT0_WORD = 2, EXPBIT0_BIT = 0.  */
    110   if (exponent == 0)
    111     return (m.word[1] >> 31);
    112   else if (exponent == EXP_MASK)
    113     return ((m.word[1] ^ 0x80000000U) | m.word[0]) != 0;
    114   else
    115     return (m.word[1] >> 31) ^ 1;
    116 #  endif
    117 # else
    118   /* Be careful to not do any floating-point operation on x, such as x == x,
    119      because x may be a signaling NaN.  */
    120 #  if defined __SUNPRO_C || defined __ICC || defined _MSC_VER \
    121       || defined __DECC || defined __TINYC__ \
    122       || (defined __sgi && !defined __GNUC__)
    123   /* The Sun C 5.0, Intel ICC 10.0, Microsoft Visual C/C++ 9.0, Compaq (ex-DEC)
    124      6.4, and TinyCC compilers don't recognize the initializers as constant
    125      expressions.  The Compaq compiler also fails when constant-folding
    126      0.0 / 0.0 even when constant-folding is not required.  The Microsoft
    127      Visual C/C++ compiler also fails when constant-folding 1.0 / 0.0 even
    128      when constant-folding is not required. The SGI MIPSpro C compiler
    129      complains about "floating-point operation result is out of range".  */
    130   static DOUBLE zero = L_(0.0);
    131   memory_double nan;
    132   DOUBLE plus_inf = L_(1.0) / zero;
    133   DOUBLE minus_inf = -L_(1.0) / zero;
    134   nan.value = zero / zero;
    135 #  else
    136   static memory_double nan = { L_(0.0) / L_(0.0) };
    137   static DOUBLE plus_inf = L_(1.0) / L_(0.0);
    138   static DOUBLE minus_inf = -L_(1.0) / L_(0.0);
    139 #  endif
    140   {
    141     memory_double m;
    142 
    143     /* A NaN can be recognized through its exponent.  But exclude +Infinity and
    144        -Infinity, which have the same exponent.  */
    145     m.value = x;
    146     if (((m.word[EXPBIT0_WORD] ^ nan.word[EXPBIT0_WORD])
    147          & (EXP_MASK << EXPBIT0_BIT))
    148         == 0)
    149       return (memcmp (&m.value, &plus_inf, SIZE) != 0
    150               && memcmp (&m.value, &minus_inf, SIZE) != 0);
    151     else
    152       return 0;
    153   }
    154 # endif
    155 #else
    156   /* The configuration did not find sufficient information.  Give up about
    157      the signaling NaNs, handle only the quiet NaNs.  */
    158   if (x == x)
    159     {
    160 # if defined USE_LONG_DOUBLE && ((defined __ia64 && LDBL_MANT_DIG == 64) || (defined __x86_64__ || defined __amd64__) || (defined __i386 || defined __i386__ || defined _I386 || defined _M_IX86 || defined _X86_)) && !HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
    161       /* Detect any special bit patterns that pass ==; see comment above.  */
    162       memory_double m1;
    163       memory_double m2;
    164 
    165       memset (&m1.value, 0, SIZE);
    166       memset (&m2.value, 0, SIZE);
    167       m1.value = x;
    168       m2.value = x + (x ? 0.0L : -0.0L);
    169       if (memcmp (&m1.value, &m2.value, SIZE) != 0)
    170         return 1;
    171 # endif
    172       return 0;
    173     }
    174   else
    175     return 1;
    176 #endif
    177 }
    178