From 030e94daf50d922278f91070cfa878c39b2ad67d Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Sat, 25 Apr 2015 09:01:11 +0200
Subject: [PATCH] Fourier: fix Bluestein for sequences with more than 46341
 samples (not power-of-two). #286

---
 .../IntegralTransforms/Fourier.Bluestein.cs   | 27 ++++++++++++++++---
 .../MatchingNaiveTransformTest.cs             | 10 +++++++
 2 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
index 020a672a..c008fe04 100644
--- a/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
+++ b/src/Numerics/IntegralTransforms/Fourier.Bluestein.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2014 Math.NET
+// Copyright (c) 2009-2015 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -42,6 +42,12 @@ namespace MathNet.Numerics.IntegralTransforms
     /// </summary>
     public static partial class Fourier
     {
+        /// <summary>
+        /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1
+        /// will cause k*k in the Bluestein sequence to overflow (GH-286).
+        /// </summary>
+        const int BluesteinSequenceLengthThreshold = 46341;
+
         /// <summary>
         /// Generate the bluestein sequence for the provided problem size.
         /// </summary>
@@ -52,10 +58,23 @@ namespace MathNet.Numerics.IntegralTransforms
             double s = Constants.Pi/n;
             var sequence = new Complex[n];
 
-            for (int k = 0; k < sequence.Length; k++)
+            // TODO: benchmark whether the second variation is significantly
+            // faster than the former one. If not just use the former one always.
+            if (n > BluesteinSequenceLengthThreshold)
             {
-                double t = s*(k*k);
-                sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
+                for (int k = 0; k < sequence.Length; k++)
+                {
+                    double t = (s*k)*k;
+                    sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
+                }
+            }
+            else
+            {
+                for (int k = 0; k < sequence.Length; k++)
+                {
+                    double t = s*(k*k);
+                    sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
+                }
             }
 
             return sequence;
diff --git a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs b/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs
index 656c0b64..a531f248 100644
--- a/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs
+++ b/src/UnitTests/IntegralTransformsTests/MatchingNaiveTransformTest.cs
@@ -167,5 +167,15 @@ namespace MathNet.Numerics.UnitTests.IntegralTransformsTests
 
             Verify(samples, 10, options, (a, b) => naive, Fourier.BluesteinForward);
         }
+
+        [Test, Explicit("Long-Running")]
+        public void AlgorithmsMatchNaive_Arbitrary_Large_GH286()
+        {
+            const FourierOptions options = FourierOptions.NoScaling;
+            var samples = Generate.RandomComplex(46500, GetUniform(1));
+            var naive = Fourier.NaiveForward(samples, options);
+
+            Verify(samples, 10, options, (a, b) => naive, Fourier.BluesteinForward);
+        }
     }
 }