libm: Handle 32-bit double/64-bit long double

There were lots of assumptions that doubles were 64-bits remaining in
the library. Handle cases where float and double are 32-bits while
long double is 64-bits.

Signed-off-by: Keith Packard <[email protected]>

Author: keith-packard

newlib/libc/include/machine/ieeefp.h

@@ -89,6 +89,14 @@ warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 # endif
 #endif
 
+#if __SIZEOF_DOUBLE__ == __SIZEOF_FLOAT__
+# define _DBL_EQ_FLT
+#endif
+
+#if __SIZEOF_LONG_DOUBLE == __SIZEOF_DOUBLE__
+# define LDBL_EQ_DBL
+#endif
+
 #if __SIZEOF_DOUBLE__ == 4
 # define _DOUBLE_IS_32BITS
 #endif
@@ -97,14 +105,6 @@ warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 # define _LONG_DOUBLE_IS_32BITS
 #endif
 
-#if defined(__SIZEOF_FLOAT__) && defined(__SIZEOF_DOUBLE__)
-# if __SIZEOF_FLOAT__ == __SIZEOF_DOUBLE__
-#  define _DBL_EQ_FLT
-# else
-#  undef _DBL_EQ_FLT
-#endif
-#endif
-
 #if (defined(__arm__) || defined(__thumb__)) && !defined(__MAVERICK__)
 /* arm with hard fp and soft dp cannot use new float code */
 # if (__ARM_FP & 4) && !(__ARM_FP & 8)
@@ -230,13 +230,14 @@ warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 
 #ifdef __sh__
-#ifdef __LITTLE_ENDIAN__
+#define _IEEE_754_2008_SNAN 0
+#if __FLOAT_WORD_ORDER__ == __ORDER_LITTLE_ENDIAN__
 #define __IEEE_LITTLE_ENDIAN
-#else
+#elif __FLOAT_WORD_ORDER__ == __ORDER_BIG_ENDIAN__
 #define __IEEE_BIG_ENDIAN
 #endif
-#if defined(__SH2E__) || defined(__SH3E__) || defined(__SH4_SINGLE_ONLY__) || defined(__SH2A_SINGLE_ONLY__)
-#define _DOUBLE_IS_32BITS
+#ifdef __SH_FPU_ANY__
+#define _SUPPORTS_ERREXCEPT
 #endif
 #endif
 

newlib/libc/include/math.h

@@ -350,11 +350,7 @@ int __issignalingf(float f);
 int __issignaling(double d);
 
 #if defined(_HAVE_LONG_DOUBLE)
-#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__ && __SIZEOF_DOUBLE__ > 0
-static inline int __issignalingl(long double d) { return __issignaling((double) d); }
-#else
 int __issignalingl(long double d);
-#endif
 #define issignaling(__x)						\
 	((sizeof(__x) == sizeof(float))  ? __issignalingf(__x) :	\
 	 (sizeof(__x) == sizeof(double)) ? __issignaling ((double) (__x)) :	\

newlib/libm/common/exp.c

@@ -61,13 +61,13 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
     {
       /* k > 0, the exponent of scale might have overflowed by <= 460.  */
       sbits -= 1009ull << 52;
-      scale = asdouble (sbits);
+      scale = asfloat64 (sbits);
       y = 0x1p1009 * (scale + scale * tmp);
       return check_oflow (y);
     }
   /* k < 0, need special care in the subnormal range.  */
   sbits += 1022ull << 52;
-  scale = asdouble (sbits);
+  scale = asfloat64 (sbits);
   y = scale + scale * tmp;
   if (y < 1.0)
     {
@@ -148,7 +148,7 @@ exp (double x)
   /* 2^(k/N) ~= scale * (1 + tail).  */
   idx = 2 * (ki % N);
   top = ki << (52 - EXP_TABLE_BITS);
-  tail = asdouble (T[idx]);
+  tail = asfloat64 (T[idx]);
   /* This is only a valid scale when -1023*N < k < 1024*N.  */
   sbits = T[idx + 1] + top;
   /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
@@ -165,7 +165,7 @@ exp (double x)
 #endif
   if (unlikely (abstop == 0))
     return specialcase (tmp, sbits, ki);
-  scale = asdouble (sbits);
+  scale = asfloat64 (sbits);
   /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-739, so there
      is no spurious underflow here even without fma.  */
   return scale + scale * tmp;

newlib/libm/common/exp2.c

@@ -59,13 +59,13 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
     {
       /* k > 0, the exponent of scale might have overflowed by 1.  */
       sbits -= 1ull << 52;
-      scale = asdouble (sbits);
+      scale = asfloat64 (sbits);
       y = 2 * (scale + scale * tmp);
       return check_oflow (y);
     }
   /* k < 0, need special care in the subnormal range.  */
   sbits += 1022ull << 52;
-  scale = asdouble (sbits);
+  scale = asfloat64 (sbits);
   y = scale + scale * tmp;
   if (y < 1.0)
     {
@@ -135,7 +135,7 @@ exp2 (double x)
   /* 2^(k/N) ~= scale * (1 + tail).  */
   idx = 2 * (ki % N);
   top = ki << (52 - EXP_TABLE_BITS);
-  tail = asdouble (T[idx]);
+  tail = asfloat64 (T[idx]);
   /* This is only a valid scale when -1023*N < k < 1024*N.  */
   sbits = T[idx + 1] + top;
   /* exp2(x) = 2^(k/N) * 2^r ~= scale + scale * (tail + 2^r - 1).  */
@@ -152,7 +152,7 @@ exp2 (double x)
 #endif
   if (unlikely (abstop == 0))
     return specialcase (tmp, sbits, ki);
-  scale = asdouble (sbits);
+  scale = asfloat64 (sbits);
   /* Note: tmp == 0 or |tmp| > 2^-65 and scale > 2^-928, so there
      is no spurious underflow here even without fma.  */
   return scale + scale * tmp;

newlib/libm/common/fdlibm.h

@@ -168,13 +168,15 @@
 
 #define X_TLOSS		1.41484755040568800000e+16 
 
+#ifdef _NEED_FLOAT64
 extern __int32_t __rem_pio2 (__float64,__float64*);
 
 /* fdlibm kernel function */
 extern __float64 __kernel_sin (__float64,__float64,int);
 extern __float64 __kernel_cos (__float64,__float64);
 extern __float64 __kernel_tan (__float64,__float64,int);
 extern int    __kernel_rem_pio2 (__float64*,__float64*,int,int,int,const __int32_t*);
+#endif
 
 extern __int32_t __rem_pio2f (float,float*);
 
@@ -267,7 +269,7 @@ do {								\
   ieee_double_shape_type iw_u;					\
   iw_u.parts.msw = (ix0);					\
   iw_u.parts.lsw = (ix1);					\
-  (d) = asdouble(iw_u.bits);					\
+  (d) = asfloat64(iw_u.bits);                                   \
 } while (0)
 
 /* Set the more significant 32 bits of a double from an int.  */
@@ -277,7 +279,7 @@ do {								\
   ieee_double_shape_type sh_u;					\
   sh_u.bits = asuint64(d);					\
   sh_u.parts.msw = (v);						\
-  (d) = asdouble(sh_u.bits);					\
+  (d) = asfloat64(sh_u.bits);                                   \
 } while (0)
 
 /* Set the less significant 32 bits of a double from an int.  */
@@ -287,7 +289,7 @@ do {								\
   ieee_double_shape_type sl_u;					\
   sl_u.bits = asuint64(d);					\
   sl_u.parts.lsw = (v);						\
-  (d) = asdouble(sl_u.bits);					\
+  (d) = asfloat64(sl_u.bits);                                   \
 } while (0)
 
 /* A union which permits us to convert between a float and a 32 bit

newlib/libm/common/log.c

@@ -142,7 +142,7 @@ log (double x)
   iz = ix - (tmp & 0xfffULL << 52);
   invc = T[i].invc;
   logc = T[i].logc;
-  z = asdouble (iz);
+  z = asfloat64 (iz);
 
   /* log(x) = log1p(z/c-1) + log(c) + k*Ln2.  */
   /* r ~= z/c - 1, |r| < 1/(2*N).  */

newlib/libm/common/log2.c

@@ -77,7 +77,7 @@ double
       lo = r * InvLn2lo + fma (r, InvLn2hi, -hi);
 #else
       double_t rhi, rlo;
-      rhi = asdouble (asuint64 (r) & -1ULL << 32);
+      rhi = asfloat64 (asuint64 (r) & -1ULL << 32);
       rlo = r - rhi;
       hi = rhi * InvLn2hi;
       lo = rlo * InvLn2hi + r * InvLn2lo;
@@ -118,7 +118,7 @@ double
   iz = ix - (tmp & 0xfffULL << 52);
   invc = T[i].invc;
   logc = T[i].logc;
-  z = asdouble (iz);
+  z = asfloat64 (iz);
   kd = (double_t) k;
 
   /* log2(x) = log2(z/c) + log2(c) + k.  */
@@ -132,7 +132,7 @@ double
   double_t rhi, rlo;
   /* rounding error: 0x1p-55/N + 0x1p-65.  */
   r = (z - T2[i].chi - T2[i].clo) * invc;
-  rhi = asdouble (asuint64 (r) & -1ULL << 32);
+  rhi = asfloat64 (asuint64 (r) & -1ULL << 32);
   rlo = r - rhi;
   t1 = rhi * InvLn2hi;
   t2 = rlo * InvLn2hi + r * InvLn2lo;