/proc/self/cwd/libfaad/cfft.c

Line	Count	Source
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	160k	{
73	160k	uint16_t i, k, ah, ac;
74
75	160k	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	160k	} else {
98	321k	for (k = 0; k < l1; k++)
99	160k	{
100	160k	ah = k*ido;
101	160k	ac = 2kido;
102
103	40.8M	for (i = 0; i < ido; i++)
104	40.6M	{
105	40.6M	complex_t t2;
106
107	40.6M	RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
108	40.6M	RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
109
110	40.6M	IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
111	40.6M	IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
112
113	40.6M	#if 1
114	40.6M	ComplexMult(&IM(ch[ah+i+l1ido]), &RE(ch[ah+i+l1ido]),
115	40.6M	IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
116		#else
117		ComplexMult(&RE(ch[ah+i+l1ido]), &IM(ch[ah+i+l1ido]),
118		RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
119		#endif
120	40.6M	}
121	160k	}
122	160k	}
123	160k	}
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	6.33k	{
128	6.33k	uint16_t i, k, ah, ac;
129
130	6.33k	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	6.33k	} else {
153	12.6k	for (k = 0; k < l1; k++)
154	6.33k	{
155	6.33k	ah = k*ido;
156	6.33k	ac = 2kido;
157
158	1.57M	for (i = 0; i < ido; i++)
159	1.56M	{
160	1.56M	complex_t t2;
161
162	1.56M	RE(ch[ah+i]) = RE(cc[ac+i]) + RE(cc[ac+i+ido]);
163	1.56M	RE(t2) = RE(cc[ac+i]) - RE(cc[ac+i+ido]);
164
165	1.56M	IM(ch[ah+i]) = IM(cc[ac+i]) + IM(cc[ac+i+ido]);
166	1.56M	IM(t2) = IM(cc[ac+i]) - IM(cc[ac+i+ido]);
167
168	1.56M	#if 1
169	1.56M	ComplexMult(&RE(ch[ah+i+l1ido]), &IM(ch[ah+i+l1ido]),
170	1.56M	RE(t2), IM(t2), RE(wa[i]), IM(wa[i]));
171		#else
172		ComplexMult(&IM(ch[ah+i+l1ido]), &RE(ch[ah+i+l1ido]),
173		IM(t2), RE(t2), RE(wa[i]), IM(wa[i]));
174		#endif
175	1.56M	}
176	6.33k	}
177	6.33k	}
178	6.33k	}
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	62.4k	{
185	62.4k	static real_t taur = FRAC_CONST(-0.5);
186	62.4k	static real_t taui = FRAC_CONST(0.866025403784439);
187	62.4k	uint16_t i, k, ac, ah;
188	62.4k	complex_t c2, c3, d2, d3, t2;
189
190	62.4k	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	62.4k	} else {
249	62.4k	if (isign == 1)
250	58.8k	{
251	144k	for (k = 0; k < l1; k++)
252	85.5k	{
253	5.04M	for (i = 0; i < ido; i++)
254	4.96M	{
255	4.96M	ac = i + (3k+1)ido;
256	4.96M	ah = i + k * ido;
257
258	4.96M	RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
259	4.96M	RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
260	4.96M	IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
261	4.96M	IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
262
263	4.96M	RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
264	4.96M	IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
265
266	4.96M	RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
267	4.96M	IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
268
269	4.96M	RE(d2) = RE(c2) - IM(c3);
270	4.96M	IM(d3) = IM(c2) - RE(c3);
271	4.96M	RE(d3) = RE(c2) + IM(c3);
272	4.96M	IM(d2) = IM(c2) + RE(c3);
273
274	4.96M	#if 1
275	4.96M	ComplexMult(&IM(ch[ah+l1ido]), &RE(ch[ah+l1ido]),
276	4.96M	IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
277	4.96M	ComplexMult(&IM(ch[ah+2l1ido]), &RE(ch[ah+2l1ido]),
278	4.96M	IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
279		#else
280		ComplexMult(&RE(ch[ah+l1ido]), &IM(ch[ah+l1ido]),
281		RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
282		ComplexMult(&RE(ch[ah+2l1ido]), &IM(ch[ah+2l1ido]),
283		RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
284		#endif
285	4.96M	}
286	85.5k	}
287	58.8k	} else {
288	10.8k	for (k = 0; k < l1; k++)
289	7.21k	{
290	584k	for (i = 0; i < ido; i++)
291	577k	{
292	577k	ac = i + (3k+1)ido;
293	577k	ah = i + k * ido;
294
295	577k	RE(t2) = RE(cc[ac]) + RE(cc[ac+ido]);
296	577k	RE(c2) = RE(cc[ac-ido]) + MUL_F(RE(t2),taur);
297	577k	IM(t2) = IM(cc[ac]) + IM(cc[ac+ido]);
298	577k	IM(c2) = IM(cc[ac-ido]) + MUL_F(IM(t2),taur);
299
300	577k	RE(ch[ah]) = RE(cc[ac-ido]) + RE(t2);
301	577k	IM(ch[ah]) = IM(cc[ac-ido]) + IM(t2);
302
303	577k	RE(c3) = MUL_F((RE(cc[ac]) - RE(cc[ac+ido])), taui);
304	577k	IM(c3) = MUL_F((IM(cc[ac]) - IM(cc[ac+ido])), taui);
305
306	577k	RE(d2) = RE(c2) + IM(c3);
307	577k	IM(d3) = IM(c2) + RE(c3);
308	577k	RE(d3) = RE(c2) - IM(c3);
309	577k	IM(d2) = IM(c2) - RE(c3);
310
311	577k	#if 1
312	577k	ComplexMult(&RE(ch[ah+l1ido]), &IM(ch[ah+l1ido]),
313	577k	RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
314	577k	ComplexMult(&RE(ch[ah+2l1ido]), &IM(ch[ah+2l1ido]),
315	577k	RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
316		#else
317		ComplexMult(&IM(ch[ah+l1ido]), &RE(ch[ah+l1ido]),
318		IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
319		ComplexMult(&IM(ch[ah+2l1ido]), &RE(ch[ah+2l1ido]),
320		IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
321		#endif
322	577k	}
323	7.21k	}
324	3.63k	}
325	62.4k	}
326	62.4k	}
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	1.10M	{
333	1.10M	uint16_t i, k, ac, ah;
334
335	1.10M	if (ido == 1)
336	293k	{
337	20.0M	for (k = 0; k < l1; k++)
338	19.7M	{
339	19.7M	complex_t t1, t2, t3, t4;
340
341	19.7M	ac = 4*k;
342	19.7M	ah = k;
343
344	19.7M	RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
345	19.7M	RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
346	19.7M	IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
347	19.7M	IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
348	19.7M	RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
349	19.7M	IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
350	19.7M	IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
351	19.7M	RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
352
353	19.7M	RE(ch[ah]) = RE(t2) + RE(t3);
354	19.7M	RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
355
356	19.7M	IM(ch[ah]) = IM(t2) + IM(t3);
357	19.7M	IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
358
359	19.7M	RE(ch[ah+l1]) = RE(t1) + RE(t4);
360	19.7M	RE(ch[ah+3*l1]) = RE(t1) - RE(t4);
361
362	19.7M	IM(ch[ah+l1]) = IM(t1) + IM(t4);
363	19.7M	IM(ch[ah+3*l1]) = IM(t1) - IM(t4);
364	19.7M	}
365	808k	} else {
366	8.15M	for (k = 0; k < l1; k++)
367	7.34M	{
368	7.34M	ac = 4kido;
369	7.34M	ah = k*ido;
370
371	71.0M	for (i = 0; i < ido; i++)
372	63.6M	{
373	63.6M	complex_t c2, c3, c4, t1, t2, t3, t4;
374
375	63.6M	RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
376	63.6M	RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
377	63.6M	IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
378	63.6M	IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
379	63.6M	RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
380	63.6M	IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
381	63.6M	IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
382	63.6M	RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
383
384	63.6M	RE(c2) = RE(t1) + RE(t4);
385	63.6M	RE(c4) = RE(t1) - RE(t4);
386
387	63.6M	IM(c2) = IM(t1) + IM(t4);
388	63.6M	IM(c4) = IM(t1) - IM(t4);
389
390	63.6M	RE(ch[ah+i]) = RE(t2) + RE(t3);
391	63.6M	RE(c3) = RE(t2) - RE(t3);
392
393	63.6M	IM(ch[ah+i]) = IM(t2) + IM(t3);
394	63.6M	IM(c3) = IM(t2) - IM(t3);
395
396	63.6M	#if 1
397	63.6M	ComplexMult(&IM(ch[ah+i+l1ido]), &RE(ch[ah+i+l1ido]),
398	63.6M	IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
399	63.6M	ComplexMult(&IM(ch[ah+i+2l1ido]), &RE(ch[ah+i+2l1ido]),
400	63.6M	IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
401	63.6M	ComplexMult(&IM(ch[ah+i+3l1ido]), &RE(ch[ah+i+3l1ido]),
402	63.6M	IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
403		#else
404		ComplexMult(&RE(ch[ah+i+l1ido]), &IM(ch[ah+i+l1ido]),
405		RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
406		ComplexMult(&RE(ch[ah+i+2l1ido]), &IM(ch[ah+i+2l1ido]),
407		RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
408		ComplexMult(&RE(ch[ah+i+3l1ido]), &IM(ch[ah+i+3l1ido]),
409		RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
410		#endif
411	63.6M	}
412	7.34M	}
413	808k	}
414	1.10M	}
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	18.9k	{
420	18.9k	uint16_t i, k, ac, ah;
421
422	18.9k	if (ido == 1)
423	2.93k	{
424	367k	for (k = 0; k < l1; k++)
425	364k	{
426	364k	complex_t t1, t2, t3, t4;
427
428	364k	ac = 4*k;
429	364k	ah = k;
430
431	364k	RE(t2) = RE(cc[ac]) + RE(cc[ac+2]);
432	364k	RE(t1) = RE(cc[ac]) - RE(cc[ac+2]);
433	364k	IM(t2) = IM(cc[ac]) + IM(cc[ac+2]);
434	364k	IM(t1) = IM(cc[ac]) - IM(cc[ac+2]);
435	364k	RE(t3) = RE(cc[ac+1]) + RE(cc[ac+3]);
436	364k	IM(t4) = RE(cc[ac+1]) - RE(cc[ac+3]);
437	364k	IM(t3) = IM(cc[ac+3]) + IM(cc[ac+1]);
438	364k	RE(t4) = IM(cc[ac+3]) - IM(cc[ac+1]);
439
440	364k	RE(ch[ah]) = RE(t2) + RE(t3);
441	364k	RE(ch[ah+2*l1]) = RE(t2) - RE(t3);
442
443	364k	IM(ch[ah]) = IM(t2) + IM(t3);
444	364k	IM(ch[ah+2*l1]) = IM(t2) - IM(t3);
445
446	364k	RE(ch[ah+l1]) = RE(t1) - RE(t4);
447	364k	RE(ch[ah+3*l1]) = RE(t1) + RE(t4);
448
449	364k	IM(ch[ah+l1]) = IM(t1) - IM(t4);
450	364k	IM(ch[ah+3*l1]) = IM(t1) + IM(t4);
451	364k	}
452	16.0k	} else {
453	243k	for (k = 0; k < l1; k++)
454	227k	{
455	227k	ac = 4kido;
456	227k	ah = k*ido;
457
458	2.18M	for (i = 0; i < ido; i++)
459	1.95M	{
460	1.95M	complex_t c2, c3, c4, t1, t2, t3, t4;
461
462	1.95M	RE(t2) = RE(cc[ac+i]) + RE(cc[ac+i+2*ido]);
463	1.95M	RE(t1) = RE(cc[ac+i]) - RE(cc[ac+i+2*ido]);
464	1.95M	IM(t2) = IM(cc[ac+i]) + IM(cc[ac+i+2*ido]);
465	1.95M	IM(t1) = IM(cc[ac+i]) - IM(cc[ac+i+2*ido]);
466	1.95M	RE(t3) = RE(cc[ac+i+ido]) + RE(cc[ac+i+3*ido]);
467	1.95M	IM(t4) = RE(cc[ac+i+ido]) - RE(cc[ac+i+3*ido]);
468	1.95M	IM(t3) = IM(cc[ac+i+3*ido]) + IM(cc[ac+i+ido]);
469	1.95M	RE(t4) = IM(cc[ac+i+3*ido]) - IM(cc[ac+i+ido]);
470
471	1.95M	RE(c2) = RE(t1) - RE(t4);
472	1.95M	RE(c4) = RE(t1) + RE(t4);
473
474	1.95M	IM(c2) = IM(t1) - IM(t4);
475	1.95M	IM(c4) = IM(t1) + IM(t4);
476
477	1.95M	RE(ch[ah+i]) = RE(t2) + RE(t3);
478	1.95M	RE(c3) = RE(t2) - RE(t3);
479
480	1.95M	IM(ch[ah+i]) = IM(t2) + IM(t3);
481	1.95M	IM(c3) = IM(t2) - IM(t3);
482
483	1.95M	#if 1
484	1.95M	ComplexMult(&RE(ch[ah+i+l1ido]), &IM(ch[ah+i+l1ido]),
485	1.95M	RE(c2), IM(c2), RE(wa1[i]), IM(wa1[i]));
486	1.95M	ComplexMult(&RE(ch[ah+i+2l1ido]), &IM(ch[ah+i+2l1ido]),
487	1.95M	RE(c3), IM(c3), RE(wa2[i]), IM(wa2[i]));
488	1.95M	ComplexMult(&RE(ch[ah+i+3l1ido]), &IM(ch[ah+i+3l1ido]),
489	1.95M	RE(c4), IM(c4), RE(wa3[i]), IM(wa3[i]));
490		#else
491		ComplexMult(&IM(ch[ah+i+l1ido]), &RE(ch[ah+i+l1ido]),
492		IM(c2), RE(c2), RE(wa1[i]), IM(wa1[i]));
493		ComplexMult(&IM(ch[ah+i+2l1ido]), &RE(ch[ah+i+2l1ido]),
494		IM(c3), RE(c3), RE(wa2[i]), IM(wa2[i]));
495		ComplexMult(&IM(ch[ah+i+3l1ido]), &RE(ch[ah+i+3l1ido]),
496		IM(c4), RE(c4), RE(wa3[i]), IM(wa3[i]));
497		#endif
498	1.95M	}
499	227k	}
500	16.0k	}
501	18.9k	}
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	62.4k	{
507	62.4k	static real_t tr11 = FRAC_CONST(0.309016994374947);
508	62.4k	static real_t ti11 = FRAC_CONST(0.951056516295154);
509	62.4k	static real_t tr12 = FRAC_CONST(-0.809016994374947);
510	62.4k	static real_t ti12 = FRAC_CONST(0.587785252292473);
511	62.4k	uint16_t k, ac, ah;
512	62.4k	complex_t c2, c3, c4, c5, t2, t3, t4, t5;
513	62.4k	(void)wa1;
514	62.4k	(void)wa2;
515	62.4k	(void)wa3;
516	62.4k	(void)wa4;
517
518	62.4k	if (ido == 1)
519	62.4k	{
520	62.4k	if (isign == 1)
521	58.8k	{
522	3.03M	for (k = 0; k < l1; k++)
523	2.97M	{
524	2.97M	ac = 5*k + 1;
525	2.97M	ah = k;
526
527	2.97M	RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
528	2.97M	IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
529	2.97M	RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
530	2.97M	IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
531	2.97M	RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
532	2.97M	IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
533	2.97M	RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
534	2.97M	IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
535
536	2.97M	RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
537	2.97M	IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
538
539	2.97M	RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
540	2.97M	IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
541	2.97M	RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
542	2.97M	IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
543
544	2.97M	ComplexMult(&RE(c5), &RE(c4),
545	2.97M	ti11, ti12, RE(t5), RE(t4));
546	2.97M	ComplexMult(&IM(c5), &IM(c4),
547	2.97M	ti11, ti12, IM(t5), IM(t4));
548
549	2.97M	RE(ch[ah+l1]) = RE(c2) - IM(c5);
550	2.97M	IM(ch[ah+l1]) = IM(c2) + RE(c5);
551	2.97M	RE(ch[ah+2*l1]) = RE(c3) - IM(c4);
552	2.97M	IM(ch[ah+2*l1]) = IM(c3) + RE(c4);
553	2.97M	RE(ch[ah+3*l1]) = RE(c3) + IM(c4);
554	2.97M	IM(ch[ah+3*l1]) = IM(c3) - RE(c4);
555	2.97M	RE(ch[ah+4*l1]) = RE(c2) + IM(c5);
556	2.97M	IM(ch[ah+4*l1]) = IM(c2) - RE(c5);
557	2.97M	}
558	58.8k	} else {
559	349k	for (k = 0; k < l1; k++)
560	346k	{
561	346k	ac = 5*k + 1;
562	346k	ah = k;
563
564	346k	RE(t2) = RE(cc[ac]) + RE(cc[ac+3]);
565	346k	IM(t2) = IM(cc[ac]) + IM(cc[ac+3]);
566	346k	RE(t3) = RE(cc[ac+1]) + RE(cc[ac+2]);
567	346k	IM(t3) = IM(cc[ac+1]) + IM(cc[ac+2]);
568	346k	RE(t4) = RE(cc[ac+1]) - RE(cc[ac+2]);
569	346k	IM(t4) = IM(cc[ac+1]) - IM(cc[ac+2]);
570	346k	RE(t5) = RE(cc[ac]) - RE(cc[ac+3]);
571	346k	IM(t5) = IM(cc[ac]) - IM(cc[ac+3]);
572
573	346k	RE(ch[ah]) = RE(cc[ac-1]) + RE(t2) + RE(t3);
574	346k	IM(ch[ah]) = IM(cc[ac-1]) + IM(t2) + IM(t3);
575
576	346k	RE(c2) = RE(cc[ac-1]) + MUL_F(RE(t2),tr11) + MUL_F(RE(t3),tr12);
577	346k	IM(c2) = IM(cc[ac-1]) + MUL_F(IM(t2),tr11) + MUL_F(IM(t3),tr12);
578	346k	RE(c3) = RE(cc[ac-1]) + MUL_F(RE(t2),tr12) + MUL_F(RE(t3),tr11);
579	346k	IM(c3) = IM(cc[ac-1]) + MUL_F(IM(t2),tr12) + MUL_F(IM(t3),tr11);
580
581	346k	ComplexMult(&RE(c4), &RE(c5),
582	346k	ti12, ti11, RE(t5), RE(t4));
583	346k	ComplexMult(&IM(c4), &IM(c5),
584	346k	ti12, ti11, IM(t5), IM(t4));
585
586	346k	RE(ch[ah+l1]) = RE(c2) + IM(c5);
587	346k	IM(ch[ah+l1]) = IM(c2) - RE(c5);
588	346k	RE(ch[ah+2*l1]) = RE(c3) + IM(c4);
589	346k	IM(ch[ah+2*l1]) = IM(c3) - RE(c4);
590	346k	RE(ch[ah+3*l1]) = RE(c3) - IM(c4);
591	346k	IM(ch[ah+3*l1]) = IM(c3) + RE(c4);
592	346k	RE(ch[ah+4*l1]) = RE(c2) - IM(c5);
593	346k	IM(ch[ah+4*l1]) = IM(c2) + RE(c5);
594	346k	}
595	3.63k	}
596	62.4k	} 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 + (k5 + 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+l1ido]), &RE(ch[ah+l1ido]),
649		IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
650		ComplexMult(&IM(ch[ah+2l1ido]), &RE(ch[ah+2l1ido]),
651		IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
652		ComplexMult(&IM(ch[ah+3l1ido]), &RE(ch[ah+3l1ido]),
653		IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
654		ComplexMult(&IM(ch[ah+4l1ido]), &RE(ch[ah+4l1ido]),
655		IM(d5), RE(d5), RE(wa4[i]), IM(wa4[i]));
656		#else
657		ComplexMult(&RE(ch[ah+l1ido]), &IM(ch[ah+l1ido]),
658		RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
659		ComplexMult(&RE(ch[ah+2l1ido]), &IM(ch[ah+2l1ido]),
660		RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
661		ComplexMult(&RE(ch[ah+3l1ido]), &IM(ch[ah+3l1ido]),
662		RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
663		ComplexMult(&RE(ch[ah+4l1ido]), &IM(ch[ah+4l1ido]),
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 + (k5 + 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+l1ido]), &IM(ch[ah+l1ido]),
709		RE(d2), IM(d2), RE(wa1[i]), IM(wa1[i]));
710		ComplexMult(&RE(ch[ah+2l1ido]), &IM(ch[ah+2l1ido]),
711		RE(d3), IM(d3), RE(wa2[i]), IM(wa2[i]));
712		ComplexMult(&RE(ch[ah+3l1ido]), &IM(ch[ah+3l1ido]),
713		RE(d4), IM(d4), RE(wa3[i]), IM(wa3[i]));
714		ComplexMult(&RE(ch[ah+4l1ido]), &IM(ch[ah+4l1ido]),
715		RE(d5), IM(d5), RE(wa4[i]), IM(wa4[i]));
716		#else
717		ComplexMult(&IM(ch[ah+l1ido]), &RE(ch[ah+l1ido]),
718		IM(d2), RE(d2), RE(wa1[i]), IM(wa1[i]));
719		ComplexMult(&IM(ch[ah+2l1ido]), &RE(ch[ah+2l1ido]),
720		IM(d3), RE(d3), RE(wa2[i]), IM(wa2[i]));
721		ComplexMult(&IM(ch[ah+3l1ido]), &RE(ch[ah+3l1ido]),
722		IM(d4), RE(d4), RE(wa3[i]), IM(wa3[i]));
723		ComplexMult(&IM(ch[ah+4l1ido]), &RE(ch[ah+4l1ido]),
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	62.4k	}
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	352k	{
744	352k	uint16_t i;
745	352k	uint16_t k1, l1, l2;
746	352k	uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido;
747
748	352k	nf = ifac[1];
749	352k	na = 0;
750	352k	l1 = 1;
751	352k	iw = 0;
752
753	1.73M	for (k1 = 2; k1 <= nf+1; k1++)
754	1.38M	{
755	1.38M	ip = ifac[k1];
756	1.38M	l2 = ip*l1;
757	1.38M	ido = n / l2;
758
759	1.38M	switch (ip)
760	1.38M	{
761	1.10M	case 4:
762	1.10M	ix2 = iw + ido;
763	1.10M	ix3 = ix2 + ido;
764
765	1.10M	if (na == 0)
766	614k	passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
767	487k	else
768	487k	passf4pos((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
769
770	1.10M	na = 1 - na;
771	1.10M	break;
772	160k	case 2:
773	160k	if (na == 0)
774	160k	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	160k	na = 1 - na;
779	160k	break;
780	58.8k	case 3:
781	58.8k	ix2 = iw + ido;
782
783	58.8k	if (na == 0)
784	32.1k	passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
785	26.7k	else
786	26.7k	passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
787
788	58.8k	na = 1 - na;
789	58.8k	break;
790	58.8k	case 5:
791	58.8k	ix2 = iw + ido;
792	58.8k	ix3 = ix2 + ido;
793	58.8k	ix4 = ix3 + ido;
794
795	58.8k	if (na == 0)
796	58.0k	passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
797	751	else
798	751	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	58.8k	na = 1 - na;
801	58.8k	break;
802	1.38M	}
803
804	1.38M	l1 = l2;
805	1.38M	iw += (ip-1) * ido;
806	1.38M	}
807
808	352k	if (na == 0)
809	1.90k	return;
810
811	93.7M	for (i = 0; i < n; i++)
812	93.4M	{
813	93.4M	RE(c[i]) = RE(ch[i]);
814	93.4M	IM(c[i]) = IM(ch[i]);
815	93.4M	}
816	350k	}
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	6.56k	{
822	6.56k	uint16_t i;
823	6.56k	uint16_t k1, l1, l2;
824	6.56k	uint16_t na, nf, ip, iw, ix2, ix3, ix4, ido;
825
826	6.56k	nf = ifac[1];
827	6.56k	na = 0;
828	6.56k	l1 = 1;
829	6.56k	iw = 0;
830
831	39.1k	for (k1 = 2; k1 <= nf+1; k1++)
832	32.6k	{
833	32.6k	ip = ifac[k1];
834	32.6k	l2 = ip*l1;
835	32.6k	ido = n / l2;
836
837	32.6k	switch (ip)
838	32.6k	{
839	18.9k	case 4:
840	18.9k	ix2 = iw + ido;
841	18.9k	ix3 = ix2 + ido;
842
843	18.9k	if (na == 0)
844	9.49k	passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3]);
845	9.49k	else
846	9.49k	passf4neg((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], &wa[ix3]);
847
848	18.9k	na = 1 - na;
849	18.9k	break;
850	6.33k	case 2:
851	6.33k	if (na == 0)
852	6.33k	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	6.33k	na = 1 - na;
857	6.33k	break;
858	3.63k	case 3:
859	3.63k	ix2 = iw + ido;
860
861	3.63k	if (na == 0)
862	59	passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], isign);
863	3.57k	else
864	3.57k	passf3((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)ch, c, &wa[iw], &wa[ix2], isign);
865
866	3.63k	na = 1 - na;
867	3.63k	break;
868	3.63k	case 5:
869	3.63k	ix2 = iw + ido;
870	3.63k	ix3 = ix2 + ido;
871	3.63k	ix4 = ix3 + ido;
872
873	3.63k	if (na == 0)
874	3.57k	passf5((const uint16_t)ido, (const uint16_t)l1, (const complex_t*)c, ch, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign);
875	59	else
876	59	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	3.63k	na = 1 - na;
879	3.63k	break;
880	32.6k	}
881
882	32.6k	l1 = l2;
883	32.6k	iw += (ip-1) * ido;
884	32.6k	}
885
886	6.56k	if (na == 0)
887	230	return;
888
889	3.13M	for (i = 0; i < n; i++)
890	3.13M	{
891	3.13M	RE(c[i]) = RE(ch[i]);
892	3.13M	IM(c[i]) = IM(ch[i]);
893	3.13M	}
894	6.33k	}
895
896		void cfftf(cfft_info cfft, complex_t c)
897	6.56k	{
898	6.56k	cfftf1neg(cfft->n, c, cfft->work, (const uint16_t)cfft->ifac, (const complex_t)cfft->tab, -1);
899	6.56k	}
900
901		void cfftb(cfft_info cfft, complex_t c)
902	352k	{
903	352k	cfftf1pos(cfft->n, c, cfft->work, (const uint16_t)cfft->ifac, (const complex_t)cfft->tab, +1);
904	352k	}
905
906		static void cffti1(uint16_t n, complex_t wa, uint16_t ifac)
907	20.0k	{
908	20.0k	static uint16_t ntryh[4] = {3, 4, 2, 5};
909	20.0k	#ifndef FIXED_POINT
910	20.0k	real_t arg, argh, argld, fi;
911	20.0k	uint16_t ido, ipm;
912	20.0k	uint16_t i1, k1, l1, l2;
913	20.0k	uint16_t ld, ii, ip;
914	20.0k	#endif
915	20.0k	uint16_t ntry = 0, i, j;
916	20.0k	uint16_t ib;
917	20.0k	uint16_t nf, nl, nq, nr;
918
919	20.0k	nl = n;
920	20.0k	nf = 0;
921	20.0k	j = 0;
922
923	52.8k	startloop:
924	52.8k	j++;
925
926	52.8k	if (j <= 4)
927	52.8k	ntry = ntryh[j-1];
928	0	else
929	0	ntry += 2;
930
931	52.8k	do
932	113k	{
933	113k	nq = nl / ntry;
934	113k	nr = nl - ntry*nq;
935
936	113k	if (nr != 0)
937	32.8k	goto startloop;
938
939	80.2k	nf++;
940	80.2k	ifac[nf+1] = ntry;
941	80.2k	nl = nq;
942
943	80.2k	if (ntry == 2 && nf != 1)
944	6.68k	{
945	32.2k	for (i = 2; i <= nf; i++)
946	25.5k	{
947	25.5k	ib = nf - i + 2;
948	25.5k	ifac[ib+1] = ifac[ib];
949	25.5k	}
950	6.68k	ifac[2] = 2;
951	6.68k	}
952	80.2k	} while (nl != 1);
953
954	20.0k	ifac[0] = n;
955	20.0k	ifac[1] = nf;
956
957	20.0k	#ifndef FIXED_POINT
958	20.0k	argh = (real_t)2.0*(real_t)M_PI / (real_t)n;
959	20.0k	i = 0;
960	20.0k	l1 = 1;
961
962	100k	for (k1 = 1; k1 <= nf; k1++)
963	80.2k	{
964	80.2k	ip = ifac[k1+1];
965	80.2k	ld = 0;
966	80.2k	l2 = l1*ip;
967	80.2k	ido = n / l2;
968	80.2k	ipm = ip - 1;
969
970	307k	for (j = 0; j < ipm; j++)
971	227k	{
972	227k	i1 = i;
973	227k	RE(wa[i]) = 1.0;
974	227k	IM(wa[i]) = 0.0;
975	227k	ld += l1;
976	227k	fi = 0;
977	227k	argld = ld*argh;
978
979	5.70M	for (ii = 0; ii < ido; ii++)
980	5.47M	{
981	5.47M	i++;
982	5.47M	fi++;
983	5.47M	arg = fi * argld;
984	5.47M	RE(wa[i]) = (real_t)cos(arg);
985	5.47M	#if 1
986	5.47M	IM(wa[i]) = (real_t)sin(arg);
987		#else
988		IM(wa[i]) = (real_t)-sin(arg);
989		#endif
990	5.47M	}
991
992	227k	if (ip > 5)
993	0	{
994	0	RE(wa[i1]) = RE(wa[i]);
995	0	IM(wa[i1]) = IM(wa[i]);
996	0	}
997	227k	}
998	80.2k	l1 = l2;
999	80.2k	}
1000		#else // FIXED_POINT
1001		(void)wa; /* whoa! does it work for FIXED_POINT? */
1002		#endif
1003	20.0k	}
1004
1005		cfft_info *cffti(uint16_t n)
1006	20.0k	{
1007	20.0k	cfft_info cfft = (cfft_info)faad_malloc(sizeof(cfft_info));
1008
1009	20.0k	cfft->n = n;
1010	20.0k	cfft->work = (complex_t)faad_malloc(nsizeof(complex_t));
1011
1012	20.0k	#ifndef FIXED_POINT
1013	20.0k	cfft->tab = (complex_t)faad_malloc(nsizeof(complex_t));
1014
1015	20.0k	cffti1(n, cfft->tab, cfft->ifac);
1016		#else
1017		cffti1(n, NULL, cfft->ifac);
1018
1019		switch (n)
1020		{
1021		case 64: cfft->tab = (complex_t*)cfft_tab_64; break;
1022		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		#ifdef ALLOW_SMALL_FRAMELENGTH
1028		case 60: cfft->tab = (complex_t*)cfft_tab_60; break;
1029		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		#endif
1034		case 128: cfft->tab = (complex_t*)cfft_tab_128; break;
1035		}
1036		#endif
1037
1038	20.0k	return cfft;
1039	20.0k	}
1040
1041		void cfftu(cfft_info *cfft)
1042	20.0k	{
1043	20.0k	if (cfft->work) faad_free(cfft->work);
1044	20.0k	#ifndef FIXED_POINT
1045	20.0k	if (cfft->tab) faad_free(cfft->tab);
1046	20.0k	#endif
1047
1048	20.0k	if (cfft) faad_free(cfft);
1049	20.0k	}
1050

Coverage Report

Created: 2025-11-16 06:18