github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/dsp/fourier/internal/fftpack/fftpack_test.go (about) 1 // Copyright ©2018 The Gonum Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 // This is a translation and extension of the FFTPACK test 6 // functions by Paul N Swarztrauber, placed in the public 7 // domain at http://www.netlib.org/fftpack/. 8 9 package fftpack 10 11 import ( 12 "math" 13 "math/cmplx" 14 "reflect" 15 "testing" 16 17 "github.com/jingcheng-WU/gonum/cmplxs" 18 "github.com/jingcheng-WU/gonum/floats" 19 "github.com/jingcheng-WU/gonum/floats/scalar" 20 ) 21 22 func TestRfft(t *testing.T) { 23 t.Parallel() 24 // Golden value cases where golden values have 25 // been obtained from the original Fortran code. 26 for _, test := range rfftTests { 27 const tol = 1e-12 28 testRfft(t, test.n, tol, test.wantiwork, test.wantiifac) 29 } 30 31 // General cases based purely on FFT behaviour. 32 for n := 1; n < 500; n++ { 33 const tol = 1e-9 34 testRfft(t, n, tol, nil, nil) 35 } 36 } 37 38 func testRfft(t *testing.T, n int, tol float64, wantiwork []float64, wantiifac []int) { 39 // Compute the work and factor slices and compare to known values. 40 work := make([]float64, 2*n) 41 ifac := make([]int, 15) 42 Rffti(n, work, ifac) 43 var failed bool 44 if wantiwork != nil && !floats.EqualApprox(work, wantiwork, 1e-6) { 45 failed = true 46 t.Errorf("unexpected work after call to rffti for n=%d", n) 47 } 48 if wantiifac != nil && !reflect.DeepEqual(ifac, wantiifac) { 49 failed = true 50 t.Errorf("unexpected ifac after call to rffti for n=%d", n) 51 } 52 if failed { 53 return 54 } 55 56 // Construct a sequence with known a frequency spectrum and compare 57 // the computed spectrum. 58 modn := n % 2 59 fn := float64(n) 60 nm1 := n - 1 61 x, y, xh := series(n) 62 63 dt := 2 * math.Pi / fn 64 ns2 := (n + 1) / 2 65 if ns2 >= 2 { 66 for k := 1; k < ns2; k++ { //eek 67 var sum1, sum2 float64 68 arg := float64(k) * dt 69 for i := 0; i < n; i++ { 70 arg1 := float64(i) * arg 71 sum1 += x[i] * math.Cos(arg1) 72 sum2 += x[i] * math.Sin(arg1) 73 } 74 y[2*k-1] = sum1 75 y[2*k] = -sum2 76 } 77 } 78 var sum1, sum2 float64 79 for i := 0; i < nm1; i += 2 { 80 sum1 += x[i] 81 sum2 += x[i+1] 82 } 83 if modn == 1 { 84 sum1 += x[n-1] 85 } 86 y[0] = sum1 + sum2 87 if modn == 0 { 88 y[n-1] = sum1 - sum2 89 } 90 91 want := naiveDFTreal(x) 92 93 Rfftf(n, x, work, ifac) 94 95 // Compare the result to a naive DFT implementation. 96 if got := packCoeffs(x); !cmplxs.EqualApprox(got, want, tol) { 97 t.Errorf("unexpected coefficients for n=%d: got:%f want:%f", n, got, want) 98 } 99 100 var rftf float64 101 for i := 0; i < n; i++ { 102 rftf = math.Max(rftf, math.Abs(x[i]-y[i])) 103 x[i] = xh[i] 104 } 105 106 rftf /= fn 107 if !scalar.EqualWithinAbsOrRel(rftf, 0, tol, tol) { 108 t.Errorf("unexpected rftf value for n=%d: got:%f want:0", n, rftf) 109 } 110 111 // Construct a frequency spectrum and compare the computed sequence. 112 for i := 0; i < n; i++ { 113 sum := x[0] / 2 114 arg := float64(i) * dt 115 if ns2 >= 2 { 116 for k := 1; k < ns2; k++ { //eek 117 arg1 := float64(k) * arg 118 sum += x[2*k-1]*math.Cos(arg1) - x[2*k]*math.Sin(arg1) 119 } 120 } 121 if modn == 0 { 122 // This is how it was written in FFTPACK. 123 sum += 0.5 * math.Pow(-1, float64(i)) * x[n-1] 124 } 125 y[i] = 2 * sum 126 } 127 Rfftb(n, x, work, ifac) 128 var rftb float64 129 for i := 0; i < n; i++ { 130 rftb = math.Max(rftb, math.Abs(x[i]-y[i])) 131 x[i] = xh[i] 132 y[i] = xh[i] 133 } 134 if !scalar.EqualWithinAbsOrRel(rftb, 0, tol, tol) { 135 t.Errorf("unexpected rftb value for n=%d: got:%g want:0", n, rftb) 136 } 137 138 // Check that Rfftb and Rfftf are inverses. 139 Rfftb(n, y, work, ifac) 140 Rfftf(n, y, work, ifac) 141 cf := 1.0 / fn 142 var rftfb float64 143 for i := 0; i < n; i++ { 144 rftfb = math.Max(rftfb, math.Abs(cf*y[i]-x[i])) 145 } 146 if !scalar.EqualWithinAbsOrRel(rftfb, 0, tol, tol) { 147 t.Errorf("unexpected rftfb value for n=%d: got:%f want:0", n, rftfb) 148 } 149 } 150 151 // naiveDFTreal is the naive O(n^2) DFT implementation. 152 func naiveDFTreal(x []float64) (y []complex128) { 153 y = make([]complex128, len(x)) 154 dt := -2 * math.Pi / float64(len(x)) 155 for i := range x { 156 arg1 := float64(i) * dt 157 for k, xv := range x { 158 arg2 := float64(k) * arg1 159 y[i] += complex(xv*math.Cos(arg2), xv*math.Sin(arg2)) 160 } 161 } 162 return y[:len(x)/2+1] 163 } 164 165 func packCoeffs(x []float64) []complex128 { 166 y := make([]complex128, len(x)/2+1) 167 y[0] = complex(x[0], 0) 168 if len(x) < 2 { 169 return y 170 } 171 if len(x)%2 == 1 { 172 y[len(y)-1] = complex(x[len(x)-2], x[len(x)-1]) 173 } else { 174 y[len(y)-1] = complex(x[len(x)-1], 0) 175 } 176 for i := 1; i < len(y)-1; i++ { 177 y[i] = complex(x[2*i-1], x[2*i]) 178 } 179 return y 180 } 181 182 var rfftTests = []struct { 183 n int 184 185 // The following two fields are added as there is no unit testing in 186 // FFTPACK for RFFTI. 187 // 188 // wantiwork is obtained from the FFTPACK test.f with modification. 189 // The W array is zeroed at each iteration and the first 2n elements 190 // of W are printed after the call to RFFTI. 191 wantiwork []float64 192 // wantiifac is obtained from the FFTPACK rffti1.f with modification. 193 // The IFAC array is zeroed at each iteration of test.f and the 15 elements 194 // of IFAC are printed before RFFTI1 returns. 195 wantiifac []int 196 }{ 197 { 198 n: 120, 199 200 wantiwork: []float64{ 201 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 202 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 203 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 204 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 205 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 206 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 207 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 208 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 209 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 210 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 211 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 212 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 213 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 214 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 215 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 216 0.9986295, 0.5233596e-01, 0.9945219, 0.1045285, 0.9876884, 0.1564345, 0.9781476, 0.2079117, 217 0.9659258, 0.2588190, 0.9510565, 0.3090170, 0.9335804, 0.3583679, 0.9135454, 0.4067366, 218 0.8910065, 0.4539905, 0.8660254, 0.5000000, 0.8386706, 0.5446391, 0.8090170, 0.5877852, 219 0.7771459, 0.6293204, 0.7431448, 0.6691306, 0.7071068, 0.7071068, 0.6691306, 0.7431449, 220 0.6293204, 0.7771460, 0.5877852, 0.8090170, 0.5446390, 0.8386706, 0.5000000, 0.8660254, 221 0.4539905, 0.8910065, 0.4067366, 0.9135455, 0.3583679, 0.9335805, 0.3090170, 0.9510565, 222 0.2588191, 0.9659258, 0.2079117, 0.9781476, 0.1564344, 0.9876884, 0.1045284, 0.9945219, 223 0.5233597e-01, 0.9986295, 0.000000, 0.000000, 0.9945219, 0.1045285, 0.9781476, 0.2079117, 224 0.9510565, 0.3090170, 0.9135454, 0.4067366, 0.8660254, 0.5000000, 0.8090170, 0.5877852, 225 0.7431448, 0.6691306, 0.000000, 0.9781476, 0.2079117, 0.9135454, 0.4067366, 0.8090170, 226 0.5877852, 0.6691306, 0.7431449, 0.5000000, 0.8660254, 0.3090170, 0.9510565, 0.1045284, 227 0.9945219, 0.000000, 0.9510565, 0.3090170, 0.8090170, 0.5877852, 0.5877852, 0.8090170, 228 0.3090170, 0.9510565, -0.4371139e-07, 1.000000, -0.3090170, 0.9510565, -0.5877852, 0.8090170, 229 0.000000, 0.9135454, 0.4067366, 0.6691306, 0.7431449, 0.000000, 0.6691306, 0.7431449, 230 -0.1045285, 0.9945219, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 231 }, 232 wantiifac: []int{120, 4, 2, 4, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 233 }, 234 { 235 n: 54, 236 237 wantiwork: []float64{ 238 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 239 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 240 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 241 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 242 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 243 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 244 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.9932383, 0.1160929, 245 0.9730449, 0.2306159, 0.9396926, 0.3420201, 0.8936327, 0.4487992, 0.8354878, 0.5495090, 246 0.7660444, 0.6427876, 0.6862416, 0.7273737, 0.5971586, 0.8021232, 0.5000000, 0.8660254, 247 0.3960797, 0.9182161, 0.2868032, 0.9579895, 0.1736482, 0.9848077, 0.5814485e-01, 0.9983082, 248 0.000000, 0.9730449, 0.2306159, 0.8936327, 0.4487992, 0.7660444, 0.6427876, 0.5971586, 249 0.8021232, 0.000000, 0.8936327, 0.4487992, 0.5971586, 0.8021232, 0.1736482, 0.9848077, 250 -0.2868032, 0.9579895, 0.000000, 0.7660444, 0.6427876, 0.000000, 0.1736482, 0.9848077, 251 0.000000, 0.000000, 0.000000, 0.000000, 252 }, 253 wantiifac: []int{54, 4, 2, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 254 }, 255 { 256 n: 49, 257 258 wantiwork: []float64{ 259 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 260 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 261 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 262 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 263 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 264 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 265 0.000000, 0.9917900, 0.1278772, 0.9672949, 0.2536546, 0.9269168, 0.3752670, 0.000000, 266 0.9672949, 0.2536546, 0.8713187, 0.4907176, 0.7183493, 0.6956826, 0.000000, 0.9269168, 267 0.3752670, 0.7183493, 0.6956826, 0.4047833, 0.9144127, 0.000000, 0.8713187, 0.4907176, 268 0.5183925, 0.8551428, 0.3205151e-01, 0.9994862, 0.000000, 0.8014136, 0.5981106, 0.2845275, 269 0.9586679, -0.3453652, 0.9384683, 0.000000, 0.7183493, 0.6956826, 0.3205151e-01, 0.9994862, 270 -0.6723010, 0.7402779, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 271 0.000000, 0.000000, 272 }, 273 wantiifac: []int{49, 2, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 274 }, 275 { 276 n: 32, 277 278 wantiwork: []float64{ 279 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 280 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 281 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 282 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 283 0.9807853, 0.1950903, 0.9238795, 0.3826835, 0.8314696, 0.5555702, 0.7071068, 0.7071068, 284 0.5555702, 0.8314697, 0.3826834, 0.9238795, 0.1950902, 0.9807853, 0.000000, 0.000000, 285 0.9238795, 0.3826835, 0.000000, 0.000000, 0.7071068, 0.7071068, 0.000000, 0.000000, 286 0.3826834, 0.9238795, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 287 }, 288 wantiifac: []int{32, 3, 2, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 289 }, 290 { 291 n: 25, 292 293 wantiwork: []float64{ 294 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 295 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 296 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 297 0.00000, 0.968583, 0.248690, 0.876307, 0.481754, 0.00000, 0.876307, 0.481754, 298 0.535827, 0.844328, 0.00000, 0.728969, 0.684547, 0.627904e-01, 0.998027, 0.00000, 299 0.535827, 0.844328, -0.425779, 0.904827, 0.00000, 0.00000, 0.00000, 0.00000, 300 0.00000, 0.00000, 301 }, 302 wantiifac: []int{25, 2, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 303 }, 304 { 305 n: 4, 306 307 wantiwork: []float64{ 308 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 309 }, 310 wantiifac: []int{4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 311 }, 312 { 313 n: 3, 314 315 wantiwork: []float64{ 316 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 317 }, 318 wantiifac: []int{3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 319 }, 320 { 321 n: 2, 322 323 wantiwork: []float64{ 324 0.000000, 0.000000, 0.000000, 0.000000, 325 }, 326 wantiifac: []int{2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 327 }, 328 } 329 330 func TestCfft(t *testing.T) { 331 t.Parallel() 332 333 const tol = 1e-12 334 335 // Golden value cases where golden values have 336 // been obtained from the original Fortran code. 337 for _, test := range cfftTests { 338 testCfft(t, test.n, tol, test.wantiwork, test.wantiifac) 339 } 340 341 // General cases based purely on FFT behaviour. 342 for n := 1; n < 500; n++ { 343 testCfft(t, n, tol, nil, nil) 344 } 345 } 346 347 func testCfft(t *testing.T, n int, tol float64, wantiwork []float64, wantiifac []int) { 348 // Compute the work and factor slices and compare to known values. 349 work := make([]float64, 4*n) 350 ifac := make([]int, 15) 351 Cffti(n, work, ifac) 352 var failed bool 353 if wantiwork != nil && !floats.EqualApprox(work, wantiwork, 1e-6) { 354 failed = true 355 t.Errorf("unexpected work after call to cffti for n=%d", n) 356 } 357 if wantiifac != nil && !reflect.DeepEqual(ifac, wantiifac) { 358 failed = true 359 t.Errorf("unexpected ifac after call to cffti for n=%d", n) 360 } 361 if failed { 362 return 363 } 364 365 // Construct a sequence with known a frequency spectrum and compare 366 // the computed spectrum. 367 fn := float64(n) 368 cn := complex(fn, 0) 369 370 x, y1 := cmplxSeries(n) 371 372 want := naiveDFT(x) 373 374 cx := cmplxAsFloat(x) 375 Cfftf(n, cx, work, ifac) 376 x = floatAsCmplx(cx) 377 378 // Compare the result to a naive DFT implementation. 379 if !cmplxs.EqualApprox(x, want, tol*10) { 380 t.Errorf("unexpected coefficients for n=%d: got:%f want:%f", n, x, want) 381 } 382 383 var cftf float64 384 for i := 0; i < n; i++ { 385 cftf = math.Max(cftf, cmplx.Abs(x[i]-y1[i])) 386 x[i] /= cn 387 } 388 cftf /= fn 389 390 if !scalar.EqualWithinAbsOrRel(cftf, 0, tol, tol) { 391 t.Errorf("unexpected cftf value for n=%d: got:%f want:0", n, cftf) 392 } 393 394 // Construct a frequency spectrum and compare the computed sequence. 395 y2 := updatedCmplxSeries(x) 396 397 cx = cmplxAsFloat(x) 398 Cfftb(n, cx, work, ifac) 399 x = floatAsCmplx(cx) 400 401 var cftb float64 402 for i := 0; i < n; i++ { 403 cftb = math.Max(cftb, cmplx.Abs(x[i]-y2[i])) 404 x[i] = y2[i] 405 } 406 407 if !scalar.EqualWithinAbsOrRel(cftb, 0, tol, tol) { 408 t.Errorf("unexpected cftb value for n=%d: got:%f want:0", n, cftb) 409 } 410 411 // Check that Cfftb and Cfftf are inverses. 412 cx = cmplxAsFloat(x) 413 Cfftf(n, cx, work, ifac) 414 Cfftb(n, cx, work, ifac) 415 x = floatAsCmplx(cx) 416 417 var cftfb float64 418 for i := 0; i < n; i++ { 419 cftfb = math.Max(cftfb, cmplx.Abs(x[i]/cn-y2[i])) 420 } 421 422 if !scalar.EqualWithinAbsOrRel(cftfb, 0, tol, tol) { 423 t.Errorf("unexpected cftfb value for n=%d: got:%f want:0", n, cftfb) 424 } 425 } 426 427 // naiveDFT is the naive O(n^2) DFT implementation. 428 func naiveDFT(x []complex128) (y []complex128) { 429 y = make([]complex128, len(x)) 430 dt := -2 * math.Pi / float64(len(x)) 431 for i := range x { 432 arg1 := float64(i) * dt 433 for k, xv := range x { 434 arg2 := float64(k) * arg1 435 y[i] += complex(math.Cos(arg2), math.Sin(arg2)) * xv 436 } 437 } 438 return y 439 } 440 func cmplxSeries(n int) (x, y []complex128) { 441 x = make([]complex128, n) 442 for i := 0; i < n; i++ { 443 x[i] = complex(math.Cos(math.Sqrt2*float64(i+1)), math.Sin(math.Sqrt2*float64((i+1)*(i+1)))) 444 } 445 446 y = make([]complex128, n) 447 dt := 2 * math.Pi / float64(n) 448 for i := 0; i < n; i++ { 449 arg1 := -float64(i) * dt 450 for k := 0; k < n; k++ { 451 arg2 := float64(k) * arg1 452 y[i] += complex(math.Cos(arg2), math.Sin(arg2)) * x[k] 453 } 454 } 455 return x, y 456 } 457 458 func updatedCmplxSeries(x []complex128) (y []complex128) { 459 y = make([]complex128, len(x)) 460 dt := 2 * math.Pi / float64(len(x)) 461 for i := range x { 462 arg1 := float64(i) * dt 463 for k, xv := range x { 464 arg2 := float64(k) * arg1 465 y[i] += complex(math.Cos(arg2), math.Sin(arg2)) * xv 466 } 467 } 468 return y 469 } 470 471 func cmplxAsFloat(c []complex128) []float64 { 472 f := make([]float64, 2*len(c)) 473 for i, v := range c { 474 f[2*i] = real(v) 475 f[2*i+1] = imag(v) 476 } 477 return f 478 } 479 480 func floatAsCmplx(f []float64) []complex128 { 481 c := make([]complex128, len(f)/2) 482 for i := range c { 483 c[i] = complex(f[2*i], f[2*i+1]) 484 } 485 return c 486 } 487 488 var cfftTests = []struct { 489 n int 490 491 // The following two fields are added as there is no unit testing in 492 // FFTPACK for RFFTI. 493 // 494 // wantiwork is obtained from the FFTPACK test.f with modification. 495 // The W array is zeroed at each iteration and the first 4n elements 496 // of W are printed after the call to RFFTI. 497 wantiwork []float64 498 // wantiifac is obtained from the FFTPACK rffti1.f with modification. 499 // The IFAC array is zeroed at each iteration of test.f and the 15 elements 500 // of IFAC are printed before RFFTI1 returns. 501 wantiifac []int 502 }{ 503 { 504 n: 120, 505 506 wantiwork: []float64{ 507 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 508 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 509 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 510 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 511 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 512 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 513 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 514 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 515 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 516 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 517 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 518 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 519 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 520 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 521 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 522 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 523 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 524 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 525 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 526 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 527 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 528 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 529 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 530 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 531 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 532 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 533 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 534 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 535 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 536 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 537 1.00000, 0.00000, 0.998630, 0.523360e-01, 0.994522, 0.104528, 0.987688, 0.156434, 538 0.978148, 0.207912, 0.965926, 0.258819, 0.951057, 0.309017, 0.933580, 0.358368, 539 0.913545, 0.406737, 0.891007, 0.453991, 0.866025, 0.500000, 0.838671, 0.544639, 540 0.809017, 0.587785, 0.777146, 0.629320, 0.743145, 0.669131, 0.707107, 0.707107, 541 0.669131, 0.743145, 0.629320, 0.777146, 0.587785, 0.809017, 0.544639, 0.838671, 542 0.500000, 0.866025, 0.453991, 0.891007, 0.406737, 0.913545, 0.358368, 0.933580, 543 0.309017, 0.951057, 0.258819, 0.965926, 0.207912, 0.978148, 0.156434, 0.987688, 544 0.104528, 0.994522, 0.523360e-01, 0.998630, -0.437114e-07, 1.00000, -0.523361e-01, 0.998630, 545 -0.104529, 0.994522, -0.156434, 0.987688, -0.207912, 0.978148, -0.258819, 0.965926, 546 -0.309017, 0.951056, -0.358368, 0.933580, -0.406737, 0.913545, -0.453991, 0.891006, 547 -0.500000, 0.866025, -0.544639, 0.838671, -0.587785, 0.809017, -0.629321, 0.777146, 548 -0.669131, 0.743145, -0.707107, 0.707107, -0.743145, 0.669130, -0.777146, 0.629320, 549 -0.809017, 0.587785, -0.838671, 0.544639, -0.866025, 0.500000, -0.891007, 0.453990, 550 -0.913545, 0.406737, -0.933580, 0.358368, -0.951057, 0.309017, -0.965926, 0.258819, 551 -0.978148, 0.207912, -0.987688, 0.156434, -0.994522, 0.104528, -0.998630, 0.523358e-01, 552 1.00000, 0.00000, 0.994522, 0.104528, 0.978148, 0.207912, 0.951057, 0.309017, 553 0.913545, 0.406737, 0.866025, 0.500000, 0.809017, 0.587785, 0.743145, 0.669131, 554 0.669131, 0.743145, 0.587785, 0.809017, 0.500000, 0.866025, 0.406737, 0.913545, 555 0.309017, 0.951057, 0.207912, 0.978148, 0.104528, 0.994522, -0.437114e-07, 1.00000, 556 -0.104529, 0.994522, -0.207912, 0.978148, -0.309017, 0.951056, -0.406737, 0.913545, 557 1.00000, 0.00000, 0.978148, 0.207912, 0.913545, 0.406737, 0.809017, 0.587785, 558 0.669131, 0.743145, 0.500000, 0.866025, 0.309017, 0.951057, 0.104528, 0.994522, 559 -0.104529, 0.994522, -0.309017, 0.951056, -0.500000, 0.866025, -0.669131, 0.743145, 560 -0.809017, 0.587785, -0.913545, 0.406737, -0.978148, 0.207912, -1.00000, -0.874228e-07, 561 -0.978148, -0.207912, -0.913545, -0.406737, -0.809017, -0.587785, -0.669131, -0.743145, 562 1.00000, 0.00000, 0.951057, 0.309017, 0.809017, 0.587785, 0.587785, 0.809017, 563 0.309017, 0.951057, 1.00000, 0.00000, 0.809017, 0.587785, 0.309017, 0.951057, 564 -0.309017, 0.951056, -0.809017, 0.587785, 1.00000, 0.00000, 0.587785, 0.809017, 565 -0.309017, 0.951056, -0.951057, 0.309017, -0.809017, -0.587785, 1.00000, 0.00000, 566 1.00000, 0.00000, 1.00000, 0.00000, 1.00000, 0.00000, 0.309017, -0.951056, 567 }, 568 wantiifac: []int{120, 4, 2, 3, 4, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 569 }, 570 { 571 n: 54, 572 573 wantiwork: []float64{ 574 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 575 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 576 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 577 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 578 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 579 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 580 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 581 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 582 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 583 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 584 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 585 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 586 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 587 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.993238, 0.116093, 588 0.973045, 0.230616, 0.939693, 0.342020, 0.893633, 0.448799, 0.835488, 0.549509, 589 0.766044, 0.642788, 0.686242, 0.727374, 0.597159, 0.802123, 0.500000, 0.866025, 590 0.396080, 0.918216, 0.286803, 0.957990, 0.173648, 0.984808, 0.581449e-01, 0.998308, 591 -0.581448e-01, 0.998308, -0.173648, 0.984808, -0.286803, 0.957990, -0.396080, 0.918216, 592 -0.500000, 0.866025, -0.597159, 0.802123, -0.686242, 0.727374, -0.766044, 0.642788, 593 -0.835488, 0.549509, -0.893633, 0.448799, -0.939693, 0.342020, -0.973045, 0.230616, 594 -0.993238, 0.116093, 1.00000, 0.00000, 0.973045, 0.230616, 0.893633, 0.448799, 595 0.766044, 0.642788, 0.597159, 0.802123, 0.396080, 0.918216, 0.173648, 0.984808, 596 -0.581448e-01, 0.998308, -0.286803, 0.957990, 1.00000, 0.00000, 0.893633, 0.448799, 597 0.597159, 0.802123, 0.173648, 0.984808, -0.286803, 0.957990, -0.686242, 0.727374, 598 -0.939693, 0.342020, -0.993238, -0.116093, -0.835488, -0.549509, 1.00000, 0.00000, 599 0.766044, 0.642788, 0.173648, 0.984808, 1.00000, 0.00000, 0.173648, 0.984808, 600 -0.939693, 0.342020, 1.00000, 0.00000, 1.00000, 0.00000, -0.500000, -0.866025, 601 }, 602 wantiifac: []int{54, 4, 2, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 603 }, 604 { 605 n: 49, 606 607 wantiwork: []float64{ 608 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 609 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 610 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 611 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 612 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 613 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 614 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 615 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 616 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 617 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 618 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 619 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 620 0.00000, 0.00000, 0.623490, 0.781832, 0.991790, 0.127877, 0.967295, 0.253655, 621 0.926917, 0.375267, 0.871319, 0.490718, 0.801414, 0.598111, 0.718349, 0.695683, 622 -0.222521, 0.974928, 0.967295, 0.253655, 0.871319, 0.490718, 0.718349, 0.695683, 623 0.518393, 0.855143, 0.284527, 0.958668, 0.320515e-01, 0.999486, -0.900969, 0.433884, 624 0.926917, 0.375267, 0.718349, 0.695683, 0.404783, 0.914413, 0.320515e-01, 0.999486, 625 -0.345365, 0.938468, -0.672301, 0.740278, -0.900969, -0.433884, 0.871319, 0.490718, 626 0.518393, 0.855143, 0.320515e-01, 0.999486, -0.462538, 0.886599, -0.838088, 0.545535, 627 -0.997945, 0.640701e-01, -0.222521, -0.974928, 0.801414, 0.598111, 0.284527, 0.958668, 628 -0.345365, 0.938468, -0.838088, 0.545535, -0.997945, -0.640705e-01, -0.761446, -0.648229, 629 0.623490, -0.781831, 0.718349, 0.695683, 0.320515e-01, 0.999486, -0.672301, 0.740278, 630 -0.997945, 0.640701e-01, -0.761446, -0.648228, -0.960227e-01, -0.995379, 0.623490, 0.781832, 631 -0.222521, 0.974928, -0.900969, 0.433884, -0.900969, -0.433884, -0.222521, -0.974928, 632 0.623490, -0.781831, 0.623490, -0.781831, 633 }, 634 wantiifac: []int{49, 2, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 635 }, 636 { 637 n: 32, 638 639 wantiwork: []float64{ 640 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 641 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 642 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 643 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 644 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 645 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 646 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 647 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 648 1.00000, 0.00000, 0.980785, 0.195090, 0.923880, 0.382683, 0.831470, 0.555570, 649 0.707107, 0.707107, 0.555570, 0.831470, 0.382683, 0.923880, 0.195090, 0.980785, 650 -0.437114e-07, 1.00000, -0.195090, 0.980785, -0.382684, 0.923880, -0.555570, 0.831470, 651 -0.707107, 0.707107, -0.831470, 0.555570, -0.923880, 0.382683, -0.980785, 0.195090, 652 1.00000, 0.00000, 0.923880, 0.382683, 0.707107, 0.707107, 0.382683, 0.923880, 653 1.00000, 0.00000, 0.707107, 0.707107, -0.437114e-07, 1.00000, -0.707107, 0.707107, 654 1.00000, 0.00000, 0.382683, 0.923880, -0.707107, 0.707107, -0.923880, -0.382683, 655 1.00000, 0.00000, 1.00000, 0.00000, 1.00000, 0.00000, 0.119249e-07, -1.00000, 656 }, 657 wantiifac: []int{32, 3, 2, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 658 }, 659 { 660 n: 25, 661 662 wantiwork: []float64{ 663 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 664 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 665 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 666 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 667 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 668 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 669 0.00000, 0.00000, 1.00000, 0.00000, 0.968583, 0.248690, 0.876307, 0.481754, 670 0.728969, 0.684547, 0.535827, 0.844328, 1.00000, 0.00000, 0.876307, 0.481754, 671 0.535827, 0.844328, 0.627904e-01, 0.998027, -0.425779, 0.904827, 1.00000, 0.00000, 672 0.728969, 0.684547, 0.627904e-01, 0.998027, -0.637424, 0.770513, -0.992115, 0.125333, 673 1.00000, 0.00000, 0.535827, 0.844328, -0.425779, 0.904827, -0.992115, 0.125333, 674 -0.637424, -0.770513, 1.00000, 0.00000, 1.00000, 0.00000, 1.00000, 0.00000, 675 1.00000, 0.00000, 0.309017, -0.951056, 676 }, 677 wantiifac: []int{25, 2, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 678 }, 679 { 680 n: 4, 681 682 wantiwork: []float64{ 683 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 684 1.00000, 0.00000, 1.00000, 0.00000, 1.00000, 0.00000, 0.119249e-07, -1.00000, 685 }, 686 wantiifac: []int{4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 687 }, 688 { 689 n: 3, 690 691 wantiwork: []float64{ 692 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 693 1.00000, 0.00000, -0.500000, -0.866025, 694 }, 695 wantiifac: []int{3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 696 }, 697 { 698 n: 2, 699 700 wantiwork: []float64{ 701 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, -1.00000, -0.874228e-07, 702 }, 703 wantiifac: []int{2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 704 }, 705 } 706 707 func TestSint(t *testing.T) { 708 t.Parallel() 709 const tol = 1e-12 710 for _, test := range sintTests { 711 // Compute the work and factor slices and compare to known values. 712 work := make([]float64, 5*(test.n)/2) 713 ifac := make([]int, 15) 714 Sinti(test.n-1, work, ifac) 715 var failed bool 716 if !floats.EqualApprox(work, test.wantiwork, 1e-6) { 717 failed = true 718 t.Errorf("unexpected work after call to sinti for n=%d", test.n) 719 } 720 if !reflect.DeepEqual(ifac, test.wantiifac) { 721 failed = true 722 t.Errorf("unexpected ifac after call to sinti for n=%d", test.n) 723 } 724 if failed { 725 continue 726 } 727 728 // Construct a sequence with known a frequency spectrum and compare 729 // the computed spectrum. 730 fn := float64(test.n) 731 x, _, xh := series(test.n - 1) 732 y := make([]float64, test.n-1) 733 734 dt := math.Pi / fn 735 for i := 0; i < test.n-1; i++ { 736 arg1 := float64(i+1) * dt 737 for k := 0; k < test.n-1; k++ { //eek 738 y[i] += x[k] * math.Sin(float64(k+1)*arg1) 739 } 740 y[i] *= 2 741 } 742 743 Sint(test.n-1, x, work, ifac) 744 745 cf := 0.5 / fn 746 var sintt float64 747 for i := 0; i < test.n-1; i++ { 748 sintt = math.Max(sintt, math.Abs(x[i]-y[i])) 749 x[i] = xh[i] 750 y[i] = x[i] 751 } 752 sintt *= cf 753 if !scalar.EqualWithinAbsOrRel(sintt, 0, tol, tol) { 754 t.Errorf("unexpected sintt value for n=%d: got:%f want:0", test.n, sintt) 755 } 756 757 // Check that the transform is its own inverse. 758 Sint(test.n-1, x, work, ifac) 759 Sint(test.n-1, x, work, ifac) 760 761 var sintfb float64 762 for i := 0; i < test.n-1; i++ { 763 sintfb = math.Max(sintfb, math.Abs(cf*x[i]-y[i])) 764 } 765 if !scalar.EqualWithinAbsOrRel(sintfb, 0, tol, tol) { 766 t.Errorf("unexpected sintfb value for n=%d: got:%f want:0", test.n, sintfb) 767 } 768 } 769 } 770 771 var sintTests = []struct { 772 n int 773 774 // The following two fields are added as there is no unit testing in 775 // FFTPACK for SINTI. 776 // 777 // wantiwork is obtained from the FFTPACK test.f with modification. 778 // The W array is zeroed at each iteration and the first 2.5n elements 779 // of W are printed after the call to RFFTI. 780 wantiwork []float64 781 // wantiifac is obtained from the FFTPACK sint.f with modification. 782 // The IFAC array is zeroed at each iteration of test.f and the 15 elements 783 // of IFAC are printed before RFFTI returns. 784 wantiifac []int 785 }{ 786 { 787 n: 120, 788 789 wantiwork: []float64{ 790 0.5235390e-01, 0.1046719, 0.1569182, 0.2090569, 0.2610524, 0.3128690, 0.3644710, 0.4158234, 791 0.4668908, 0.5176381, 0.5680307, 0.6180340, 0.6676137, 0.7167359, 0.7653669, 0.8134733, 792 0.8610222, 0.9079810, 0.9543176, 1.000000, 1.044997, 1.089278, 1.132812, 1.175570, 793 1.217523, 1.258641, 1.298896, 1.338261, 1.376709, 1.414214, 1.450749, 1.486290, 794 1.520812, 1.554292, 1.586707, 1.618034, 1.648252, 1.677341, 1.705280, 1.732051, 795 1.757634, 1.782013, 1.805171, 1.827091, 1.847759, 1.867161, 1.885283, 1.902113, 796 1.917639, 1.931852, 1.944740, 1.956295, 1.966510, 1.975377, 1.982890, 1.989044, 797 1.993835, 1.997259, 1.999315, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 798 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 799 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 800 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 801 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 802 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 803 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 804 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 805 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 806 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 807 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 808 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 809 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 810 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 811 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 812 0.000000, 0.000000, 0.000000, 0.9986295, 0.5233596e-01, 0.9945219, 0.1045285, 0.9876884, 813 0.1564345, 0.9781476, 0.2079117, 0.9659258, 0.2588190, 0.9510565, 0.3090170, 0.9335804, 814 0.3583679, 0.9135454, 0.4067366, 0.8910065, 0.4539905, 0.8660254, 0.5000000, 0.8386706, 815 0.5446391, 0.8090170, 0.5877852, 0.7771459, 0.6293204, 0.7431448, 0.6691306, 0.7071068, 816 0.7071068, 0.6691306, 0.7431449, 0.6293204, 0.7771460, 0.5877852, 0.8090170, 0.5446390, 817 0.8386706, 0.5000000, 0.8660254, 0.4539905, 0.8910065, 0.4067366, 0.9135455, 0.3583679, 818 0.9335805, 0.3090170, 0.9510565, 0.2588191, 0.9659258, 0.2079117, 0.9781476, 0.1564344, 819 0.9876884, 0.1045284, 0.9945219, 0.5233597e-01, 0.9986295, 0.000000, 0.000000, 0.9945219, 820 0.1045285, 0.9781476, 0.2079117, 0.9510565, 0.3090170, 0.9135454, 0.4067366, 0.8660254, 821 0.5000000, 0.8090170, 0.5877852, 0.7431448, 0.6691306, 0.000000, 0.9781476, 0.2079117, 822 0.9135454, 0.4067366, 0.8090170, 0.5877852, 0.6691306, 0.7431449, 0.5000000, 0.8660254, 823 0.3090170, 0.9510565, 0.1045284, 0.9945219, 0.000000, 0.9510565, 0.3090170, 0.8090170, 824 0.5877852, 0.5877852, 0.8090170, 0.3090170, 0.9510565, -0.4371139e-07, 1.000000, -0.3090170, 825 0.9510565, -0.5877852, 0.8090170, 0.000000, 0.9135454, 0.4067366, 0.6691306, 0.7431449, 826 0.000000, 0.6691306, 0.7431449, -0.1045285, 0.9945219, 0.000000, 0.000000, 0.000000, 827 0.000000, 0.000000, 0.000000, 0.1681558e-42, 828 }, 829 wantiifac: []int{120, 4, 2, 4, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 830 }, 831 { 832 n: 54, 833 834 wantiwork: []float64{ 835 0.1162897, 0.2321858, 0.3472964, 0.4612317, 0.5736065, 0.6840402, 0.7921596, 0.8975984, 836 1.000000, 1.099018, 1.194317, 1.285575, 1.372483, 1.454747, 1.532089, 1.604246, 837 1.670976, 1.732051, 1.787265, 1.836432, 1.879385, 1.915979, 1.946090, 1.969615, 838 1.986477, 1.996616, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 839 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 840 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 841 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 842 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 843 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 844 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 845 0.9932383, 0.1160929, 0.9730449, 0.2306159, 0.9396926, 0.3420201, 0.8936327, 0.4487992, 846 0.8354878, 0.5495090, 0.7660444, 0.6427876, 0.6862416, 0.7273737, 0.5971586, 0.8021232, 847 0.5000000, 0.8660254, 0.3960797, 0.9182161, 0.2868032, 0.9579895, 0.1736482, 0.9848077, 848 0.5814485e-01, 0.9983082, 0.000000, 0.9730449, 0.2306159, 0.8936327, 0.4487992, 0.7660444, 849 0.6427876, 0.5971586, 0.8021232, 0.000000, 0.8936327, 0.4487992, 0.5971586, 0.8021232, 850 0.1736482, 0.9848077, -0.2868032, 0.9579895, 0.000000, 0.7660444, 0.6427876, 0.000000, 851 0.1736482, 0.9848077, 0.000000, 0.000000, 0.000000, 0.000000, 0.7567012e-43, 852 }, 853 wantiifac: []int{54, 4, 2, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 854 }, 855 { 856 n: 49, 857 858 wantiwork: []float64{ 859 0.1281404, 0.2557543, 0.3823173, 0.5073092, 0.6302165, 0.7505341, 0.8677675, 0.9814351, 860 1.091070, 1.196221, 1.296457, 1.391365, 1.480556, 1.563663, 1.640345, 1.710286, 861 1.773199, 1.828825, 1.876937, 1.917336, 1.949856, 1.974364, 1.990758, 1.998972, 862 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 863 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 864 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 865 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 866 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 867 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 868 0.000000, 0.9917900, 0.1278772, 0.9672949, 0.2536546, 0.9269168, 0.3752670, 0.000000, 869 0.9672949, 0.2536546, 0.8713187, 0.4907176, 0.7183493, 0.6956826, 0.000000, 0.9269168, 870 0.3752670, 0.7183493, 0.6956826, 0.4047833, 0.9144127, 0.000000, 0.8713187, 0.4907176, 871 0.5183925, 0.8551428, 0.3205151e-01, 0.9994862, 0.000000, 0.8014136, 0.5981106, 0.2845275, 872 0.9586679, -0.3453652, 0.9384683, 0.000000, 0.7183493, 0.6956826, 0.3205151e-01, 0.9994862, 873 -0.6723010, 0.7402779, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 874 0.000000, 0.000000, 875 }, 876 wantiifac: []int{49, 2, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 877 }, 878 { 879 n: 32, 880 881 wantiwork: []float64{ 882 0.1960343, 0.3901806, 0.5805693, 0.7653669, 0.9427935, 1.111140, 1.268787, 1.414214, 883 1.546021, 1.662939, 1.763843, 1.847759, 1.913881, 1.961571, 1.990369, 0.000000, 884 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 885 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 886 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 887 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.9807853, 888 0.1950903, 0.9238795, 0.3826835, 0.8314696, 0.5555702, 0.7071068, 0.7071068, 0.5555702, 889 0.8314697, 0.3826834, 0.9238795, 0.1950902, 0.9807853, 0.000000, 0.000000, 0.9238795, 890 0.3826835, 0.000000, 0.000000, 0.7071068, 0.7071068, 0.000000, 0.000000, 0.3826834, 891 0.9238795, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.4484155e-43, 892 }, 893 wantiifac: []int{32, 3, 2, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 894 }, 895 { 896 n: 25, 897 898 wantiwork: []float64{ 899 0.2506665, 0.4973798, 0.7362491, 0.9635074, 1.175570, 1.369094, 1.541027, 1.688656, 900 1.809654, 1.902113, 1.964575, 1.996053, 0.000000, 0.000000, 0.000000, 0.000000, 901 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 902 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 903 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.9685832, 0.2486899, 0.8763067, 904 0.4817537, 0.000000, 0.8763067, 0.4817537, 0.5358267, 0.8443279, 0.000000, 0.7289686, 905 0.6845471, 0.6279038e-01, 0.9980267, 0.000000, 0.5358267, 0.8443279, -0.4257794, 0.9048270, 906 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 907 }, 908 wantiifac: []int{25, 2, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 909 }, 910 { 911 n: 4, 912 913 wantiwork: []float64{ 914 1.414214, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 915 0.000000, 0.000000, 916 }, 917 wantiifac: []int{4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 918 }, 919 { 920 n: 3, 921 922 wantiwork: []float64{ 923 1.732051, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 924 }, 925 wantiifac: []int{3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 926 }, 927 { 928 n: 2, 929 930 wantiwork: []float64{ 931 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 932 }, 933 wantiifac: []int{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 934 }, 935 } 936 937 func TestCost(t *testing.T) { 938 t.Parallel() 939 const tol = 1e-12 940 for _, test := range costTests { 941 // Compute the work and factor slices and compare to known values. 942 work := make([]float64, 3*(test.n+1)) 943 ifac := make([]int, 15) 944 Costi(test.n+1, work, ifac) 945 var failed bool 946 if !floats.EqualApprox(work, test.wantiwork, 1e-6) { 947 failed = true 948 t.Errorf("unexpected work after call to costi for n=%d", test.n) 949 } 950 if !reflect.DeepEqual(ifac, test.wantiifac) { 951 failed = true 952 t.Errorf("unexpected ifac after call to costi for n=%d", test.n) 953 } 954 if failed { 955 continue 956 } 957 958 // Construct a sequence with known a frequency spectrum and compare 959 // the computed spectrum. 960 fn := float64(test.n) 961 x, _, xh := series(test.n + 1) 962 y := make([]float64, test.n+1) 963 964 dt := math.Pi / fn 965 for i := 0; i < test.n+1; i++ { 966 y[i] = 0.5 * (x[0] + math.Pow(-1, float64(i))*x[test.n]) 967 arg1 := float64(i) * dt 968 for k := 1; k < test.n; k++ { //eek 969 y[i] += x[k] * math.Cos(float64(k)*arg1) 970 } 971 y[i] *= 2 972 } 973 974 Cost(test.n+1, x, work, ifac) 975 976 cf := 0.5 / fn 977 var costt float64 978 for i := 0; i < test.n; i++ { 979 costt = math.Max(costt, math.Abs(x[i]-y[i])) 980 x[i] = xh[i] 981 y[i] = x[i] 982 } 983 costt *= cf 984 if !scalar.EqualWithinAbsOrRel(costt, 0, tol, tol) { 985 t.Errorf("unexpected costt value for n=%d: got:%f want:0", test.n, costt) 986 } 987 988 // Check that the transform is its own inverse. 989 Cost(test.n+1, x, work, ifac) 990 Cost(test.n+1, x, work, ifac) 991 992 var costfb float64 993 for i := 0; i < test.n-1; i++ { 994 costfb = math.Max(costfb, math.Abs(cf*x[i]-y[i])) 995 } 996 if !scalar.EqualWithinAbsOrRel(costfb, 0, tol, tol) { 997 t.Errorf("unexpected costfb value for n=%d: got:%f want:0", test.n, costfb) 998 } 999 } 1000 } 1001 1002 var costTests = []struct { 1003 n int 1004 1005 // The following two fields are added as there is no unit testing in 1006 // FFTPACK for SINTI. 1007 // 1008 // wantiwork is obtained from the FFTPACK test.f with modification. 1009 // The W array is zeroed at each iteration and the first 3n elements 1010 // of W are printed after the call to RFFTI. 1011 wantiwork []float64 1012 // wantiifac is obtained from the FFTPACK sint.f with modification. 1013 // The IFAC array is zeroed at each iteration of test.f and the 15 elements 1014 // of IFAC are printed before RFFTI returns. 1015 wantiifac []int 1016 }{ 1017 { 1018 n: 120, 1019 1020 wantiwork: []float64{ 1021 0.000000, 0.5235390e-01, 0.1046719, 0.1569182, 0.2090569, 0.2610524, 0.3128690, 0.3644710, 1022 0.4158234, 0.4668908, 0.5176381, 0.5680307, 0.6180340, 0.6676137, 0.7167359, 0.7653669, 1023 0.8134733, 0.8610222, 0.9079810, 0.9543176, 1.000000, 1.044997, 1.089278, 1.132812, 1024 1.175570, 1.217523, 1.258641, 1.298896, 1.338261, 1.376709, 1.414214, 1.450749, 1025 1.486290, 1.520812, 1.554292, 1.586707, 1.618034, 1.648252, 1.677341, 1.705280, 1026 1.732051, 1.757634, 1.782013, 1.805171, 1.827091, 1.847759, 1.867161, 1.885283, 1027 1.902113, 1.917639, 1.931852, 1.944740, 1.956295, 1.966510, 1.975377, 1.982890, 1028 1.989044, 1.993835, 1.997259, 1.999315, 0.000000, 0.5235375e-01, 0.1046719, 0.1569182, 1029 0.2090568, 0.2610523, 0.3128687, 0.3644710, 0.4158233, 0.4668906, 0.5176381, 0.5680307, 1030 0.6180339, 0.6676136, 0.7167357, 0.7653669, 0.8134732, 0.8610221, 0.9079810, 0.9543175, 1031 1.000000, 1.044997, 1.089278, 1.132812, 1.175570, 1.217523, 1.258641, 1.298896, 1032 1.338261, 1.376709, 1.414214, 1.450749, 1.486290, 1.520812, 1.554292, 1.586707, 1033 1.618034, 1.648252, 1.677341, 1.705280, 1.732051, 1.757634, 1.782013, 1.805171, 1034 1.827091, 1.847759, 1.867161, 1.885283, 1.902113, 1.917639, 1.931852, 1.944740, 1035 1.956295, 1.966510, 1.975377, 1.982890, 1.989044, 1.993835, 1.997259, 1.999315, 1036 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1037 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1038 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1039 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1040 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1041 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1042 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1043 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1044 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1045 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1046 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1047 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1048 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1049 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1050 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1051 0.000000, 0.9986295, 0.5233596e-01, 0.9945219, 0.1045285, 0.9876884, 0.1564345, 0.9781476, 1052 0.2079117, 0.9659258, 0.2588190, 0.9510565, 0.3090170, 0.9335804, 0.3583679, 0.9135454, 1053 0.4067366, 0.8910065, 0.4539905, 0.8660254, 0.5000000, 0.8386706, 0.5446391, 0.8090170, 1054 0.5877852, 0.7771459, 0.6293204, 0.7431448, 0.6691306, 0.7071068, 0.7071068, 0.6691306, 1055 0.7431449, 0.6293204, 0.7771460, 0.5877852, 0.8090170, 0.5446390, 0.8386706, 0.5000000, 1056 0.8660254, 0.4539905, 0.8910065, 0.4067366, 0.9135455, 0.3583679, 0.9335805, 0.3090170, 1057 0.9510565, 0.2588191, 0.9659258, 0.2079117, 0.9781476, 0.1564344, 0.9876884, 0.1045284, 1058 0.9945219, 0.5233597e-01, 0.9986295, 0.000000, 0.000000, 0.9945219, 0.1045285, 0.9781476, 1059 0.2079117, 0.9510565, 0.3090170, 0.9135454, 0.4067366, 0.8660254, 0.5000000, 0.8090170, 1060 0.5877852, 0.7431448, 0.6691306, 0.000000, 0.9781476, 0.2079117, 0.9135454, 0.4067366, 1061 0.8090170, 0.5877852, 0.6691306, 0.7431449, 0.5000000, 0.8660254, 0.3090170, 0.9510565, 1062 0.1045284, 0.9945219, 0.000000, 0.9510565, 0.3090170, 0.8090170, 0.5877852, 0.5877852, 1063 0.8090170, 0.3090170, 0.9510565, -0.4371139e-07, 1.000000, -0.3090170, 0.9510565, -0.5877852, 1064 0.8090170, 0.000000, 0.9135454, 0.4067366, 0.6691306, 0.7431449, 0.000000, 0.6691306, 1065 0.7431449, -0.1045285, 0.9945219, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1066 0.000000, 0.1681558e-42, 0.5605194e-44, 1067 }, 1068 wantiifac: []int{120, 4, 2, 4, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1069 }, 1070 { 1071 n: 54, 1072 1073 wantiwork: []float64{ 1074 0.000000, 0.1162897, 0.2321858, 0.3472964, 0.4612317, 0.5736065, 0.6840402, 0.7921596, 1075 0.8975984, 1.000000, 1.099018, 1.194317, 1.285575, 1.372483, 1.454747, 1.532089, 1076 1.604246, 1.670976, 1.732051, 1.787265, 1.836432, 1.879385, 1.915979, 1.946090, 1077 1.969615, 1.986477, 1.996616, 0.000000, 0.1162897, 0.2321858, 0.3472964, 0.4612317, 1078 0.5736064, 0.6840403, 0.7921594, 0.8975984, 1.000000, 1.099018, 1.194317, 1.285575, 1079 1.372483, 1.454747, 1.532089, 1.604246, 1.670976, 1.732051, 1.787265, 1.836432, 1080 1.879385, 1.915979, 1.946090, 1.969615, 1.986477, 1.996616, 0.000000, 0.000000, 1081 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1082 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1083 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1084 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1085 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1086 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1087 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.9932383, 0.1160929, 0.9730449, 1088 0.2306159, 0.9396926, 0.3420201, 0.8936327, 0.4487992, 0.8354878, 0.5495090, 0.7660444, 1089 0.6427876, 0.6862416, 0.7273737, 0.5971586, 0.8021232, 0.5000000, 0.8660254, 0.3960797, 1090 0.9182161, 0.2868032, 0.9579895, 0.1736482, 0.9848077, 0.5814485e-01, 0.9983082, 0.000000, 1091 0.9730449, 0.2306159, 0.8936327, 0.4487992, 0.7660444, 0.6427876, 0.5971586, 0.8021232, 1092 0.000000, 0.8936327, 0.4487992, 0.5971586, 0.8021232, 0.1736482, 0.9848077, -0.2868032, 1093 0.9579895, 0.000000, 0.7660444, 0.6427876, 0.000000, 0.1736482, 0.9848077, 0.000000, 1094 0.000000, 0.000000, 0.000000, 0.7567012e-43, 0.5605194e-44, 1095 }, 1096 wantiifac: []int{54, 4, 2, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1097 }, 1098 { 1099 n: 49, 1100 1101 wantiwork: []float64{ 1102 0.000000, 0.1281404, 0.2557543, 0.3823173, 0.5073092, 0.6302165, 0.7505341, 0.8677675, 1103 0.9814351, 1.091070, 1.196221, 1.296457, 1.391365, 1.480556, 1.563663, 1.640345, 1104 1.710286, 1.773199, 1.828825, 1.876937, 1.917336, 1.949856, 1.974364, 1.990758, 1105 1.998972, 0.6410302e-01, 0.1920458, 0.3191997, 0.4450417, 0.5690550, 0.6907300, 0.8095666, 1106 0.9250766, 1.036785, 1.144233, 1.246980, 1.344602, 1.436699, 1.522892, 1.602827, 1107 1.676176, 1.742637, 1.801938, 1.853834, 1.898111, 1.934590, 1.963118, 1.983580, 1108 1.995891, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1109 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1110 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1111 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1112 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1113 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1114 0.000000, 0.000000, 0.000000, 0.9917900, 0.1278772, 0.9672949, 0.2536546, 0.9269168, 1115 0.3752670, 0.000000, 0.9672949, 0.2536546, 0.8713187, 0.4907176, 0.7183493, 0.6956826, 1116 0.000000, 0.9269168, 0.3752670, 0.7183493, 0.6956826, 0.4047833, 0.9144127, 0.000000, 1117 0.8713187, 0.4907176, 0.5183925, 0.8551428, 0.3205151e-01, 0.9994862, 0.000000, 0.8014136, 1118 0.5981106, 0.2845275, 0.9586679, -0.3453652, 0.9384683, 0.000000, 0.7183493, 0.6956826, 1119 0.3205151e-01, 0.9994862, -0.6723010, 0.7402779, 0.000000, 0.000000, 0.000000, 0.000000, 1120 0.000000, 0.000000, 0.000000, 0.000000, 0.6866362e-43, 0.2802597e-44, 1121 }, 1122 wantiifac: []int{49, 2, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1123 }, 1124 { 1125 n: 32, 1126 1127 wantiwork: []float64{ 1128 0.000000, 0.1960343, 0.3901806, 0.5805693, 0.7653669, 0.9427935, 1.111140, 1.268787, 1129 1.414214, 1.546021, 1.662939, 1.763843, 1.847759, 1.913881, 1.961571, 1.990369, 1130 0.000000, 0.1960343, 0.3901805, 0.5805693, 0.7653669, 0.9427933, 1.111140, 1.268787, 1131 1.414214, 1.546021, 1.662939, 1.763842, 1.847759, 1.913881, 1.961571, 1.990369, 1132 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1133 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1134 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1135 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1136 0.000000, 0.9807853, 0.1950903, 0.9238795, 0.3826835, 0.8314696, 0.5555702, 0.7071068, 1137 0.7071068, 0.5555702, 0.8314697, 0.3826834, 0.9238795, 0.1950902, 0.9807853, 0.000000, 1138 0.000000, 0.9238795, 0.3826835, 0.000000, 0.000000, 0.7071068, 0.7071068, 0.000000, 1139 0.000000, 0.3826834, 0.9238795, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1140 0.000000, 0.4484155e-43, 0.4203895e-44, 1141 }, 1142 wantiifac: []int{32, 3, 2, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1143 }, 1144 { 1145 n: 25, 1146 1147 wantiwork: []float64{ 1148 0.000000, 0.2506665, 0.4973798, 0.7362491, 0.9635074, 1.175570, 1.369094, 1.541027, 1149 1.688656, 1.809654, 1.902113, 1.964575, 1.996053, 0.1255808, 0.3747624, 0.6180339, 1150 0.8515584, 1.071653, 1.274848, 1.457937, 1.618034, 1.752613, 1.859553, 1.937166, 1151 1.984229, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1152 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1153 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1154 0.000000, 0.000000, 0.000000, 0.9685832, 0.2486899, 0.8763067, 0.4817537, 0.000000, 1155 0.8763067, 0.4817537, 0.5358267, 0.8443279, 0.000000, 0.7289686, 0.6845471, 0.6279038e-01, 1156 0.9980267, 0.000000, 0.5358267, 0.8443279, -0.4257794, 0.9048270, 0.000000, 0.000000, 1157 0.000000, 0.000000, 0.000000, 0.000000, 0.3503246e-43, 0.2802597e-44, 1158 }, 1159 wantiifac: []int{25, 2, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1160 }, 1161 { 1162 n: 4, 1163 1164 wantiwork: []float64{ 1165 0.000000, 1.414214, 0.000000, 1.414214, 0.000000, 0.000000, 0.000000, 0.000000, 1166 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.5605194e-44, 0.1401298e-44, 1167 }, 1168 wantiifac: []int{4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1169 }, 1170 { 1171 n: 3, 1172 1173 wantiwork: []float64{ 1174 0.000000, 1.732051, 1.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1175 0.000000, 0.000000, 0.4203895e-44, 0.1401298e-44, 1176 }, 1177 wantiifac: []int{3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1178 }, 1179 { 1180 n: 2, 1181 1182 wantiwork: []float64{ 1183 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1184 0.000000, 1185 }, 1186 wantiifac: []int{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1187 }, 1188 } 1189 1190 func TestCosq(t *testing.T) { 1191 t.Parallel() 1192 const tol = 1e-12 1193 for _, test := range sincosqTests { 1194 // Compute the work and factor slices and compare to known values. 1195 work := make([]float64, 3*test.n) 1196 ifac := make([]int, 15) 1197 Cosqi(test.n, work, ifac) 1198 var failed bool 1199 if !floats.EqualApprox(work, test.wantiwork, 1e-6) { 1200 failed = true 1201 t.Errorf("unexpected work after call to cosqi for n=%d", test.n) 1202 t.Logf("\n%v\n%v", work, test.wantiwork) 1203 } 1204 if !reflect.DeepEqual(ifac, test.wantiifac) { 1205 failed = true 1206 t.Errorf("unexpected ifac after call to cosqi for n=%d", test.n) 1207 } 1208 if failed { 1209 continue 1210 } 1211 1212 // Construct a sequence with known a frequency spectrum and compare 1213 // the computed spectrum. 1214 fn := float64(test.n) 1215 x := make([]float64, test.n) 1216 y, _, xh := series(test.n) 1217 1218 dt := math.Pi / (2 * fn) 1219 for i := 0; i < test.n; i++ { 1220 arg := float64(i) * dt 1221 for k := 0; k < test.n; k++ { 1222 x[i] += y[k] * math.Cos(float64(2*k+1)*arg) 1223 } 1224 x[i] *= 4 1225 } 1226 1227 Cosqb(test.n, y, work, ifac) 1228 1229 cf := 0.25 / fn 1230 var cosqbt float64 1231 for i := 0; i < test.n; i++ { 1232 cosqbt = math.Max(cosqbt, math.Abs(x[i]-y[i])) 1233 x[i] = xh[i] 1234 } 1235 cosqbt *= cf 1236 if !scalar.EqualWithinAbsOrRel(cosqbt, 0, tol, tol) { 1237 t.Errorf("unexpected cosqbt value for n=%d: got:%f want:0", test.n, cosqbt) 1238 } 1239 1240 // Construct a frequency spectrum and compare the computed sequence. 1241 for i := 0; i < test.n; i++ { 1242 y[i] = 0.5 * x[0] 1243 arg := float64(2*i+1) * dt 1244 for k := 1; k < test.n; k++ { 1245 y[i] += x[k] * math.Cos(float64(k)*arg) 1246 } 1247 y[i] *= 2 1248 } 1249 1250 Cosqf(test.n, x, work, ifac) 1251 1252 var cosqft float64 1253 for i := 0; i < test.n; i++ { 1254 cosqft = math.Max(cosqft, math.Abs(x[i]-y[i])) 1255 x[i] = xh[i] 1256 y[i] = xh[i] 1257 } 1258 cosqft *= cf 1259 if !scalar.EqualWithinAbsOrRel(cosqft, 0, tol, tol) { 1260 t.Errorf("unexpected cosqft value for n=%d: got:%f want:0", test.n, cosqft) 1261 } 1262 1263 // Check that Cosqb and Cosqf are inverses. 1264 Cosqb(test.n, x, work, ifac) 1265 Cosqf(test.n, x, work, ifac) 1266 1267 var cosqfb float64 1268 for i := 0; i < test.n; i++ { 1269 cosqfb = math.Max(cosqfb, math.Abs(cf*x[i]-y[i])) 1270 } 1271 if !scalar.EqualWithinAbsOrRel(cosqfb, 0, tol, tol) { 1272 t.Errorf("unexpected cosqfb value for n=%d: got:%f want:0", test.n, cosqfb) 1273 } 1274 } 1275 } 1276 1277 func TestSinq(t *testing.T) { 1278 t.Parallel() 1279 const tol = 1e-12 1280 for _, test := range sincosqTests { 1281 // Compute the work and factor slices and compare to known values. 1282 work := make([]float64, 3*test.n) 1283 ifac := make([]int, 15) 1284 Sinqi(test.n, work, ifac) 1285 var failed bool 1286 if !floats.EqualApprox(work, test.wantiwork, 1e-6) { 1287 failed = true 1288 t.Errorf("unexpected work after call to sinqi for n=%d", test.n) 1289 t.Logf("\n%v\n%v", work, test.wantiwork) 1290 } 1291 if !reflect.DeepEqual(ifac, test.wantiifac) { 1292 failed = true 1293 t.Errorf("unexpected ifac after call to sinqi for n=%d", test.n) 1294 } 1295 if failed { 1296 continue 1297 } 1298 1299 // Construct a sequence with known a frequency spectrum and compare 1300 // the computed spectrum. 1301 fn := float64(test.n) 1302 x := make([]float64, test.n) 1303 y, _, xh := series(test.n) 1304 1305 dt := math.Pi / (2 * fn) 1306 for i := 0; i < test.n; i++ { 1307 arg := float64(i+1) * dt 1308 for k := 0; k < test.n; k++ { 1309 x[i] += y[k] * math.Sin(float64(2*k+1)*arg) 1310 } 1311 x[i] *= 4 1312 } 1313 1314 Sinqb(test.n, y, work, ifac) 1315 1316 cf := 0.25 / fn 1317 var sinqbt float64 1318 for i := 0; i < test.n; i++ { 1319 sinqbt = math.Max(sinqbt, math.Abs(x[i]-y[i])) 1320 x[i] = xh[i] 1321 } 1322 sinqbt *= cf 1323 if !scalar.EqualWithinAbsOrRel(sinqbt, 0, tol, tol) { 1324 t.Errorf("unexpected sinqbt value for n=%d: got:%f want:0", test.n, sinqbt) 1325 } 1326 1327 // Construct a frequency spectrum and compare the computed sequence. 1328 for i := 0; i < test.n; i++ { 1329 arg := float64(2*i+1) * dt 1330 y[i] = 0.5 * math.Pow(-1, float64(i)) * x[test.n-1] 1331 for k := 0; k < test.n-1; k++ { 1332 y[i] += x[k] * math.Sin(float64(k+1)*arg) 1333 } 1334 y[i] *= 2 1335 } 1336 1337 Sinqf(test.n, x, work, ifac) 1338 1339 var sinqft float64 1340 for i := 0; i < test.n; i++ { 1341 sinqft = math.Max(sinqft, math.Abs(x[i]-y[i])) 1342 x[i] = xh[i] 1343 y[i] = xh[i] 1344 } 1345 sinqft *= cf 1346 if !scalar.EqualWithinAbsOrRel(sinqft, 0, tol, tol) { 1347 t.Errorf("unexpected sinqft value for n=%d: got:%f want:0", test.n, sinqft) 1348 } 1349 1350 // Check that Sinqb and Sinqf are inverses. 1351 Sinqb(test.n, x, work, ifac) 1352 Sinqf(test.n, x, work, ifac) 1353 1354 var sinqfb float64 1355 for i := 0; i < test.n; i++ { 1356 sinqfb = math.Max(sinqfb, math.Abs(cf*x[i]-y[i])) 1357 } 1358 if !scalar.EqualWithinAbsOrRel(sinqfb, 0, tol, tol) { 1359 t.Errorf("unexpected sinqfb value for n=%d: got:%f want:0", test.n, sinqfb) 1360 } 1361 } 1362 } 1363 1364 var sincosqTests = []struct { 1365 n int 1366 1367 // The following two fields are added as there is no unit testing in 1368 // FFTPACK for SINTI. 1369 // 1370 // wantiwork is obtained from the FFTPACK test.f with modification. 1371 // The W array is zeroed at each iteration and the first 3n elements 1372 // of W are printed after the call to RFFTI. 1373 wantiwork []float64 1374 // wantiifac is obtained from the FFTPACK sint.f with modification. 1375 // The IFAC array is zeroed at each iteration of test.f and the 15 elements 1376 // of IFAC are printed before RFFTI returns. 1377 wantiifac []int 1378 }{ 1379 { 1380 n: 120, 1381 1382 wantiwork: []float64{ 1383 0.9999143, 0.9996573, 0.9992290, 0.9986295, 0.9978589, 0.9969173, 0.9958049, 0.9945219, 1384 0.9930685, 0.9914449, 0.9896514, 0.9876884, 0.9855561, 0.9832549, 0.9807853, 0.9781476, 1385 0.9753423, 0.9723699, 0.9692309, 0.9659258, 0.9624552, 0.9588197, 0.9550200, 0.9510565, 1386 0.9469301, 0.9426415, 0.9381914, 0.9335804, 0.9288096, 0.9238795, 0.9187912, 0.9135454, 1387 0.9081432, 0.9025853, 0.8968728, 0.8910065, 0.8849877, 0.8788171, 0.8724960, 0.8660254, 1388 0.8594064, 0.8526402, 0.8457278, 0.8386706, 0.8314696, 0.8241262, 0.8166415, 0.8090170, 1389 0.8012538, 0.7933533, 0.7853169, 0.7771459, 0.7688418, 0.7604060, 0.7518398, 0.7431448, 1390 0.7343225, 0.7253744, 0.7163019, 0.7071068, 0.6977904, 0.6883546, 0.6788007, 0.6691306, 1391 0.6593458, 0.6494480, 0.6394390, 0.6293204, 0.6190940, 0.6087614, 0.5983246, 0.5877852, 1392 0.5771452, 0.5664062, 0.5555702, 0.5446390, 0.5336145, 0.5224985, 0.5112931, 0.5000000, 1393 0.4886212, 0.4771588, 0.4656145, 0.4539905, 0.4422887, 0.4305110, 0.4186597, 0.4067366, 1394 0.3947438, 0.3826834, 0.3705574, 0.3583679, 0.3461170, 0.3338068, 0.3214395, 0.3090170, 1395 0.2965415, 0.2840153, 0.2714404, 0.2588191, 0.2461533, 0.2334453, 0.2206974, 0.2079117, 1396 0.1950902, 0.1822355, 0.1693494, 0.1564344, 0.1434926, 0.1305261, 0.1175374, 0.1045284, 1397 0.9150153e-01, 0.7845908e-01, 0.6540307e-01, 0.5233597e-01, 0.3925979e-01, 0.2617688e-01, 0.1308960e-01, -0.4371139e-07, 1398 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1399 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1400 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1401 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1402 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1403 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1404 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1405 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1406 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1407 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1408 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1409 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1410 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1411 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1412 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1413 0.9986295, 0.5233596e-01, 0.9945219, 0.1045285, 0.9876884, 0.1564345, 0.9781476, 0.2079117, 1414 0.9659258, 0.2588190, 0.9510565, 0.3090170, 0.9335804, 0.3583679, 0.9135454, 0.4067366, 1415 0.8910065, 0.4539905, 0.8660254, 0.5000000, 0.8386706, 0.5446391, 0.8090170, 0.5877852, 1416 0.7771459, 0.6293204, 0.7431448, 0.6691306, 0.7071068, 0.7071068, 0.6691306, 0.7431449, 1417 0.6293204, 0.7771460, 0.5877852, 0.8090170, 0.5446390, 0.8386706, 0.5000000, 0.8660254, 1418 0.4539905, 0.8910065, 0.4067366, 0.9135455, 0.3583679, 0.9335805, 0.3090170, 0.9510565, 1419 0.2588191, 0.9659258, 0.2079117, 0.9781476, 0.1564344, 0.9876884, 0.1045284, 0.9945219, 1420 0.5233597e-01, 0.9986295, 0.000000, 0.000000, 0.9945219, 0.1045285, 0.9781476, 0.2079117, 1421 0.9510565, 0.3090170, 0.9135454, 0.4067366, 0.8660254, 0.5000000, 0.8090170, 0.5877852, 1422 0.7431448, 0.6691306, 0.000000, 0.9781476, 0.2079117, 0.9135454, 0.4067366, 0.8090170, 1423 0.5877852, 0.6691306, 0.7431449, 0.5000000, 0.8660254, 0.3090170, 0.9510565, 0.1045284, 1424 0.9945219, 0.000000, 0.9510565, 0.3090170, 0.8090170, 0.5877852, 0.5877852, 0.8090170, 1425 0.3090170, 0.9510565, -0.4371139e-07, 1.000000, -0.3090170, 0.9510565, -0.5877852, 0.8090170, 1426 0.000000, 0.9135454, 0.4067366, 0.6691306, 0.7431449, 0.000000, 0.6691306, 0.7431449, 1427 -0.1045285, 0.9945219, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1428 }, 1429 wantiifac: []int{120, 4, 2, 4, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1430 }, 1431 { 1432 n: 54, 1433 1434 wantiwork: []float64{ 1435 0.9995769, 0.9983082, 0.9961947, 0.9932383, 0.9894416, 0.9848077, 0.9793406, 0.9730449, 1436 0.9659258, 0.9579895, 0.9492427, 0.9396926, 0.9293475, 0.9182161, 0.9063078, 0.8936327, 1437 0.8802014, 0.8660254, 0.8511167, 0.8354878, 0.8191521, 0.8021232, 0.7844157, 0.7660444, 1438 0.7470251, 0.7273737, 0.7071068, 0.6862416, 0.6647958, 0.6427876, 0.6202354, 0.5971586, 1439 0.5735765, 0.5495090, 0.5249766, 0.5000000, 0.4746003, 0.4487992, 0.4226182, 0.3960797, 1440 0.3692062, 0.3420202, 0.3145447, 0.2868032, 0.2588191, 0.2306159, 0.2022175, 0.1736482, 1441 0.1449319, 0.1160929, 0.8715568e-01, 0.5814485e-01, 0.2908471e-01, -0.4371139e-07, 0.000000, 0.000000, 1442 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1443 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1444 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1445 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1446 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1447 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1448 0.000000, 0.000000, 0.000000, 0.000000, 0.9932383, 0.1160929, 0.9730449, 0.2306159, 1449 0.9396926, 0.3420201, 0.8936327, 0.4487992, 0.8354878, 0.5495090, 0.7660444, 0.6427876, 1450 0.6862416, 0.7273737, 0.5971586, 0.8021232, 0.5000000, 0.8660254, 0.3960797, 0.9182161, 1451 0.2868032, 0.9579895, 0.1736482, 0.9848077, 0.5814485e-01, 0.9983082, 0.000000, 0.9730449, 1452 0.2306159, 0.8936327, 0.4487992, 0.7660444, 0.6427876, 0.5971586, 0.8021232, 0.000000, 1453 0.8936327, 0.4487992, 0.5971586, 0.8021232, 0.1736482, 0.9848077, -0.2868032, 0.9579895, 1454 0.000000, 0.7660444, 0.6427876, 0.000000, 0.1736482, 0.9848077, 0.000000, 0.000000, 1455 0.000000, 0.000000, 1456 }, 1457 wantiifac: []int{54, 4, 2, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1458 }, 1459 { 1460 n: 49, 1461 1462 wantiwork: []float64{ 1463 0.9994862, 0.9979454, 0.9953791, 0.9917900, 0.9871818, 0.9815592, 0.9749279, 0.9672949, 1464 0.9586679, 0.9490557, 0.9384684, 0.9269168, 0.9144126, 0.9009688, 0.8865993, 0.8713187, 1465 0.8551428, 0.8380881, 0.8201723, 0.8014136, 0.7818314, 0.7614459, 0.7402779, 0.7183493, 1466 0.6956825, 0.6723009, 0.6482283, 0.6234898, 0.5981105, 0.5721166, 0.5455348, 0.5183925, 1467 0.4907175, 0.4625383, 0.4338836, 0.4047833, 0.3752669, 0.3453650, 0.3151082, 0.2845275, 1468 0.2536545, 0.2225209, 0.1911586, 0.1595999, 0.1278771, 0.9602292e-01, 0.6407014e-01, 0.3205151e-01, 1469 -0.4371139e-07, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1470 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1471 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1472 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1473 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1474 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1475 0.000000, 0.000000, 0.9917900, 0.1278772, 0.9672949, 0.2536546, 0.9269168, 0.3752670, 1476 0.000000, 0.9672949, 0.2536546, 0.8713187, 0.4907176, 0.7183493, 0.6956826, 0.000000, 1477 0.9269168, 0.3752670, 0.7183493, 0.6956826, 0.4047833, 0.9144127, 0.000000, 0.8713187, 1478 0.4907176, 0.5183925, 0.8551428, 0.3205151e-01, 0.9994862, 0.000000, 0.8014136, 0.5981106, 1479 0.2845275, 0.9586679, -0.3453652, 0.9384683, 0.000000, 0.7183493, 0.6956826, 0.3205151e-01, 1480 0.9994862, -0.6723010, 0.7402779, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1481 0.000000, 0.000000, 0.000000, 1482 }, 1483 wantiifac: []int{49, 2, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1484 }, 1485 { 1486 n: 32, 1487 1488 wantiwork: []float64{ 1489 0.9987954, 0.9951847, 0.9891765, 0.9807853, 0.9700313, 0.9569404, 0.9415441, 0.9238795, 1490 0.9039893, 0.8819212, 0.8577286, 0.8314696, 0.8032075, 0.7730104, 0.7409511, 0.7071068, 1491 0.6715589, 0.6343933, 0.5956993, 0.5555702, 0.5141027, 0.4713967, 0.4275551, 0.3826834, 1492 0.3368898, 0.2902846, 0.2429801, 0.1950902, 0.1467305, 0.9801713e-01, 0.4906765e-01, -0.4371139e-07, 1493 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1494 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1495 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1496 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1497 0.9807853, 0.1950903, 0.9238795, 0.3826835, 0.8314696, 0.5555702, 0.7071068, 0.7071068, 1498 0.5555702, 0.8314697, 0.3826834, 0.9238795, 0.1950902, 0.9807853, 0.000000, 0.000000, 1499 0.9238795, 0.3826835, 0.000000, 0.000000, 0.7071068, 0.7071068, 0.000000, 0.000000, 1500 0.3826834, 0.9238795, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1501 }, 1502 wantiifac: []int{32, 3, 2, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1503 }, 1504 { 1505 n: 25, 1506 1507 wantiwork: []float64{ 1508 0.9980267, 0.9921147, 0.9822872, 0.9685832, 0.9510565, 0.9297765, 0.9048271, 0.8763067, 1509 0.8443279, 0.8090170, 0.7705132, 0.7289686, 0.6845471, 0.6374239, 0.5877852, 0.5358267, 1510 0.4817536, 0.4257792, 0.3681245, 0.3090170, 0.2486898, 0.1873812, 0.1253331, 0.6279038e-01, 1511 -0.4371139e-07, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1512 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1513 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1514 0.000000, 0.000000, 0.9685832, 0.2486899, 0.8763067, 0.4817537, 0.000000, 0.8763067, 1515 0.4817537, 0.5358267, 0.8443279, 0.000000, 0.7289686, 0.6845471, 0.6279038e-01, 0.9980267, 1516 0.000000, 0.5358267, 0.8443279, -0.4257794, 0.9048270, 0.000000, 0.000000, 0.000000, 1517 0.000000, 0.000000, 0.000000, 1518 }, 1519 wantiifac: []int{25, 2, 5, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1520 }, 1521 { 1522 n: 4, 1523 1524 wantiwork: []float64{ 1525 0.9238795, 0.7071068, 0.3826834, -0.4371139e-07, 0.000000, 0.000000, 0.000000, 0.000000, 1526 0.000000, 0.000000, 0.000000, 0.000000, 1527 }, 1528 wantiifac: []int{4, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1529 }, 1530 { 1531 n: 3, 1532 1533 wantiwork: []float64{ 1534 0.8660254, 0.5000000, -0.4371139e-07, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 1535 0.000000, 1536 }, 1537 wantiifac: []int{3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1538 }, 1539 { 1540 n: 2, 1541 1542 wantiwork: []float64{ 1543 0.7071068, -0.4371139e-07, 0.000000, 0.000000, 0.000000, 0.000000, 1544 }, 1545 wantiifac: []int{2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 1546 }, 1547 } 1548 1549 // series returns three copies of a sinusoidal sequence n samples long. 1550 func series(n int) (x, y, xh []float64) { 1551 x = make([]float64, n) 1552 y = make([]float64, n) 1553 xh = make([]float64, n) 1554 for i := 0; i < n; i++ { 1555 x[i] = math.Sin(float64(i+1) * math.Sqrt2) 1556 y[i] = x[i] 1557 xh[i] = x[i] 1558 } 1559 return x, y, xh 1560 }