Add a much better sinc() implementation.

diff --git a/qscale.c b/qscale.c
index 6d529bcb4d1d73c49b314eb036b54915fc7dd258..1b3ab57c0b9bd20a76785d66caac745b445a774c 100644 (file)
--- a/qscale.c
+++ b/qscale.c
@@ -9,11 +9,18 @@
 
 double sinc(double x)
 {
-       // This is bad for very small x, should use power series instead.
-       if (x == 0.0)
-               return 1.0;
-       else
+       static const double eps = 2.22045e-016;
+       static const double cutoff = sqrt(sqrt(eps));
+
+       if (abs(x) < cutoff) {
+               // For small |x|, use Taylor series instead
+               const double x2 = x * x;
+               const double x4 = x2 * x2;
+
+               return 1.0 - x2 / 6.0 + x4 / 120.0;
+       } else {
                return sin(x) / x;
+       }
 }
 
 double lanczos_tap(double x)