>From ef5ffe4324472e80b3c81a8323172d552a2448bd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marcus=20M=C3=BCller?= <address@hidden>
Date: Tue, 27 Jun 2017 11:53:14 +0200
Subject: [PATCH] Fixed floating point math bug (#1348)

https://github.com/gnuradio/gnuradio/issues/1348

We were doing floating point math wrong.
---
 gr-fft/lib/window.cc | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/gr-fft/lib/window.cc b/gr-fft/lib/window.cc
index 126b28978..e80d469e9 100644
--- a/gr-fft/lib/window.cc
+++ b/gr-fft/lib/window.cc
@@ -30,9 +30,10 @@
 namespace gr {
   namespace fft {
 
-#define IzeroEPSILON 1E-21               /* Max error acceptable in Izero */
+#define BESSEL_EPSILON 1E-21               /* Max error acceptable in bessel_i_zero */
 
-    static double Izero(double x)
+    /*! \brief Iterative first kind Bessel function approximation */
+    static double bessel_i_zero(double x)
     {
       double sum, u, halfx, temp;
       int n;
@@ -45,7 +46,7 @@ namespace gr {
         temp *= temp;
         u *= temp;
         sum += u;
-      } while (u >= IzeroEPSILON*sum);
+      } while (u >= BESSEL_EPSILON*sum);
       return(sum);
     }
 
@@ -257,13 +258,17 @@ namespace gr {
 
       std::vector<float> taps(ntaps);
 
-      double IBeta = 1.0/Izero(beta);
+      double bessel_i_beta = 1.0/bessel_i_zero(beta);
       double inm1 = 1.0/((double)(ntaps-1));
       double temp;
 
-      for(int i = 0; i < ntaps; i++) {
+      /* dragging taps[0] out of the loop, because ((double)(-1))*((double)(-1))
+      might be 1.0 + epsilon and then we get sqrt() of something < 0, i.e. NaN.
+      https://github.com/gnuradio/gnuradio/issues/1348 */
+      taps[0] = bessel_i_zero(0) * bessel_i_beta;
+      for(int i = 1; i < ntaps; i++) {
         temp = 2*i*inm1 - 1;
-        taps[i] = Izero(beta*sqrt(1.0-temp*temp)) * IBeta;
+        taps[i] = bessel_i_zero(beta*sqrt(1.0-temp*temp)) * bessel_i_beta;
       }
       return taps;
     }
-- 
2.13.0