Coverage Report

Created: 2025-08-29 06:11

/proc/self/cwd/libfaad/cfft.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3
** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com
4
**
5
** This program is free software; you can redistribute it and/or modify
6
** it under the terms of the GNU General Public License as published by
7
** the Free Software Foundation; either version 2 of the License, or
8
** (at your option) any later version.
9
**
10
** This program is distributed in the hope that it will be useful,
11
** but WITHOUT ANY WARRANTY; without even the implied warranty of
12
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
** GNU General Public License for more details.
14
**
15
** You should have received a copy of the GNU General Public License
16
** along with this program; if not, write to the Free Software
17
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18
**
19
** Any non-GPL usage of this software or parts of this software is strictly
20
** forbidden.
21
**
22
** The "appropriate copyright message" mentioned in section 2c of the GPLv2
23
** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"
24
**
25
** Commercial non-GPL licensing of this software is possible.
26
** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
27
**
28
** $Id: cfft.c,v 1.35 2007/11/01 12:33:29 menno Exp $
29
**/
30
31
/*
32
 * Algorithmically based on Fortran-77 FFTPACK
33
 * by Paul N. Swarztrauber(Version 4, 1985).
34
 *
35
 * Does even sized fft only
36
 */
37
38
/* isign is +1 for backward and -1 for forward transforms */
39
40
#include "common.h"
41
#include "structs.h"
42
43
#include <stdlib.h>
44
45
#include "cfft.h"
46
#include "cfft_tab.h"
47
48
49
/* static function declarations */
50
static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
51
                      complex_t *ch, const complex_t *wa);
52
static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
53
                      complex_t *ch, const complex_t *wa);
54
static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
55
                   complex_t *ch, const complex_t *wa1, const complex_t *wa2, const int8_t isign);
56
static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
57
                      const complex_t *wa1, const complex_t *wa2, const complex_t *wa3);
58
static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
59
                      const complex_t *wa1, const complex_t *wa2, const complex_t *wa3);
60
static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc, complex_t *ch,
61
                   const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
62
                   const complex_t *wa4, const int8_t isign);
63
static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac);
64
65
66
/*----------------------------------------------------------------------
67
   passf2, passf3, passf4, passf5. Complex FFT passes fwd and bwd.
68
  ----------------------------------------------------------------------*/
69
70
static void passf2pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
71
                      complex_t *ch, const complex_t *wa)
72
533k
{
73
533k
    uint16_t i, k, ah, ac;
74
75
533k
    if (ido == 1)
76
0
    {
77
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
78
        // TODO: remove this code once fuzzer proves it is totally unreahable
79
        // For supported frame lengths that have odd number of factor 2 it is
80
        // never the last factor; consequently `ido` should never be 1.
81
0
        __builtin_trap();
82
        /*
83
#endif
84
        for (k = 0; k < l1; k++)
85
        {
86
            ah = 2*k;
87
            ac = 4*k;
88
89
            RE(ch[ah])    = RE(cc[ac]) + RE(cc[ac+1]);
90
            RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
91
            IM(ch[ah])    = IM(cc[ac]) + IM(cc[ac+1]);
92
            IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
93
        }
94
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
95
        */
96
0
#endif
97
533k
    } else {
98
1.06M
        for (k = 0; k < l1; k++)
99
533k
        {
100
533k
            ah = k*ido;
101
533k
            ac = 2*k*ido;
102
103
135M
            for (i = 0; i < ido; i++)
104
135M
            {
105
135M
                complex_t t2;
106
107
135M
                RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
108
135M
                RE(t2)       = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
109
110
135M
                IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
111
135M
                IM(t2)       = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
112
113
135M
#if 1
114
135M
                ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
115
135M
                    IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
116
#else
117
                ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
118
                    RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
119
#endif
120
135M
            }
121
533k
        }
122
533k
    }
123
533k
}
124
125
static void passf2neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
126
                      complex_t *ch, const complex_t *wa)
127
12.6k
{
128
12.6k
    uint16_t i, k, ah, ac;
129
130
12.6k
    if (ido == 1)
131
0
    {
132
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
133
        // TODO: remove this code once fuzzer proves it is totally unreahable
134
        // For supported frame lengths that have odd number of factor 2 it is
135
        // never the last factor; consequently `ido` should never be 1.
136
0
        __builtin_trap();
137
        /*
138
#endif
139
        for (k = 0; k < l1; k++)
140
        {
141
            ah = 2*k;
142
            ac = 4*k;
143
144
            RE(ch[ah])    = RE(cc[ac]) + RE(cc[ac+1]);
145
            RE(ch[ah+l1]) = RE(cc[ac]) - RE(cc[ac+1]);
146
            IM(ch[ah])    = IM(cc[ac]) + IM(cc[ac+1]);
147
            IM(ch[ah+l1]) = IM(cc[ac]) - IM(cc[ac+1]);
148
        }
149
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
150
        */
151
0
#endif
152
12.6k
    } else {
153
25.3k
        for (k = 0; k < l1; k++)
154
12.6k
        {
155
12.6k
            ah = k*ido;
156
12.6k
            ac = 2*k*ido;
157
158
3.15M
            for (i = 0; i < ido; i++)
159
3.13M
            {
160
3.13M
                complex_t t2;
161
162
3.13M
                RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
163
3.13M
                RE(t2)       = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
164
165
3.13M
                IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
166
3.13M
                IM(t2)       = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
167
168
3.13M
#if 1
169
3.13M
                ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
170
3.13M
                    RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
171
#else
172
                ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
173
                    IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
174
#endif
175
3.13M
            }
176
12.6k
        }
177
12.6k
    }
178
12.6k
}
179
180
181
static void passf3(const uint16_t ido, const uint16_t l1, const complex_t *cc,
182
                   complex_t *ch, const complex_t *wa1, const complex_t *wa2,
183
                   const int8_t isign)
184
201k
{
185
201k
    static real_t taur = FRAC_CONST(-0.5);
186
201k
    static real_t taui = FRAC_CONST(0.866025403784439);
187
201k
    uint16_t i, k, ac, ah;
188
201k
    complex_t c2, c3, d2, d3, t2;
189
190
201k
    if (ido == 1)
191
0
    {
192
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
193
        // TODO: remove this code once fuzzer proves it is totally unreahable
194
        // 3 is never the the biggest factor for supported frame lengths;
195
        // consequently `ido` should never be 1.
196
0
        __builtin_trap();
197
        /*
198
#endif
199
        if (isign == 1)
200
        {
201
            for (k = 0; k < l1; k++)
202
            {
203
                ac = 3*k+1;
204
                ah = k;
205
206
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
207
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
208
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
209
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
210
211
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
212
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
213
214
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
215
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
216
217
                RE(ch[ah+l1]) = RE(c2) - IM(c3);
218
                IM(ch[ah+l1]) = IM(c2) + RE(c3);
219
                RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
220
                IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
221
            }
222
        } else {
223
            for (k = 0; k < l1; k++)
224
            {
225
                ac = 3*k+1;
226
                ah = k;
227
228
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
229
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
230
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
231
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
232
233
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
234
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
235
236
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
237
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
238
239
                RE(ch[ah+l1]) = RE(c2) + IM(c3);
240
                IM(ch[ah+l1]) = IM(c2) - RE(c3);
241
                RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
242
                IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
243
            }
244
        }
245
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
246
        */
247
0
#endif
248
201k
    } else {
249
201k
        if (isign == 1)
250
194k
        {
251
477k
            for (k = 0; k < l1; k++)
252
283k
            {
253
16.7M
                for (i = 0; i < ido; i++)
254
16.4M
                {
255
16.4M
                    ac = i + (3*k+1)*ido;
256
16.4M
                    ah = i + k * ido;
257
258
16.4M
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
259
16.4M
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
260
16.4M
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
261
16.4M
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
262
263
16.4M
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
264
16.4M
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
265
266
16.4M
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
267
16.4M
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
268
269
16.4M
                    RE(d2) = RE(c2) - IM(c3);
270
16.4M
                    IM(d3) = IM(c2) - RE(c3);
271
16.4M
                    RE(d3) = RE(c2) + IM(c3);
272
16.4M
                    IM(d2) = IM(c2) + RE(c3);
273
274
16.4M
#if 1
275
16.4M
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
276
16.4M
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
277
16.4M
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
278
16.4M
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
279
#else
280
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
281
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
282
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
283
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
284
#endif
285
16.4M
                }
286
283k
            }
287
194k
        } else {
288
20.7k
            for (k = 0; k < l1; k++)
289
13.8k
            {
290
1.11M
                for (i = 0; i < ido; i++)
291
1.10M
                {
292
1.10M
                    ac = i + (3*k+1)*ido;
293
1.10M
                    ah = i + k * ido;
294
295
1.10M
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
296
1.10M
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
297
1.10M
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
298
1.10M
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
299
300
1.10M
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
301
1.10M
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
302
303
1.10M
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
304
1.10M
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
305
306
1.10M
                    RE(d2) = RE(c2) + IM(c3);
307
1.10M
                    IM(d3) = IM(c2) + RE(c3);
308
1.10M
                    RE(d3) = RE(c2) - IM(c3);
309
1.10M
                    IM(d2) = IM(c2) - RE(c3);
310
311
1.10M
#if 1
312
1.10M
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
313
1.10M
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
314
1.10M
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
315
1.10M
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
316
#else
317
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
318
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
319
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
320
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
321
#endif
322
1.10M
                }
323
13.8k
            }
324
6.95k
        }
325
201k
    }
326
201k
}
cfft.c:passf3
Line
Count
Source
184
79.3k
{
185
79.3k
    static real_t taur = FRAC_CONST(-0.5);
186
79.3k
    static real_t taui = FRAC_CONST(0.866025403784439);
187
79.3k
    uint16_t i, k, ac, ah;
188
79.3k
    complex_t c2, c3, d2, d3, t2;
189
190
79.3k
    if (ido == 1)
191
0
    {
192
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
193
        // TODO: remove this code once fuzzer proves it is totally unreahable
194
        // 3 is never the the biggest factor for supported frame lengths;
195
        // consequently `ido` should never be 1.
196
0
        __builtin_trap();
197
        /*
198
#endif
199
        if (isign == 1)
200
        {
201
            for (k = 0; k < l1; k++)
202
            {
203
                ac = 3*k+1;
204
                ah = k;
205
206
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
207
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
208
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
209
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
210
211
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
212
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
213
214
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
215
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
216
217
                RE(ch[ah+l1]) = RE(c2) - IM(c3);
218
                IM(ch[ah+l1]) = IM(c2) + RE(c3);
219
                RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
220
                IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
221
            }
222
        } else {
223
            for (k = 0; k < l1; k++)
224
            {
225
                ac = 3*k+1;
226
                ah = k;
227
228
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
229
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
230
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
231
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
232
233
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
234
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
235
236
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
237
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
238
239
                RE(ch[ah+l1]) = RE(c2) + IM(c3);
240
                IM(ch[ah+l1]) = IM(c2) - RE(c3);
241
                RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
242
                IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
243
            }
244
        }
245
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
246
        */
247
0
#endif
248
79.3k
    } else {
249
79.3k
        if (isign == 1)
250
76.2k
        {
251
184k
            for (k = 0; k < l1; k++)
252
108k
            {
253
6.15M
                for (i = 0; i < ido; i++)
254
6.04M
                {
255
6.04M
                    ac = i + (3*k+1)*ido;
256
6.04M
                    ah = i + k * ido;
257
258
6.04M
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
259
6.04M
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
260
6.04M
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
261
6.04M
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
262
263
6.04M
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
264
6.04M
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
265
266
6.04M
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
267
6.04M
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
268
269
6.04M
                    RE(d2) = RE(c2) - IM(c3);
270
6.04M
                    IM(d3) = IM(c2) - RE(c3);
271
6.04M
                    RE(d3) = RE(c2) + IM(c3);
272
6.04M
                    IM(d2) = IM(c2) + RE(c3);
273
274
6.04M
#if 1
275
6.04M
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
276
6.04M
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
277
6.04M
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
278
6.04M
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
279
#else
280
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
281
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
282
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
283
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
284
#endif
285
6.04M
                }
286
108k
            }
287
76.2k
        } else {
288
9.21k
            for (k = 0; k < l1; k++)
289
6.14k
            {
290
497k
                for (i = 0; i < ido; i++)
291
491k
                {
292
491k
                    ac = i + (3*k+1)*ido;
293
491k
                    ah = i + k * ido;
294
295
491k
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
296
491k
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
297
491k
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
298
491k
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
299
300
491k
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
301
491k
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
302
303
491k
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
304
491k
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
305
306
491k
                    RE(d2) = RE(c2) + IM(c3);
307
491k
                    IM(d3) = IM(c2) + RE(c3);
308
491k
                    RE(d3) = RE(c2) - IM(c3);
309
491k
                    IM(d2) = IM(c2) - RE(c3);
310
311
491k
#if 1
312
491k
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
313
491k
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
314
491k
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
315
491k
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
316
#else
317
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
318
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
319
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
320
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
321
#endif
322
491k
                }
323
6.14k
            }
324
3.07k
        }
325
79.3k
    }
326
79.3k
}
cfft.c:passf3
Line
Count
Source
184
122k
{
185
122k
    static real_t taur = FRAC_CONST(-0.5);
186
122k
    static real_t taui = FRAC_CONST(0.866025403784439);
187
122k
    uint16_t i, k, ac, ah;
188
122k
    complex_t c2, c3, d2, d3, t2;
189
190
122k
    if (ido == 1)
191
0
    {
192
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
193
        // TODO: remove this code once fuzzer proves it is totally unreahable
194
        // 3 is never the the biggest factor for supported frame lengths;
195
        // consequently `ido` should never be 1.
196
0
        __builtin_trap();
197
        /*
198
#endif
199
        if (isign == 1)
200
        {
201
            for (k = 0; k < l1; k++)
202
            {
203
                ac = 3*k+1;
204
                ah = k;
205
206
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
207
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
208
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
209
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
210
211
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
212
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
213
214
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
215
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
216
217
                RE(ch[ah+l1]) = RE(c2) - IM(c3);
218
                IM(ch[ah+l1]) = IM(c2) + RE(c3);
219
                RE(ch[ah+2*l1]) = RE(c2) + IM(c3);
220
                IM(ch[ah+2*l1]) = IM(c2) - RE(c3);
221
            }
222
        } else {
223
            for (k = 0; k < l1; k++)
224
            {
225
                ac = 3*k+1;
226
                ah = k;
227
228
                RE(t2) = RE(cc[ac]) + RE(cc[ac+1]);
229
                IM(t2) = IM(cc[ac]) + IM(cc[ac+1]);
230
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),taur);
231
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),taur);
232
233
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2);
234
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2);
235
236
                RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+1])), taui);
237
                IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+1])), taui);
238
239
                RE(ch[ah+l1]) = RE(c2) + IM(c3);
240
                IM(ch[ah+l1]) = IM(c2) - RE(c3);
241
                RE(ch[ah+2*l1]) = RE(c2) - IM(c3);
242
                IM(ch[ah+2*l1]) = IM(c2) + RE(c3);
243
            }
244
        }
245
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
246
        */
247
0
#endif
248
122k
    } else {
249
122k
        if (isign == 1)
250
118k
        {
251
293k
            for (k = 0; k < l1; k++)
252
175k
            {
253
10.5M
                for (i = 0; i < ido; i++)
254
10.4M
                {
255
10.4M
                    ac = i + (3*k+1)*ido;
256
10.4M
                    ah = i + k * ido;
257
258
10.4M
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
259
10.4M
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
260
10.4M
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
261
10.4M
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
262
263
10.4M
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
264
10.4M
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
265
266
10.4M
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
267
10.4M
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
268
269
10.4M
                    RE(d2) = RE(c2) - IM(c3);
270
10.4M
                    IM(d3) = IM(c2) - RE(c3);
271
10.4M
                    RE(d3) = RE(c2) + IM(c3);
272
10.4M
                    IM(d2) = IM(c2) + RE(c3);
273
274
10.4M
#if 1
275
10.4M
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
276
10.4M
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
277
10.4M
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
278
10.4M
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
279
#else
280
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
281
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
282
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
283
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
284
#endif
285
10.4M
                }
286
175k
            }
287
118k
        } else {
288
11.5k
            for (k = 0; k < l1; k++)
289
7.66k
            {
290
620k
                for (i = 0; i < ido; i++)
291
613k
                {
292
613k
                    ac = i + (3*k+1)*ido;
293
613k
                    ah = i + k * ido;
294
295
613k
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
296
613k
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
297
613k
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
298
613k
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
299
300
613k
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
301
613k
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
302
303
613k
                    RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
304
613k
                    IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
305
306
613k
                    RE(d2) = RE(c2) + IM(c3);
307
613k
                    IM(d3) = IM(c2) + RE(c3);
308
613k
                    RE(d3) = RE(c2) - IM(c3);
309
613k
                    IM(d2) = IM(c2) - RE(c3);
310
311
613k
#if 1
312
613k
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
313
613k
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
314
613k
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
315
613k
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
316
#else
317
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
318
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
319
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
320
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
321
#endif
322
613k
                }
323
7.66k
            }
324
3.87k
        }
325
122k
    }
326
122k
}
327
328
329
static void passf4pos(const uint16_t ido, const uint16_t l1, const complex_t *cc,
330
                      complex_t *ch, const complex_t *wa1, const complex_t *wa2,
331
                      const complex_t *wa3)
332
3.43M
{
333
3.43M
    uint16_t i, k, ac, ah;
334
335
3.43M
    if (ido == 1)
336
899k
    {
337
65.2M
        for (k = 0; k < l1; k++)
338
64.3M
        {
339
64.3M
            complex_t t1, t2, t3, t4;
340
341
64.3M
            ac = 4*k;
342
64.3M
            ah = k;
343
344
64.3M
            RE(t2) = RE(cc[ac])   + RE(cc[ac+2]);
345
64.3M
            RE(t1) = RE(cc[ac])   - RE(cc[ac+2]);
346
64.3M
            IM(t2) = IM(cc[ac])   + IM(cc[ac+2]);
347
64.3M
            IM(t1) = IM(cc[ac])   - IM(cc[ac+2]);
348
64.3M
            RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
349
64.3M
            IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
350
64.3M
            IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
351
64.3M
            RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
352
353
64.3M
            RE(ch[ah])      = RE(t2) + RE(t3);
354
64.3M
            RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
355
356
64.3M
            IM(ch[ah])      = IM(t2) + IM(t3);
357
64.3M
            IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
358
359
64.3M
            RE(ch[ah+l1])   = RE(t1) + RE(t4);
360
64.3M
            RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
361
362
64.3M
            IM(ch[ah+l1])   = IM(t1) + IM(t4);
363
64.3M
            IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
364
64.3M
        }
365
2.53M
    } else {
366
26.5M
        for (k = 0; k < l1; k++)
367
24.0M
        {
368
24.0M
            ac = 4*k*ido;
369
24.0M
            ah = k*ido;
370
371
232M
            for (i = 0; i < ido; i++)
372
208M
            {
373
208M
                complex_t c2, c3, c4, t1, t2, t3, t4;
374
375
208M
                RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
376
208M
                RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
377
208M
                IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
378
208M
                IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
379
208M
                RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
380
208M
                IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
381
208M
                IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
382
208M
                RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
383
384
208M
                RE(c2) = RE(t1) + RE(t4);
385
208M
                RE(c4) = RE(t1) - RE(t4);
386
387
208M
                IM(c2) = IM(t1) + IM(t4);
388
208M
                IM(c4) = IM(t1) - IM(t4);
389
390
208M
                RE(ch[ah+i]) = RE(t2) + RE(t3);
391
208M
                RE(c3)       = RE(t2) - RE(t3);
392
393
208M
                IM(ch[ah+i]) = IM(t2) + IM(t3);
394
208M
                IM(c3)       = IM(t2) - IM(t3);
395
396
208M
#if 1
397
208M
                ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
398
208M
                    IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
399
208M
                ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
400
208M
                    IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
401
208M
                ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
402
208M
                    IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
403
#else
404
                ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
405
                    RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
406
                ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
407
                    RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
408
                ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
409
                    RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
410
#endif
411
208M
            }
412
24.0M
        }
413
2.53M
    }
414
3.43M
}
415
416
static void passf4neg(const uint16_t ido, const uint16_t l1, const complex_t *cc,
417
                      complex_t *ch, const complex_t *wa1, const complex_t *wa2,
418
                      const complex_t *wa3)
419
38.0k
{
420
38.0k
    uint16_t i, k, ac, ah;
421
422
38.0k
    if (ido == 1)
423
6.04k
    {
424
766k
        for (k = 0; k < l1; k++)
425
760k
        {
426
760k
            complex_t t1, t2, t3, t4;
427
428
760k
            ac = 4*k;
429
760k
            ah = k;
430
431
760k
            RE(t2) = RE(cc[ac])   + RE(cc[ac+2]);
432
760k
            RE(t1) = RE(cc[ac])   - RE(cc[ac+2]);
433
760k
            IM(t2) = IM(cc[ac])   + IM(cc[ac+2]);
434
760k
            IM(t1) = IM(cc[ac])   - IM(cc[ac+2]);
435
760k
            RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
436
760k
            IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
437
760k
            IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
438
760k
            RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
439
440
760k
            RE(ch[ah])      = RE(t2) + RE(t3);
441
760k
            RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
442
443
760k
            IM(ch[ah])      = IM(t2) + IM(t3);
444
760k
            IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
445
446
760k
            RE(ch[ah+l1])   = RE(t1) - RE(t4);
447
760k
            RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
448
449
760k
            IM(ch[ah+l1])   = IM(t1) - IM(t4);
450
760k
            IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
451
760k
        }
452
32.0k
    } else {
453
488k
        for (k = 0; k < l1; k++)
454
456k
        {
455
456k
            ac = 4*k*ido;
456
456k
            ah = k*ido;
457
458
4.39M
            for (i = 0; i < ido; i++)
459
3.93M
            {
460
3.93M
                complex_t c2, c3, c4, t1, t2, t3, t4;
461
462
3.93M
                RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
463
3.93M
                RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
464
3.93M
                IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
465
3.93M
                IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
466
3.93M
                RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
467
3.93M
                IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
468
3.93M
                IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
469
3.93M
                RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
470
471
3.93M
                RE(c2) = RE(t1) - RE(t4);
472
3.93M
                RE(c4) = RE(t1) + RE(t4);
473
474
3.93M
                IM(c2) = IM(t1) - IM(t4);
475
3.93M
                IM(c4) = IM(t1) + IM(t4);
476
477
3.93M
                RE(ch[ah+i]) = RE(t2) + RE(t3);
478
3.93M
                RE(c3)       = RE(t2) - RE(t3);
479
480
3.93M
                IM(ch[ah+i]) = IM(t2) + IM(t3);
481
3.93M
                IM(c3)       = IM(t2) - IM(t3);
482
483
3.93M
#if 1
484
3.93M
                ComplexMult(&RE(ch[ah+i+l1*ido]), &IM(ch[ah+i+l1*ido]),
485
3.93M
                    RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
486
3.93M
                ComplexMult(&RE(ch[ah+i+2*l1*ido]), &IM(ch[ah+i+2*l1*ido]),
487
3.93M
                    RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
488
3.93M
                ComplexMult(&RE(ch[ah+i+3*l1*ido]), &IM(ch[ah+i+3*l1*ido]),
489
3.93M
                    RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
490
#else
491
                ComplexMult(&IM(ch[ah+i+l1*ido]), &RE(ch[ah+i+l1*ido]),
492
                    IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
493
                ComplexMult(&IM(ch[ah+i+2*l1*ido]), &RE(ch[ah+i+2*l1*ido]),
494
                    IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
495
                ComplexMult(&IM(ch[ah+i+3*l1*ido]), &RE(ch[ah+i+3*l1*ido]),
496
                    IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
497
#endif
498
3.93M
            }
499
456k
        }
500
32.0k
    }
501
38.0k
}
502
503
static void passf5(const uint16_t ido, const uint16_t l1, const complex_t *cc,
504
                   complex_t *ch, const complex_t *wa1, const complex_t *wa2, const complex_t *wa3,
505
                   const complex_t *wa4, const int8_t isign)
506
201k
{
507
201k
    static real_t tr11 = FRAC_CONST(0.309016994374947);
508
201k
    static real_t ti11 = FRAC_CONST(0.951056516295154);
509
201k
    static real_t tr12 = FRAC_CONST(-0.809016994374947);
510
201k
    static real_t ti12 = FRAC_CONST(0.587785252292473);
511
201k
    uint16_t k, ac, ah;
512
201k
    complex_t c2, c3, c4, c5, t2, t3, t4, t5;
513
201k
    (void)wa1;
514
201k
    (void)wa2;
515
201k
    (void)wa3;
516
201k
    (void)wa4;
517
518
201k
    if (ido == 1)
519
201k
    {
520
201k
        if (isign == 1)
521
194k
        {
522
10.0M
            for (k = 0; k < l1; k++)
523
9.87M
            {
524
9.87M
                ac = 5*k + 1;
525
9.87M
                ah = k;
526
527
9.87M
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
528
9.87M
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
529
9.87M
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
530
9.87M
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
531
9.87M
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
532
9.87M
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
533
9.87M
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
534
9.87M
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
535
536
9.87M
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
537
9.87M
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
538
539
9.87M
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
540
9.87M
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
541
9.87M
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
542
9.87M
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
543
544
9.87M
                ComplexMult(&RE(c5), &RE(c4),
545
9.87M
                    ti11, ti12, RE(t5), RE(t4));
546
9.87M
                ComplexMult(&IM(c5), &IM(c4),
547
9.87M
                    ti11, ti12, IM(t5), IM(t4));
548
549
9.87M
                RE(ch[ah+l1]) = RE(c2) - IM(c5);
550
9.87M
                IM(ch[ah+l1]) = IM(c2) + RE(c5);
551
9.87M
                RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
552
9.87M
                IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
553
9.87M
                RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
554
9.87M
                IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
555
9.87M
                RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
556
9.87M
                IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
557
9.87M
            }
558
194k
        } else {
559
669k
            for (k = 0; k < l1; k++)
560
662k
            {
561
662k
                ac = 5*k + 1;
562
662k
                ah = k;
563
564
662k
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
565
662k
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
566
662k
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
567
662k
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
568
662k
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
569
662k
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
570
662k
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
571
662k
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
572
573
662k
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
574
662k
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
575
576
662k
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
577
662k
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
578
662k
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
579
662k
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
580
581
662k
                ComplexMult(&RE(c4), &RE(c5),
582
662k
                    ti12, ti11, RE(t5), RE(t4));
583
662k
                ComplexMult(&IM(c4), &IM(c5),
584
662k
                    ti12, ti11, IM(t5), IM(t4));
585
586
662k
                RE(ch[ah+l1]) = RE(c2) + IM(c5);
587
662k
                IM(ch[ah+l1]) = IM(c2) - RE(c5);
588
662k
                RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
589
662k
                IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
590
662k
                RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
591
662k
                IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
592
662k
                RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
593
662k
                IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
594
662k
            }
595
6.95k
        }
596
201k
    } else {
597
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
598
        // TODO: remove this code once fuzzer proves it is totally unreahable
599
        // 5 if the biggest factor and it never repeated for supported frame
600
        // lengths; consequently `ido` should always be 1.
601
0
        __builtin_trap();
602
        /*
603
#else
604
        uint16_t i;
605
        complex_t d3, d4, d5, d2;
606
#endif
607
        if (isign == 1)
608
        {
609
            for (k = 0; k < l1; k++)
610
            {
611
                for (i = 0; i < ido; i++)
612
                {
613
                    ac = i + (k*5 + 1) * ido;
614
                    ah = i + k * ido;
615
616
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
617
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
618
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
619
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
620
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
621
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
622
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
623
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
624
625
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
626
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
627
628
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
629
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
630
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
631
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
632
633
                    ComplexMult(&RE(c5), &RE(c4),
634
                        ti11, ti12, RE(t5), RE(t4));
635
                    ComplexMult(&IM(c5), &IM(c4),
636
                        ti11, ti12, IM(t5), IM(t4));
637
638
                    IM(d2) = IM(c2) + RE(c5);
639
                    IM(d3) = IM(c3) + RE(c4);
640
                    RE(d4) = RE(c3) + IM(c4);
641
                    RE(d5) = RE(c2) + IM(c5);
642
                    RE(d2) = RE(c2) - IM(c5);
643
                    IM(d5) = IM(c2) - RE(c5);
644
                    RE(d3) = RE(c3) - IM(c4);
645
                    IM(d4) = IM(c3) - RE(c4);
646
647
#if 1
648
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
649
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
650
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
651
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
652
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
653
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
654
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
655
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
656
#else
657
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
658
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
659
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
660
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
661
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
662
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
663
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
664
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
665
#endif
666
                }
667
            }
668
        } else {
669
            for (k = 0; k < l1; k++)
670
            {
671
                for (i = 0; i < ido; i++)
672
                {
673
                    ac = i + (k*5 + 1) * ido;
674
                    ah = i + k * ido;
675
676
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
677
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
678
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
679
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
680
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
681
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
682
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
683
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
684
685
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
686
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
687
688
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
689
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
690
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
691
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
692
693
                    ComplexMult(&RE(c4), &RE(c5),
694
                        ti12, ti11, RE(t5), RE(t4));
695
                    ComplexMult(&IM(c4), &IM(c5),
696
                        ti12, ti11, IM(t5), IM(t4));
697
698
                    IM(d2) = IM(c2) - RE(c5);
699
                    IM(d3) = IM(c3) - RE(c4);
700
                    RE(d4) = RE(c3) - IM(c4);
701
                    RE(d5) = RE(c2) - IM(c5);
702
                    RE(d2) = RE(c2) + IM(c5);
703
                    IM(d5) = IM(c2) + RE(c5);
704
                    RE(d3) = RE(c3) + IM(c4);
705
                    IM(d4) = IM(c3) + RE(c4);
706
707
#if 1
708
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
709
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
710
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
711
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
712
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
713
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
714
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
715
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
716
#else
717
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
718
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
719
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
720
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
721
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
722
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
723
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
724
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
725
#endif
726
                }
727
            }
728
        }
729
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
730
        */
731
0
#endif
732
0
    }
733
201k
}
cfft.c:passf5
Line
Count
Source
506
79.3k
{
507
79.3k
    static real_t tr11 = FRAC_CONST(0.309016994374947);
508
79.3k
    static real_t ti11 = FRAC_CONST(0.951056516295154);
509
79.3k
    static real_t tr12 = FRAC_CONST(-0.809016994374947);
510
79.3k
    static real_t ti12 = FRAC_CONST(0.587785252292473);
511
79.3k
    uint16_t k, ac, ah;
512
79.3k
    complex_t c2, c3, c4, c5, t2, t3, t4, t5;
513
79.3k
    (void)wa1;
514
79.3k
    (void)wa2;
515
79.3k
    (void)wa3;
516
79.3k
    (void)wa4;
517
518
79.3k
    if (ido == 1)
519
79.3k
    {
520
79.3k
        if (isign == 1)
521
76.2k
        {
522
3.70M
            for (k = 0; k < l1; k++)
523
3.62M
            {
524
3.62M
                ac = 5*k + 1;
525
3.62M
                ah = k;
526
527
3.62M
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
528
3.62M
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
529
3.62M
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
530
3.62M
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
531
3.62M
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
532
3.62M
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
533
3.62M
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
534
3.62M
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
535
536
3.62M
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
537
3.62M
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
538
539
3.62M
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
540
3.62M
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
541
3.62M
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
542
3.62M
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
543
544
3.62M
                ComplexMult(&RE(c5), &RE(c4),
545
3.62M
                    ti11, ti12, RE(t5), RE(t4));
546
3.62M
                ComplexMult(&IM(c5), &IM(c4),
547
3.62M
                    ti11, ti12, IM(t5), IM(t4));
548
549
3.62M
                RE(ch[ah+l1]) = RE(c2) - IM(c5);
550
3.62M
                IM(ch[ah+l1]) = IM(c2) + RE(c5);
551
3.62M
                RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
552
3.62M
                IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
553
3.62M
                RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
554
3.62M
                IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
555
3.62M
                RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
556
3.62M
                IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
557
3.62M
            }
558
76.2k
        } else {
559
297k
            for (k = 0; k < l1; k++)
560
294k
            {
561
294k
                ac = 5*k + 1;
562
294k
                ah = k;
563
564
294k
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
565
294k
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
566
294k
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
567
294k
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
568
294k
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
569
294k
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
570
294k
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
571
294k
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
572
573
294k
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
574
294k
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
575
576
294k
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
577
294k
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
578
294k
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
579
294k
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
580
581
294k
                ComplexMult(&RE(c4), &RE(c5),
582
294k
                    ti12, ti11, RE(t5), RE(t4));
583
294k
                ComplexMult(&IM(c4), &IM(c5),
584
294k
                    ti12, ti11, IM(t5), IM(t4));
585
586
294k
                RE(ch[ah+l1]) = RE(c2) + IM(c5);
587
294k
                IM(ch[ah+l1]) = IM(c2) - RE(c5);
588
294k
                RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
589
294k
                IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
590
294k
                RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
591
294k
                IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
592
294k
                RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
593
294k
                IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
594
294k
            }
595
3.07k
        }
596
79.3k
    } else {
597
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
598
        // TODO: remove this code once fuzzer proves it is totally unreahable
599
        // 5 if the biggest factor and it never repeated for supported frame
600
        // lengths; consequently `ido` should always be 1.
601
0
        __builtin_trap();
602
        /*
603
#else
604
        uint16_t i;
605
        complex_t d3, d4, d5, d2;
606
#endif
607
        if (isign == 1)
608
        {
609
            for (k = 0; k < l1; k++)
610
            {
611
                for (i = 0; i < ido; i++)
612
                {
613
                    ac = i + (k*5 + 1) * ido;
614
                    ah = i + k * ido;
615
616
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
617
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
618
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
619
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
620
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
621
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
622
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
623
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
624
625
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
626
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
627
628
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
629
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
630
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
631
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
632
633
                    ComplexMult(&RE(c5), &RE(c4),
634
                        ti11, ti12, RE(t5), RE(t4));
635
                    ComplexMult(&IM(c5), &IM(c4),
636
                        ti11, ti12, IM(t5), IM(t4));
637
638
                    IM(d2) = IM(c2) + RE(c5);
639
                    IM(d3) = IM(c3) + RE(c4);
640
                    RE(d4) = RE(c3) + IM(c4);
641
                    RE(d5) = RE(c2) + IM(c5);
642
                    RE(d2) = RE(c2) - IM(c5);
643
                    IM(d5) = IM(c2) - RE(c5);
644
                    RE(d3) = RE(c3) - IM(c4);
645
                    IM(d4) = IM(c3) - RE(c4);
646
647
#if 1
648
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
649
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
650
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
651
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
652
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
653
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
654
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
655
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
656
#else
657
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
658
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
659
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
660
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
661
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
662
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
663
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
664
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
665
#endif
666
                }
667
            }
668
        } else {
669
            for (k = 0; k < l1; k++)
670
            {
671
                for (i = 0; i < ido; i++)
672
                {
673
                    ac = i + (k*5 + 1) * ido;
674
                    ah = i + k * ido;
675
676
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
677
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
678
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
679
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
680
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
681
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
682
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
683
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
684
685
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
686
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
687
688
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
689
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
690
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
691
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
692
693
                    ComplexMult(&RE(c4), &RE(c5),
694
                        ti12, ti11, RE(t5), RE(t4));
695
                    ComplexMult(&IM(c4), &IM(c5),
696
                        ti12, ti11, IM(t5), IM(t4));
697
698
                    IM(d2) = IM(c2) - RE(c5);
699
                    IM(d3) = IM(c3) - RE(c4);
700
                    RE(d4) = RE(c3) - IM(c4);
701
                    RE(d5) = RE(c2) - IM(c5);
702
                    RE(d2) = RE(c2) + IM(c5);
703
                    IM(d5) = IM(c2) + RE(c5);
704
                    RE(d3) = RE(c3) + IM(c4);
705
                    IM(d4) = IM(c3) + RE(c4);
706
707
#if 1
708
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
709
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
710
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
711
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
712
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
713
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
714
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
715
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
716
#else
717
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
718
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
719
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
720
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
721
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
722
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
723
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
724
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
725
#endif
726
                }
727
            }
728
        }
729
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
730
        */
731
0
#endif
732
0
    }
733
79.3k
}
cfft.c:passf5
Line
Count
Source
506
122k
{
507
122k
    static real_t tr11 = FRAC_CONST(0.309016994374947);
508
122k
    static real_t ti11 = FRAC_CONST(0.951056516295154);
509
122k
    static real_t tr12 = FRAC_CONST(-0.809016994374947);
510
122k
    static real_t ti12 = FRAC_CONST(0.587785252292473);
511
122k
    uint16_t k, ac, ah;
512
122k
    complex_t c2, c3, c4, c5, t2, t3, t4, t5;
513
122k
    (void)wa1;
514
122k
    (void)wa2;
515
122k
    (void)wa3;
516
122k
    (void)wa4;
517
518
122k
    if (ido == 1)
519
122k
    {
520
122k
        if (isign == 1)
521
118k
        {
522
6.36M
            for (k = 0; k < l1; k++)
523
6.24M
            {
524
6.24M
                ac = 5*k + 1;
525
6.24M
                ah = k;
526
527
6.24M
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
528
6.24M
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
529
6.24M
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
530
6.24M
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
531
6.24M
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
532
6.24M
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
533
6.24M
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
534
6.24M
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
535
536
6.24M
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
537
6.24M
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
538
539
6.24M
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
540
6.24M
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
541
6.24M
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
542
6.24M
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
543
544
6.24M
                ComplexMult(&RE(c5), &RE(c4),
545
6.24M
                    ti11, ti12, RE(t5), RE(t4));
546
6.24M
                ComplexMult(&IM(c5), &IM(c4),
547
6.24M
                    ti11, ti12, IM(t5), IM(t4));
548
549
6.24M
                RE(ch[ah+l1]) = RE(c2) - IM(c5);
550
6.24M
                IM(ch[ah+l1]) = IM(c2) + RE(c5);
551
6.24M
                RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
552
6.24M
                IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
553
6.24M
                RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
554
6.24M
                IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
555
6.24M
                RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
556
6.24M
                IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
557
6.24M
            }
558
118k
        } else {
559
371k
            for (k = 0; k < l1; k++)
560
367k
            {
561
367k
                ac = 5*k + 1;
562
367k
                ah = k;
563
564
367k
                RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
565
367k
                IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
566
367k
                RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
567
367k
                IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
568
367k
                RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
569
367k
                IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
570
367k
                RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
571
367k
                IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
572
573
367k
                RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
574
367k
                IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
575
576
367k
                RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
577
367k
                IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
578
367k
                RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
579
367k
                IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
580
581
367k
                ComplexMult(&RE(c4), &RE(c5),
582
367k
                    ti12, ti11, RE(t5), RE(t4));
583
367k
                ComplexMult(&IM(c4), &IM(c5),
584
367k
                    ti12, ti11, IM(t5), IM(t4));
585
586
367k
                RE(ch[ah+l1]) = RE(c2) + IM(c5);
587
367k
                IM(ch[ah+l1]) = IM(c2) - RE(c5);
588
367k
                RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
589
367k
                IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
590
367k
                RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
591
367k
                IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
592
367k
                RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
593
367k
                IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
594
367k
            }
595
3.87k
        }
596
122k
    } else {
597
0
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
598
        // TODO: remove this code once fuzzer proves it is totally unreahable
599
        // 5 if the biggest factor and it never repeated for supported frame
600
        // lengths; consequently `ido` should always be 1.
601
0
        __builtin_trap();
602
        /*
603
#else
604
        uint16_t i;
605
        complex_t d3, d4, d5, d2;
606
#endif
607
        if (isign == 1)
608
        {
609
            for (k = 0; k < l1; k++)
610
            {
611
                for (i = 0; i < ido; i++)
612
                {
613
                    ac = i + (k*5 + 1) * ido;
614
                    ah = i + k * ido;
615
616
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
617
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
618
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
619
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
620
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
621
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
622
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
623
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
624
625
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
626
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
627
628
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
629
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
630
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
631
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
632
633
                    ComplexMult(&RE(c5), &RE(c4),
634
                        ti11, ti12, RE(t5), RE(t4));
635
                    ComplexMult(&IM(c5), &IM(c4),
636
                        ti11, ti12, IM(t5), IM(t4));
637
638
                    IM(d2) = IM(c2) + RE(c5);
639
                    IM(d3) = IM(c3) + RE(c4);
640
                    RE(d4) = RE(c3) + IM(c4);
641
                    RE(d5) = RE(c2) + IM(c5);
642
                    RE(d2) = RE(c2) - IM(c5);
643
                    IM(d5) = IM(c2) - RE(c5);
644
                    RE(d3) = RE(c3) - IM(c4);
645
                    IM(d4) = IM(c3) - RE(c4);
646
647
#if 1
648
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
649
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
650
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
651
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
652
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
653
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
654
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
655
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
656
#else
657
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
658
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
659
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
660
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
661
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
662
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
663
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
664
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
665
#endif
666
                }
667
            }
668
        } else {
669
            for (k = 0; k < l1; k++)
670
            {
671
                for (i = 0; i < ido; i++)
672
                {
673
                    ac = i + (k*5 + 1) * ido;
674
                    ah = i + k * ido;
675
676
                    RE(t2) = RE(cc[ac]) + RE(cc[ac+3*ido]);
677
                    IM(t2) = IM(cc[ac]) + IM(cc[ac+3*ido]);
678
                    RE(t3) = RE(cc[ac+ido]) + RE(cc[ac+2*ido]);
679
                    IM(t3) = IM(cc[ac+ido]) + IM(cc[ac+2*ido]);
680
                    RE(t4) = RE(cc[ac+ido]) - RE(cc[ac+2*ido]);
681
                    IM(t4) = IM(cc[ac+ido]) - IM(cc[ac+2*ido]);
682
                    RE(t5) = RE(cc[ac]) - RE(cc[ac+3*ido]);
683
                    IM(t5) = IM(cc[ac]) - IM(cc[ac+3*ido]);
684
685
                    RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2) + RE(t3);
686
                    IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2) + IM(t3);
687
688
                    RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
689
                    IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
690
                    RE(c3) = RE(cc[ac-ido]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
691
                    IM(c3) = IM(cc[ac-ido]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
692
693
                    ComplexMult(&RE(c4), &RE(c5),
694
                        ti12, ti11, RE(t5), RE(t4));
695
                    ComplexMult(&IM(c4), &IM(c5),
696
                        ti12, ti11, IM(t5), IM(t4));
697
698
                    IM(d2) = IM(c2) - RE(c5);
699
                    IM(d3) = IM(c3) - RE(c4);
700
                    RE(d4) = RE(c3) - IM(c4);
701
                    RE(d5) = RE(c2) - IM(c5);
702
                    RE(d2) = RE(c2) + IM(c5);
703
                    IM(d5) = IM(c2) + RE(c5);
704
                    RE(d3) = RE(c3) + IM(c4);
705
                    IM(d4) = IM(c3) + RE(c4);
706
707
#if 1
708
                    ComplexMult(&RE(ch[ah+l1*ido]), &IM(ch[ah+l1*ido]),
709
                        RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
710
                    ComplexMult(&RE(ch[ah+2*l1*ido]), &IM(ch[ah+2*l1*ido]),
711
                        RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
712
                    ComplexMult(&RE(ch[ah+3*l1*ido]), &IM(ch[ah+3*l1*ido]),
713
                        RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
714
                    ComplexMult(&RE(ch[ah+4*l1*ido]), &IM(ch[ah+4*l1*ido]),
715
                        RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
716
#else
717
                    ComplexMult(&IM(ch[ah+l1*ido]), &RE(ch[ah+l1*ido]),
718
                        IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
719
                    ComplexMult(&IM(ch[ah+2*l1*ido]), &RE(ch[ah+2*l1*ido]),
720
                        IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
721
                    ComplexMult(&IM(ch[ah+3*l1*ido]), &RE(ch[ah+3*l1*ido]),
722
                        IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
723
                    ComplexMult(&IM(ch[ah+4*l1*ido]), &RE(ch[ah+4*l1*ido]),
724
                        IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
725
#endif
726
                }
727
            }
728
        }
729
#ifdef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
730
        */
731
0
#endif
732
0
    }
733
122k
}
734
735
736
/*----------------------------------------------------------------------
737
   cfftf1, cfftf, cfftb, cffti1, cffti. Complex FFTs.
738
  ----------------------------------------------------------------------*/
739
740
static INLINE void cfftf1pos(uint16_t n, complex_t *c, complex_t *ch,
741
                             const uint16_t *ifac, const complex_t *wa,
742
                             const int8_t isign)
743
1.09M
{
744
1.09M
    uint16_t i;
745
1.09M
    uint16_t k1, l1, l2;
746
1.09M
    uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido;
747
748
1.09M
    nf = ifac[1];
749
1.09M
    na = 0;
750
1.09M
    l1 = 1;
751
1.09M
    iw = 0;
752
753
5.44M
    for (k1 = 2; k1 <= nf+1; k1++)
754
4.35M
    {
755
4.35M
        ip = ifac[k1];
756
4.35M
        l2 = ip*l1;
757
4.35M
        ido = n / l2;
758
759
4.35M
        switch (ip)
760
4.35M
        {
761
3.43M
        case 4:
762
3.43M
            ix2 = iw + ido;
763
3.43M
            ix3 = ix2 + ido;
764
765
3.43M
            if (na == 0)
766
1.89M
                passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
767
1.54M
            else
768
1.54M
                passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
769
770
3.43M
            na = 1 - na;
771
3.43M
            break;
772
533k
        case 2:
773
533k
            if (na == 0)
774
533k
                passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
775
0
            else
776
0
                passf2pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
777
778
533k
            na = 1 - na;
779
533k
            break;
780
194k
        case 3:
781
194k
            ix2 = iw + ido;
782
783
194k
            if (na == 0)
784
105k
                passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
785
88.7k
            else
786
88.7k
                passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
787
788
194k
            na = 1 - na;
789
194k
            break;
790
194k
        case 5:
791
194k
            ix2 = iw + ido;
792
194k
            ix3 = ix2 + ido;
793
194k
            ix4 = ix3 + ido;
794
795
194k
            if (na == 0)
796
192k
                passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
797
2.13k
            else
798
2.13k
                passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
799
800
194k
            na = 1 - na;
801
194k
            break;
802
4.35M
        }
803
804
4.35M
        l1 = l2;
805
4.35M
        iw += (ip-1) * ido;
806
4.35M
    }
807
808
1.09M
    if (na == 0)
809
3.97k
        return;
810
811
306M
    for (i = 0; i < n; i++)
812
305M
    {
813
305M
        RE(c[i]) = RE(ch[i]);
814
305M
        IM(c[i]) = IM(ch[i]);
815
305M
    }
816
1.09M
}
817
818
static INLINE void cfftf1neg(uint16_t n, complex_t *c, complex_t *ch,
819
                             const uint16_t *ifac, const complex_t *wa,
820
                             const int8_t isign)
821
12.9k
{
822
12.9k
    uint16_t i;
823
12.9k
    uint16_t k1, l1, l2;
824
12.9k
    uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido;
825
826
12.9k
    nf = ifac[1];
827
12.9k
    na = 0;
828
12.9k
    l1 = 1;
829
12.9k
    iw = 0;
830
831
77.6k
    for (k1 = 2; k1 <= nf+1; k1++)
832
64.6k
    {
833
64.6k
        ip = ifac[k1];
834
64.6k
        l2 = ip*l1;
835
64.6k
        ido = n / l2;
836
837
64.6k
        switch (ip)
838
64.6k
        {
839
38.0k
        case 4:
840
38.0k
            ix2 = iw + ido;
841
38.0k
            ix3 = ix2 + ido;
842
843
38.0k
            if (na == 0)
844
19.0k
                passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
845
19.0k
            else
846
19.0k
                passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
847
848
38.0k
            na = 1 - na;
849
38.0k
            break;
850
12.6k
        case 2:
851
12.6k
            if (na == 0)
852
12.6k
                passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw]);
853
0
            else
854
0
                passf2neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw]);
855
856
12.6k
            na = 1 - na;
857
12.6k
            break;
858
6.95k
        case 3:
859
6.95k
            ix2 = iw + ido;
860
861
6.95k
            if (na == 0)
862
101
                passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
863
6.85k
            else
864
6.85k
                passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
865
866
6.95k
            na = 1 - na;
867
6.95k
            break;
868
6.95k
        case 5:
869
6.95k
            ix2 = iw + ido;
870
6.95k
            ix3 = ix2 + ido;
871
6.95k
            ix4 = ix3 + ido;
872
873
6.95k
            if (na == 0)
874
6.85k
                passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
875
101
            else
876
101
                passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
877
878
6.95k
            na = 1 - na;
879
6.95k
            break;
880
64.6k
        }
881
882
64.6k
        l1 = l2;
883
64.6k
        iw += (ip-1) * ido;
884
64.6k
    }
885
886
12.9k
    if (na == 0)
887
301
        return;
888
889
6.29M
    for (i = 0; i < n; i++)
890
6.27M
    {
891
6.27M
        RE(c[i]) = RE(ch[i]);
892
6.27M
        IM(c[i]) = IM(ch[i]);
893
6.27M
    }
894
12.6k
}
895
896
void cfftf(cfft_info *cfft, complex_t *c)
897
12.9k
{
898
12.9k
    cfftf1neg(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, -1);
899
12.9k
}
900
901
void cfftb(cfft_info *cfft, complex_t *c)
902
1.09M
{
903
1.09M
    cfftf1pos(cfft->n, c, cfft->work, (const uint16_t*)cfft->ifac, (const complex_t*)cfft->tab, +1);
904
1.09M
}
905
906
static void cffti1(uint16_t n, complex_t *wa, uint16_t *ifac)
907
64.0k
{
908
64.0k
    static uint16_t ntryh[4] = {3, 4, 2, 5};
909
#ifndef FIXED_POINT
910
    real_t arg, argh, argld, fi;
911
    uint16_t ido, ipm;
912
    uint16_t i1, k1, l1, l2;
913
    uint16_t ld, ii, ip;
914
#endif
915
64.0k
    uint16_t ntry = 0, i, j;
916
64.0k
    uint16_t ib;
917
64.0k
    uint16_t nf, nl, nq, nr;
918
919
64.0k
    nl = n;
920
64.0k
    nf = 0;
921
64.0k
    j = 0;
922
923
174k
startloop:
924
174k
    j++;
925
926
174k
    if (j <= 4)
927
174k
        ntry = ntryh[j-1];
928
0
    else
929
0
        ntry += 2;
930
931
174k
    do
932
366k
    {
933
366k
        nq = nl / ntry;
934
366k
        nr = nl - ntry*nq;
935
936
366k
        if (nr != 0)
937
110k
            goto startloop;
938
939
256k
        nf++;
940
256k
        ifac[nf+1] = ntry;
941
256k
        nl = nq;
942
943
256k
        if (ntry == 2 && nf != 1)
944
24.7k
        {
945
118k
            for (i = 2; i <= nf; i++)
946
93.7k
            {
947
93.7k
                ib = nf - i + 2;
948
93.7k
                ifac[ib+1] = ifac[ib];
949
93.7k
            }
950
24.7k
            ifac[2] = 2;
951
24.7k
        }
952
256k
    } while (nl != 1);
953
954
64.0k
    ifac[0] = n;
955
64.0k
    ifac[1] = nf;
956
957
#ifndef FIXED_POINT
958
    argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
959
    i = 0;
960
    l1 = 1;
961
962
213k
    for (k1 = 1; k1 <= nf; k1++)
963
170k
    {
964
170k
        ip = ifac[k1+1];
965
170k
        ld = 0;
966
170k
        l2 = l1*ip;
967
170k
        ido = n / l2;
968
170k
        ipm = ip - 1;
969
970
649k
        for (j = 0; j < ipm; j++)
971
478k
        {
972
478k
            i1 = i;
973
478k
            RE(wa[i]) = 1.0;
974
478k
            IM(wa[i]) = 0.0;
975
            ld += l1;
976
            fi = 0;
977
            argld = ld*argh;
978
979
12.3M
            for (ii = 0; ii < ido; ii++)
980
11.8M
            {
981
11.8M
                i++;
982
11.8M
                fi++;
983
11.8M
                arg = fi * argld;
984
11.8M
                RE(wa[i]) = (real_t)cos(arg);
985
11.8M
#if 1
986
11.8M
                IM(wa[i]) = (real_t)sin(arg);
987
#else
988
                IM(wa[i]) = (real_t)-sin(arg);
989
#endif
990
11.8M
            }
991
992
478k
            if (ip > 5)
993
0
            {
994
0
                RE(wa[i1]) = RE(wa[i]);
995
0
                IM(wa[i1]) = IM(wa[i]);
996
0
            }
997
478k
        }
998
170k
        l1 = l2;
999
170k
    }
1000
#else  // FIXED_POINT
1001
    (void)wa;  /* whoa! does it work for FIXED_POINT? */
1002
#endif
1003
64.0k
}
cfft.c:cffti1
Line
Count
Source
907
21.2k
{
908
21.2k
    static uint16_t ntryh[4] = {3, 4, 2, 5};
909
#ifndef FIXED_POINT
910
    real_t arg, argh, argld, fi;
911
    uint16_t ido, ipm;
912
    uint16_t i1, k1, l1, l2;
913
    uint16_t ld, ii, ip;
914
#endif
915
21.2k
    uint16_t ntry = 0, i, j;
916
21.2k
    uint16_t ib;
917
21.2k
    uint16_t nf, nl, nq, nr;
918
919
21.2k
    nl = n;
920
21.2k
    nf = 0;
921
21.2k
    j = 0;
922
923
58.0k
startloop:
924
58.0k
    j++;
925
926
58.0k
    if (j <= 4)
927
58.0k
        ntry = ntryh[j-1];
928
0
    else
929
0
        ntry += 2;
930
931
58.0k
    do
932
121k
    {
933
121k
        nq = nl / ntry;
934
121k
        nr = nl - ntry*nq;
935
936
121k
        if (nr != 0)
937
36.7k
            goto startloop;
938
939
85.1k
        nf++;
940
85.1k
        ifac[nf+1] = ntry;
941
85.1k
        nl = nq;
942
943
85.1k
        if (ntry == 2 && nf != 1)
944
7.64k
        {
945
36.4k
            for (i = 2; i <= nf; i++)
946
28.8k
            {
947
28.8k
                ib = nf - i + 2;
948
28.8k
                ifac[ib+1] = ifac[ib];
949
28.8k
            }
950
7.64k
            ifac[2] = 2;
951
7.64k
        }
952
85.1k
    } while (nl != 1);
953
954
21.2k
    ifac[0] = n;
955
21.2k
    ifac[1] = nf;
956
957
#ifndef FIXED_POINT
958
    argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
959
    i = 0;
960
    l1 = 1;
961
962
    for (k1 = 1; k1 <= nf; k1++)
963
    {
964
        ip = ifac[k1+1];
965
        ld = 0;
966
        l2 = l1*ip;
967
        ido = n / l2;
968
        ipm = ip - 1;
969
970
        for (j = 0; j < ipm; j++)
971
        {
972
            i1 = i;
973
            RE(wa[i]) = 1.0;
974
            IM(wa[i]) = 0.0;
975
            ld += l1;
976
            fi = 0;
977
            argld = ld*argh;
978
979
            for (ii = 0; ii < ido; ii++)
980
            {
981
                i++;
982
                fi++;
983
                arg = fi * argld;
984
                RE(wa[i]) = (real_t)cos(arg);
985
#if 1
986
                IM(wa[i]) = (real_t)sin(arg);
987
#else
988
                IM(wa[i]) = (real_t)-sin(arg);
989
#endif
990
            }
991
992
            if (ip > 5)
993
            {
994
                RE(wa[i1]) = RE(wa[i]);
995
                IM(wa[i1]) = IM(wa[i]);
996
            }
997
        }
998
        l1 = l2;
999
    }
1000
#else  // FIXED_POINT
1001
21.2k
    (void)wa;  /* whoa! does it work for FIXED_POINT? */
1002
21.2k
#endif
1003
21.2k
}
cfft.c:cffti1
Line
Count
Source
907
42.7k
{
908
42.7k
    static uint16_t ntryh[4] = {3, 4, 2, 5};
909
42.7k
#ifndef FIXED_POINT
910
42.7k
    real_t arg, argh, argld, fi;
911
42.7k
    uint16_t ido, ipm;
912
42.7k
    uint16_t i1, k1, l1, l2;
913
42.7k
    uint16_t ld, ii, ip;
914
42.7k
#endif
915
42.7k
    uint16_t ntry = 0, i, j;
916
42.7k
    uint16_t ib;
917
42.7k
    uint16_t nf, nl, nq, nr;
918
919
42.7k
    nl = n;
920
42.7k
    nf = 0;
921
42.7k
    j = 0;
922
923
116k
startloop:
924
116k
    j++;
925
926
116k
    if (j <= 4)
927
116k
        ntry = ntryh[j-1];
928
0
    else
929
0
        ntry += 2;
930
931
116k
    do
932
244k
    {
933
244k
        nq = nl / ntry;
934
244k
        nr = nl - ntry*nq;
935
936
244k
        if (nr != 0)
937
74.1k
            goto startloop;
938
939
170k
        nf++;
940
170k
        ifac[nf+1] = ntry;
941
170k
        nl = nq;
942
943
170k
        if (ntry == 2 && nf != 1)
944
17.1k
        {
945
82.0k
            for (i = 2; i <= nf; i++)
946
64.8k
            {
947
64.8k
                ib = nf - i + 2;
948
64.8k
                ifac[ib+1] = ifac[ib];
949
64.8k
            }
950
17.1k
            ifac[2] = 2;
951
17.1k
        }
952
170k
    } while (nl != 1);
953
954
42.7k
    ifac[0] = n;
955
42.7k
    ifac[1] = nf;
956
957
42.7k
#ifndef FIXED_POINT
958
42.7k
    argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
959
42.7k
    i = 0;
960
42.7k
    l1 = 1;
961
962
213k
    for (k1 = 1; k1 <= nf; k1++)
963
170k
    {
964
170k
        ip = ifac[k1+1];
965
170k
        ld = 0;
966
170k
        l2 = l1*ip;
967
170k
        ido = n / l2;
968
170k
        ipm = ip - 1;
969
970
649k
        for (j = 0; j < ipm; j++)
971
478k
        {
972
478k
            i1 = i;
973
478k
            RE(wa[i]) = 1.0;
974
478k
            IM(wa[i]) = 0.0;
975
478k
            ld += l1;
976
478k
            fi = 0;
977
478k
            argld = ld*argh;
978
979
12.3M
            for (ii = 0; ii < ido; ii++)
980
11.8M
            {
981
11.8M
                i++;
982
11.8M
                fi++;
983
11.8M
                arg = fi * argld;
984
11.8M
                RE(wa[i]) = (real_t)cos(arg);
985
11.8M
#if 1
986
11.8M
                IM(wa[i]) = (real_t)sin(arg);
987
#else
988
                IM(wa[i]) = (real_t)-sin(arg);
989
#endif
990
11.8M
            }
991
992
478k
            if (ip > 5)
993
0
            {
994
0
                RE(wa[i1]) = RE(wa[i]);
995
0
                IM(wa[i1]) = IM(wa[i]);
996
0
            }
997
478k
        }
998
170k
        l1 = l2;
999
170k
    }
1000
#else  // FIXED_POINT
1001
    (void)wa;  /* whoa! does it work for FIXED_POINT? */
1002
#endif
1003
42.7k
}
1004
1005
cfft_info *cffti(uint16_t n)
1006
3.29k
{
1007
3.29k
    cfft_info *cfft = (cfft_info*)faad_malloc(sizeof(cfft_info));
1008
1009
3.29k
    cfft->n = n;
1010
3.29k
    cfft->work = (complex_t*)faad_malloc(n*sizeof(complex_t));
1011
1012
#ifndef FIXED_POINT
1013
    cfft->tab = (complex_t*)faad_malloc(n*sizeof(complex_t));
1014
1015
    cffti1(n, cfft->tab, cfft->ifac);
1016
#else
1017
3.29k
    cffti1(n, NULL, cfft->ifac);
1018
1019
3.29k
    switch (n)
1020
3.29k
    {
1021
1.12k
    case 64: cfft->tab = (complex_t*)cfft_tab_64; break;
1022
1.12k
    case 512: cfft->tab = (complex_t*)cfft_tab_512; break;
1023
#ifdef LD_DEC
1024
    case 256: cfft->tab = (complex_t*)cfft_tab_256; break;
1025
#endif
1026
1027
0
#ifdef ALLOW_SMALL_FRAMELENGTH
1028
523
    case 60: cfft->tab = (complex_t*)cfft_tab_60; break;
1029
523
    case 480: cfft->tab = (complex_t*)cfft_tab_480; break;
1030
#ifdef LD_DEC
1031
    case 240: cfft->tab = (complex_t*)cfft_tab_240; break;
1032
#endif
1033
0
#endif
1034
0
    case 128: cfft->tab = (complex_t*)cfft_tab_128; break;
1035
3.29k
    }
1036
3.29k
#endif
1037
1038
3.29k
    return cfft;
1039
3.29k
}
1040
1041
void cfftu(cfft_info *cfft)
1042
21.2k
{
1043
21.2k
    if (cfft->work) faad_free(cfft->work);
1044
#ifndef FIXED_POINT
1045
    if (cfft->tab) faad_free(cfft->tab);
1046
#endif
1047
1048
21.2k
    if (cfft) faad_free(cfft);
1049
21.2k
}
1050