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
///
public static partial class Fourier
{
+ ///
+ /// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1
+ /// will cause k*k in the Bluestein sequence to overflow (GH-286).
+ ///
+ const int BluesteinSequenceLengthThreshold = 46341;
+
///
/// Generate the bluestein sequence for the provided problem size.
///
@@ -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);
+ }
}
}