-
Notifications
You must be signed in to change notification settings - Fork 58
/
test_fft_factors.c
142 lines (111 loc) · 3.96 KB
/
test_fft_factors.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#ifdef PFFFT_ENABLE_FLOAT
#include "pffft.h"
#endif
#ifdef PFFFT_ENABLE_DOUBLE
#include "pffft_double.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#ifdef PFFFT_ENABLE_FLOAT
int test_float(int TL)
{
PFFFT_Setup * S;
for (int dir_i = 0; dir_i <= 1; ++dir_i)
{
for (int cplx_i = 0; cplx_i <= 1; ++cplx_i)
{
const pffft_direction_t dir = (!dir_i) ? PFFFT_FORWARD : PFFFT_BACKWARD;
const pffft_transform_t cplx = (!cplx_i) ? PFFFT_REAL : PFFFT_COMPLEX;
const int N_min = pffft_min_fft_size(cplx);
const int N_max = N_min * 11 + N_min;
int NTL = pffft_nearest_transform_size(TL, cplx, (!dir_i));
double near_off = (NTL - TL) * 100.0 / (double)TL;
fprintf(stderr, "testing float, %s, %s ..\tminimum transform %d; nearest transform for %d is %d (%.2f%% off)\n",
(!dir_i) ? "FORWARD" : "BACKWARD", (!cplx_i) ? "REAL" : "COMPLEX", N_min, TL, NTL, near_off );
for (int N = (N_min/2); N <= N_max; N += (N_min/2))
{
int R = N, f2 = 0, f3 = 0, f5 = 0, tmp_f;
const int factorizable = pffft_is_valid_size(N, cplx);
while (R >= 5*N_min && (R % 5) == 0) { R /= 5; ++f5; }
while (R >= 3*N_min && (R % 3) == 0) { R /= 3; ++f3; }
while (R >= 2*N_min && (R % 2) == 0) { R /= 2; ++f2; }
tmp_f = (R == N_min) ? 1 : 0;
assert( factorizable == tmp_f );
S = pffft_new_setup(N, cplx);
if ( S && !factorizable )
{
fprintf(stderr, "fft setup successful, but NOT factorizable into min(=%d), 2^%d, 3^%d, 5^%d for N = %d (R = %d)\n", N_min, f2, f3, f5, N, R);
return 1;
}
else if ( !S && factorizable)
{
fprintf(stderr, "fft setup UNsuccessful, but factorizable into min(=%d), 2^%d, 3^%d, 5^%d for N = %d (R = %d)\n", N_min, f2, f3, f5, N, R);
return 1;
}
if (S)
pffft_destroy_setup(S);
}
}
}
return 0;
}
#endif
#ifdef PFFFT_ENABLE_DOUBLE
int test_double(int TL)
{
PFFFTD_Setup * S;
for (int dir_i = 0; dir_i <= 1; ++dir_i)
{
for (int cplx_i = 0; cplx_i <= 1; ++cplx_i)
{
const pffft_direction_t dir = (!dir_i) ? PFFFT_FORWARD : PFFFT_BACKWARD;
const pffft_transform_t cplx = (!cplx_i) ? PFFFT_REAL : PFFFT_COMPLEX;
const int N_min = pffftd_min_fft_size(cplx);
const int N_max = N_min * 11 + N_min;
int NTL = pffftd_nearest_transform_size(TL, cplx, (!dir_i));
double near_off = (NTL - TL) * 100.0 / (double)TL;
fprintf(stderr, "testing double, %s, %s ..\tminimum transform %d; nearest transform for %d is %d (%.2f%% off)\n",
(!dir_i) ? "FORWARD" : "BACKWARD", (!cplx_i) ? "REAL" : "COMPLEX", N_min, TL, NTL, near_off );
for (int N = (N_min/2); N <= N_max; N += (N_min/2))
{
int R = N, f2 = 0, f3 = 0, f5 = 0, tmp_f;
const int factorizable = pffftd_is_valid_size(N, cplx);
while (R >= 5*N_min && (R % 5) == 0) { R /= 5; ++f5; }
while (R >= 3*N_min && (R % 3) == 0) { R /= 3; ++f3; }
while (R >= 2*N_min && (R % 2) == 0) { R /= 2; ++f2; }
tmp_f = (R == N_min) ? 1 : 0;
assert( factorizable == tmp_f );
S = pffftd_new_setup(N, cplx);
if ( S && !factorizable )
{
fprintf(stderr, "fft setup successful, but NOT factorizable into min(=%d), 2^%d, 3^%d, 5^%d for N = %d (R = %d)\n", N_min, f2, f3, f5, N, R);
return 1;
}
else if ( !S && factorizable)
{
fprintf(stderr, "fft setup UNsuccessful, but factorizable into min(=%d), 2^%d, 3^%d, 5^%d for N = %d (R = %d)\n", N_min, f2, f3, f5, N, R);
return 1;
}
if (S)
pffftd_destroy_setup(S);
}
}
}
return 0;
}
#endif
int main(int argc, char *argv[])
{
int N = (1 < argc) ? atoi(argv[1]) : 2;
int r = 0;
#ifdef PFFFT_ENABLE_FLOAT
r = test_float(N);
if (r)
return r;
#endif
#ifdef PFFFT_ENABLE_DOUBLE
r = test_double(N);
#endif
return r;
}