This change should significantly improve the performance of single precision floating point math library functions. The vast majority of changes have to do with preventing the compiler from needlessly promoting floats to doubles, performing the calculation with doubles, only to demote the result to float. These changes only affect the math lib functions that return float.

include/nuttx/math.h

@@ -93,10 +93,15 @@
 #define NAN         (0.0/0.0)
 #define HUGE_VAL    INFINITY
 
+#define INFINITY_F  (1.0F/0.0F)
+#define NAN_F       (0.0F/0.0F)
+
 #define isnan(x)    ((x) != (x))
 #define isinf(x)    (((x) == INFINITY) || ((x) == -INFINITY))
 #define isfinite(x) (!(isinf(x)) && (x != NAN))
 
+#define isinf_f(x)  (((x) == INFINITY_F) || ((x) == -INFINITY_F))
+
 /* Exponential and Logarithmic constants ************************************/
 
 #define M_E        2.7182818284590452353602874713526625
@@ -116,6 +121,9 @@
 #define M_2_PI     0.6366197723675813430755350534900574
 #define M_2_SQRTPI 1.1283791670955125738961589031215452
 
+#define M_PI_F     ((float)M_PI)
+#define M_PI_2_F   ((float)M_PI_2)
+
 /****************************************************************************
  * Public Function Prototypes
  ****************************************************************************/

libc/libc.h

@@ -1,7 +1,7 @@
 /****************************************************************************
  * libc/libc.h
  *
- *   Copyright (C) 2007-2014 Gregory Nutt. All rights reserved.
+ *   Copyright (C) 2007-2014, 2016 Gregory Nutt. All rights reserved.
  *   Author: Gregory Nutt <gnutt@nuttx.org>
  *
  * Redistribution and use in source and binary forms, with or without
@@ -211,6 +211,7 @@ int lib_checkbase(int base, const char **pptr);
 /* Defined in lib_expi.c */
 
 #ifdef CONFIG_LIBM
+float  lib_expif(size_t n);
 double lib_expi(size_t n);
 #endif
 

libc/math/Make.defs

@@ -59,6 +59,7 @@ CSRCS += lib_tanhl.c lib_asinhl.c lib_acoshl.c lib_atanhl.c lib_erfl.c lib_copys
 CSRCS += lib_truncl.c
 
 CSRCS += lib_libexpi.c lib_libsqrtapprox.c
+CSRCS += lib_libexpif.c
 
 # Add the floating point math directory to the build
 

libc/math/lib_acosf.c

@@ -37,5 +37,5 @@
 
 float acosf(float x)
 {
-  return (M_PI_2 - asinf(x));
+  return (M_PI_2_F - asinf(x));
 }

libc/math/lib_acoshf.c

@@ -48,5 +48,5 @@
 
 float acoshf(float x)
 {
-  return logf(x + sqrtf(x * x - 1));
+  return logf(x + sqrtf(x * x - 1.0F));
 }

libc/math/lib_asinf.c

@@ -47,9 +47,9 @@ float asinf(float x)
       y_sin = sinf(y);
       y_cos = cosf(y);
 
-      if (y > M_PI_2 || y < -M_PI_2)
+      if (y > M_PI_2_F || y < -M_PI_2_F)
         {
-          y = fmodf(y, M_PI);
+          y = fmodf(y, M_PI_F);
         }
 
       if (y_sin + FLT_EPSILON >= x && y_sin - FLT_EPSILON <= x)

libc/math/lib_asinhf.c

@@ -48,5 +48,5 @@
 
 float asinhf(float x)
 {
-  return logf(x + sqrtf(x * x + 1));
+  return logf(x + sqrtf(x * x + 1.0F));
 }

libc/math/lib_atan2f.c

@@ -43,22 +43,22 @@ float atan2f(float y, float x)
     }
   else if (y >= 0 && x < 0)
     {
-      return atanf(y / x) + M_PI;
+      return atanf(y / x) + M_PI_F;
     }
   else if (y < 0)
     {
       if (x == 0)
         {
-          return -M_PI_2;
+          return -M_PI_2_F;
         }
       else /* Can only be x < 0 */
         {
-          return atanf(y / x) - M_PI;
+          return atanf(y / x) - M_PI_F;
         }
     }
   else if (y > 0 && x == 0)
     {
-      return M_PI_2;
+      return M_PI_2_F;
     }
   else /* if (y == 0 && x == 0) Undefined but returns normally 0 */
     {

libc/math/lib_atanf.c

@@ -39,5 +39,5 @@
 
 float atanf(float x)
 {
-  return asinf(x / sqrtf(x * x + 1));
+  return asinf(x / sqrtf(x * x + 1.0F));
 }

libc/math/lib_atanhf.c

@@ -48,5 +48,5 @@
 
 float atanhf(float x)
 {
-  return 0.5 * logf((1 + x) / (1 - x));
+  return 0.5F * logf((1.0F + x) / (1.0F - x));
 }

libc/math/lib_ceilf.c

@@ -38,9 +38,9 @@
 float ceilf(float x)
 {
   modff(x, &x);
-  if (x > 0.0)
+  if (x > 0.0F)
     {
-      x += 1.0;
+      x += 1.0F;
     }
 
   return x;

libc/math/lib_cosf.c

@@ -37,5 +37,5 @@
 
 float cosf(float x)
 {
-  return sinf(x + M_PI_2);
+  return sinf(x + M_PI_2_F);
 }

libc/math/lib_coshf.c

@@ -38,5 +38,5 @@
 float coshf(float x)
 {
   x = expf(x);
-  return ((x + (1.0 / x)) / 2.0);
+  return ((x + (1.0F / x)) / 2.0F);
 }

libc/math/lib_erff.c

@@ -57,14 +57,14 @@ float erff(float x)
   float t;
   float a1, a2, a3, a4, a5, p;
 
-  a1 =  0.254829592;
-  a2 = -0.284496736;
-  a3 =  1.421413741;
-  a4 = -1.453152027;
-  a5 =  1.061405429;
-  p  =  0.3275911;
+  a1 =  0.254829592F;
+  a2 = -0.284496736F;
+  a3 =  1.421413741F;
+  a4 = -1.453152027F;
+  a5 =  1.061405429F;
+  p  =  0.3275911F;
 
   sign = (x >= 0 ? 1 : -1);
-  t = 1.0/(1.0 + p*x);
-  return sign * (1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * expf(-x * x));
+  t = 1.0F/(1.0F + p*x);
+  return sign * (1.0F - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * expf(-x * x));
 }

libc/math/lib_expf.c

@@ -89,8 +89,8 @@ float expf(float x)
 
   /* Perform Taylor series approximation with eleven terms */
 
-  value = 0.0;
-  x0 = 1.0;
+  value = 0.0F;
+  x0 = 1.0F;
   for (i = 0; i < 10; i++)
     {
       value += x0 * _flt_inv_fact[i];
@@ -99,11 +99,11 @@ float expf(float x)
 
   /* Multiply by exp of the integer component */
 
-  value *= lib_expi(int_part);
+  value *= lib_expif(int_part);
 
   if (invert)
     {
-      return (1.0 / value);
+      return (1.0F / value);
     }
   else
     {

libc/math/lib_floorf.c

@@ -38,9 +38,9 @@
 float floorf(float x)
 {
   modff(x, &x);
-  if (x < 0.0)
+  if (x < 0.0F)
     {
-      x -= 1.0;
+      x -= 1.0F;
     }
 
   return x;

libc/math/lib_frexpf.c

@@ -38,5 +38,5 @@
 float frexpf(float x, int *exponent)
 {
   *exponent = (int)ceilf(log2f(x));
-  return x / ldexpf(1.0, *exponent);
+  return x / ldexpf(1.0F, *exponent);
 }

libc/math/lib_ldexpf.c

@@ -37,5 +37,5 @@
 
 float ldexpf(float x, int n)
 {
-  return (x * powf(2.0, (float)n));
+  return (x * powf(2.0F, (float)n));
 }

libc/math/lib_log10f.c

@@ -37,5 +37,5 @@
 
 float log10f(float x)
 {
-  return (logf(x) / M_LN10);
+  return (logf(x) / (float)M_LN10);
 }

libc/math/lib_log2f.c

@@ -37,5 +37,5 @@
 
 float log2f(float x)
 {
-  return (logf(x) / M_LN2);
+  return (logf(x) / (float)M_LN2);
 }

libc/math/lib_logf.c

@@ -40,8 +40,8 @@ float logf(float x)
 {
   float y, y_old, ey, epsilon;
 
-  y = 0.0;
-  y_old = 1.0;
+  y = 0.0F;
+  y_old = 1.0F;
   epsilon = FLT_EPSILON;
 
   while (y > y_old + epsilon || y < y_old - epsilon)
@@ -50,25 +50,25 @@ float logf(float x)
       ey = exp(y);
       y -= (ey - x) / ey;
 
-      if (y > 700.0)
+      if (y > 700.0F)
         {
-          y = 700.0;
+          y = 700.0F;
         }
 
-      if (y < -700.0)
+      if (y < -700.0F)
         {
-          y = -700.0;
+          y = -700.0F;
         }
 
-      epsilon = (fabs(y) > 1.0) ? fabs(y) * FLT_EPSILON : FLT_EPSILON;
+      epsilon = (fabsf(y) > 1.0F) ? fabsf(y) * FLT_EPSILON : FLT_EPSILON;
     }
 
-  if (y == 700.0)
+  if (y == 700.0F)
     {
       return INFINITY;
     }
 
-  if (y == -700.0)
+  if (y == -700.0F)
     {
       return INFINITY;
     }

libc/math/lib_modff.c

@@ -37,14 +37,14 @@
 
 float modff(float x, float *iptr)
 {
-  if (fabsf(x) >= 8388608.0)
+  if (fabsf(x) >= 8388608.0F)
     {
       *iptr = x;
-      return 0.0;
+      return 0.0F;
     }
-  else if (fabs(x) < 1.0)
+  else if (fabsf(x) < 1.0F)
     {
-      *iptr = 0.0;
+      *iptr = 0.0F;
       return x;
     }
   else

libc/math/lib_rintf.c

@@ -105,15 +105,15 @@ float rintf(float x)
   linteger  = (long)x;
   fremainder = x - (float)linteger;
 
-  if (x < 0.0)
+  if (x < 0.0F)
     {
       /* fremainder should be in range 0 .. -1 */
 
-      if (fremainder == -0.5)
+      if (fremainder == -0.5F)
         {
           linteger = ((linteger + 1) & ~1);
         }
-      else if (fremainder < -0.5)
+      else if (fremainder < -0.5F)
         {
           linteger--;
         }
@@ -122,11 +122,11 @@ float rintf(float x)
     {
       /* fremainder should be in range 0 .. 1 */
 
-      if (fremainder == 0.5)
+      if (fremainder == 0.5F)
         {
           linteger = ((linteger + 1) & ~1);
         }
-      else if (fremainder > 0.5)
+      else if (fremainder > 0.5F)
         {
           linteger++;
         }

libc/math/lib_sinf.c

@@ -58,31 +58,31 @@ float sinf(float x)
 
   /* Move x to [-pi, pi) */
 
-  x = fmodf(x, 2 * M_PI);
-  if (x >= M_PI)
+  x = fmodf(x, 2 * M_PI_F);
+  if (x >= M_PI_F)
     {
-      x -= 2 * M_PI;
+      x -= 2 * M_PI_F;
     }
 
-  if (x < -M_PI)
+  if (x < -M_PI_F)
     {
-      x += 2 * M_PI;
+      x += 2 * M_PI_F;
     }
 
   /* Move x to [-pi/2, pi/2) */
 
-  if (x >= M_PI_2)
+  if (x >= M_PI_2_F)
     {
-      x = M_PI - x;
+      x = M_PI_F - x;
     }
 
-  if (x < -M_PI_2)
+  if (x < -M_PI_2_F)
     {
-      x = -M_PI - x;
+      x = -M_PI_F - x;
     }
 
   x_squared = x * x;
-  sin_x = 0.0;
+  sin_x = 0.0F;
 
   /* Perform Taylor series approximation for sin(x) with six terms */
 

libc/math/lib_sinhf.c

@@ -38,5 +38,5 @@
 float sinhf(float x)
 {
   x = expf(x);
-  return ((x - (1.0 / x)) / 2.0);
+  return ((x - (1.0F / x)) / 2.0F);
 }

libc/math/lib_sqrtf.c

@@ -47,25 +47,25 @@ float sqrtf(float x)
 
   /* Filter out invalid/trivial inputs */
 
-  if (x < 0.0)
+  if (x < 0.0F)
     {
       set_errno(EDOM);
-      return NAN;
+      return NAN_F;
     }
 
   if (isnan(x))
     {
-      return NAN;
+      return NAN_F;
     }
 
-  if (isinf(x))
+  if (isinf_f(x))
     {
-      return INFINITY;
+      return INFINITY_F;
     }
 
-  if (x == 0.0)
+  if (x == 0.0F)
     {
-      return 0.0;
+      return 0.0F;
     }
 
   /* Guess square root (using bit manipulation) */
@@ -76,9 +76,9 @@ float sqrtf(float x)
    * definitely optimal
    */
 
-  y = 0.5 * (y + x / y);
-  y = 0.5 * (y + x / y);
-  y = 0.5 * (y + x / y);
+  y = 0.5F * (y + x / y);
+  y = 0.5F * (y + x / y);
+  y = 0.5F * (y + x / y);
 
   return y;
 }

libc/math/lib_tanhf.c

@@ -38,7 +38,7 @@
 float tanhf(float x)
 {
   float x0 = expf(x);
-  float x1 = 1.0 / x0;
+  float x1 = 1.0F / x0;
 
   return ((x0 + x1) / (x0 - x1));
 }