/proc/self/cwd/external/fft2d/fftsg.c

Line	Count	Source (jump to first uncovered line)
1		/*
2		Fast Fourier/Cosine/Sine Transform
3		dimension :one
4		data length :power of 2
5		decimation :frequency
6		radix :split-radix
7		data :inplace
8		table :use
9		functions
10		cdft: Complex Discrete Fourier Transform
11		rdft: Real Discrete Fourier Transform
12		ddct: Discrete Cosine Transform
13		ddst: Discrete Sine Transform
14		dfct: Cosine Transform of RDFT (Real Symmetric DFT)
15		dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
16		function prototypes
17		void cdft(int, int, double , int , double *);
18		void rdft(int, int, double , int , double *);
19		void ddct(int, int, double , int , double *);
20		void ddst(int, int, double , int , double *);
21		void dfct(int, double , double , int , double );
22		void dfst(int, double , double , int , double );
23		macro definitions
24		USE_CDFT_PTHREADS : default=not defined
25		CDFT_THREADS_BEGIN_N : must be >= 512, default=8192
26		CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
27		USE_CDFT_WINTHREADS : default=not defined
28		CDFT_THREADS_BEGIN_N : must be >= 512, default=32768
29		CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288
30
31
32		-------- Complex DFT (Discrete Fourier Transform) --------
33		[definition]
34		<case1>
35		X[k] = sum_j=0^n-1 x[j]exp(2piij*k/n), 0<=k<n
36		<case2>
37		X[k] = sum_j=0^n-1 x[j]exp(-2piij*k/n), 0<=k<n
38		(notes: sum_j=0^n-1 is a summation from j=0 to n-1)
39		[usage]
40		<case1>
41		ip[0] = 0; // first time only
42		cdft(2*n, 1, a, ip, w);
43		<case2>
44		ip[0] = 0; // first time only
45		cdft(2*n, -1, a, ip, w);
46		[parameters]
47		2*n :data length (int)
48		n >= 1, n = power of 2
49		a[0...2n-1] :input/output data (double )
50		input data
51		a[2*j] = Re(x[j]),
52		a[2*j+1] = Im(x[j]), 0<=j<n
53		output data
54		a[2*k] = Re(X[k]),
55		a[2*k+1] = Im(X[k]), 0<=k<n
56		ip[0...] :work area for bit reversal (int )
57		length of ip >= 2+sqrt(n)
58		strictly,
59		length of ip >=
60		2+(1<<(int)(log(n+0.5)/log(2))/2).
61		ip[0],ip[1] are pointers of the cos/sin table.
62		w[0...n/2-1] :cos/sin table (double *)
63		w[],ip[] are initialized if ip[0] == 0.
64		[remark]
65		Inverse of
66		cdft(2*n, -1, a, ip, w);
67		is
68		cdft(2*n, 1, a, ip, w);
69		for (j = 0; j <= 2 * n - 1; j++) {
70		a[j] *= 1.0 / n;
71		}
72		.
73
74
75		-------- Real DFT / Inverse of Real DFT --------
76		[definition]
77		<case1> RDFT
78		R[k] = sum_j=0^n-1 a[j]cos(2pijk/n), 0<=k<=n/2
79		I[k] = sum_j=0^n-1 a[j]sin(2pijk/n), 0<k<n/2
80		<case2> IRDFT (excluding scale)
81		a[k] = (R[0] + R[n/2]cos(pik))/2 +
82		sum_j=1^n/2-1 R[j]cos(2pijk/n) +
83		sum_j=1^n/2-1 I[j]sin(2pijk/n), 0<=k<n
84		[usage]
85		<case1>
86		ip[0] = 0; // first time only
87		rdft(n, 1, a, ip, w);
88		<case2>
89		ip[0] = 0; // first time only
90		rdft(n, -1, a, ip, w);
91		[parameters]
92		n :data length (int)
93		n >= 2, n = power of 2
94		a[0...n-1] :input/output data (double *)
95		<case1>
96		output data
97		a[2*k] = R[k], 0<=k<n/2
98		a[2*k+1] = I[k], 0<k<n/2
99		a[1] = R[n/2]
100		<case2>
101		input data
102		a[2*j] = R[j], 0<=j<n/2
103		a[2*j+1] = I[j], 0<j<n/2
104		a[1] = R[n/2]
105		ip[0...] :work area for bit reversal (int )
106		length of ip >= 2+sqrt(n/2)
107		strictly,
108		length of ip >=
109		2+(1<<(int)(log(n/2+0.5)/log(2))/2).
110		ip[0],ip[1] are pointers of the cos/sin table.
111		w[0...n/2-1] :cos/sin table (double *)
112		w[],ip[] are initialized if ip[0] == 0.
113		[remark]
114		Inverse of
115		rdft(n, 1, a, ip, w);
116		is
117		rdft(n, -1, a, ip, w);
118		for (j = 0; j <= n - 1; j++) {
119		a[j] *= 2.0 / n;
120		}
121		.
122
123
124		-------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
125		[definition]
126		<case1> IDCT (excluding scale)
127		C[k] = sum_j=0^n-1 a[j]cos(pij*(k+1/2)/n), 0<=k<n
128		<case2> DCT
129		C[k] = sum_j=0^n-1 a[j]cos(pi(j+1/2)*k/n), 0<=k<n
130		[usage]
131		<case1>
132		ip[0] = 0; // first time only
133		ddct(n, 1, a, ip, w);
134		<case2>
135		ip[0] = 0; // first time only
136		ddct(n, -1, a, ip, w);
137		[parameters]
138		n :data length (int)
139		n >= 2, n = power of 2
140		a[0...n-1] :input/output data (double *)
141		output data
142		a[k] = C[k], 0<=k<n
143		ip[0...] :work area for bit reversal (int )
144		length of ip >= 2+sqrt(n/2)
145		strictly,
146		length of ip >=
147		2+(1<<(int)(log(n/2+0.5)/log(2))/2).
148		ip[0],ip[1] are pointers of the cos/sin table.
149		w[0...n5/4-1] :cos/sin table (double )
150		w[],ip[] are initialized if ip[0] == 0.
151		[remark]
152		Inverse of
153		ddct(n, -1, a, ip, w);
154		is
155		a[0] *= 0.5;
156		ddct(n, 1, a, ip, w);
157		for (j = 0; j <= n - 1; j++) {
158		a[j] *= 2.0 / n;
159		}
160		.
161
162
163		-------- DST (Discrete Sine Transform) / Inverse of DST --------
164		[definition]
165		<case1> IDST (excluding scale)
166		S[k] = sum_j=1^n A[j]sin(pij*(k+1/2)/n), 0<=k<n
167		<case2> DST
168		S[k] = sum_j=0^n-1 a[j]sin(pi(j+1/2)*k/n), 0<k<=n
169		[usage]
170		<case1>
171		ip[0] = 0; // first time only
172		ddst(n, 1, a, ip, w);
173		<case2>
174		ip[0] = 0; // first time only
175		ddst(n, -1, a, ip, w);
176		[parameters]
177		n :data length (int)
178		n >= 2, n = power of 2
179		a[0...n-1] :input/output data (double *)
180		<case1>
181		input data
182		a[j] = A[j], 0<j<n
183		a[0] = A[n]
184		output data
185		a[k] = S[k], 0<=k<n
186		<case2>
187		output data
188		a[k] = S[k], 0<k<n
189		a[0] = S[n]
190		ip[0...] :work area for bit reversal (int )
191		length of ip >= 2+sqrt(n/2)
192		strictly,
193		length of ip >=
194		2+(1<<(int)(log(n/2+0.5)/log(2))/2).
195		ip[0],ip[1] are pointers of the cos/sin table.
196		w[0...n5/4-1] :cos/sin table (double )
197		w[],ip[] are initialized if ip[0] == 0.
198		[remark]
199		Inverse of
200		ddst(n, -1, a, ip, w);
201		is
202		a[0] *= 0.5;
203		ddst(n, 1, a, ip, w);
204		for (j = 0; j <= n - 1; j++) {
205		a[j] *= 2.0 / n;
206		}
207		.
208
209
210		-------- Cosine Transform of RDFT (Real Symmetric DFT) --------
211		[definition]
212		C[k] = sum_j=0^n a[j]cos(pij*k/n), 0<=k<=n
213		[usage]
214		ip[0] = 0; // first time only
215		dfct(n, a, t, ip, w);
216		[parameters]
217		n :data length - 1 (int)
218		n >= 2, n = power of 2
219		a[0...n] :input/output data (double *)
220		output data
221		a[k] = C[k], 0<=k<=n
222		t[0...n/2] :work area (double *)
223		ip[0...] :work area for bit reversal (int )
224		length of ip >= 2+sqrt(n/4)
225		strictly,
226		length of ip >=
227		2+(1<<(int)(log(n/4+0.5)/log(2))/2).
228		ip[0],ip[1] are pointers of the cos/sin table.
229		w[0...n5/8-1] :cos/sin table (double )
230		w[],ip[] are initialized if ip[0] == 0.
231		[remark]
232		Inverse of
233		a[0] *= 0.5;
234		a[n] *= 0.5;
235		dfct(n, a, t, ip, w);
236		is
237		a[0] *= 0.5;
238		a[n] *= 0.5;
239		dfct(n, a, t, ip, w);
240		for (j = 0; j <= n; j++) {
241		a[j] *= 2.0 / n;
242		}
243		.
244
245
246		-------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
247		[definition]
248		S[k] = sum_j=1^n-1 a[j]sin(pij*k/n), 0<k<n
249		[usage]
250		ip[0] = 0; // first time only
251		dfst(n, a, t, ip, w);
252		[parameters]
253		n :data length + 1 (int)
254		n >= 2, n = power of 2
255		a[0...n-1] :input/output data (double *)
256		output data
257		a[k] = S[k], 0<k<n
258		(a[0] is used for work area)
259		t[0...n/2-1] :work area (double *)
260		ip[0...] :work area for bit reversal (int )
261		length of ip >= 2+sqrt(n/4)
262		strictly,
263		length of ip >=
264		2+(1<<(int)(log(n/4+0.5)/log(2))/2).
265		ip[0],ip[1] are pointers of the cos/sin table.
266		w[0...n5/8-1] :cos/sin table (double )
267		w[],ip[] are initialized if ip[0] == 0.
268		[remark]
269		Inverse of
270		dfst(n, a, t, ip, w);
271		is
272		dfst(n, a, t, ip, w);
273		for (j = 1; j <= n - 1; j++) {
274		a[j] *= 2.0 / n;
275		}
276		.
277
278
279		Appendix :
280		The cos/sin table is recalculated when the larger table required.
281		w[] and ip[] are compatible with all routines.
282		*/
283
284
285		void cdft(int n, int isgn, double a, int ip, double *w)
286	0	{
287	0	void makewt(int nw, int ip, double w);
288	0	void cftfsub(int n, double a, int ip, int nw, double *w);
289	0	void cftbsub(int n, double a, int ip, int nw, double *w);
290	0	int nw;
291
292	0	nw = ip[0];
293	0	if (n > (nw << 2)) {
294	0	nw = n >> 2;
295	0	makewt(nw, ip, w);
296	0	}
297	0	if (isgn >= 0) {
298	0	cftfsub(n, a, ip, nw, w);
299	0	} else {
300	0	cftbsub(n, a, ip, nw, w);
301	0	}
302	0	}
303
304
305		void rdft(int n, int isgn, double a, int ip, double *w)
306	0	{
307	0	void makewt(int nw, int ip, double w);
308	0	void makect(int nc, int ip, double c);
309	0	void cftfsub(int n, double a, int ip, int nw, double *w);
310	0	void cftbsub(int n, double a, int ip, int nw, double *w);
311	0	void rftfsub(int n, double a, int nc, double c);
312	0	void rftbsub(int n, double a, int nc, double c);
313	0	int nw, nc;
314	0	double xi;
315
316	0	nw = ip[0];
317	0	if (n > (nw << 2)) {
318	0	nw = n >> 2;
319	0	makewt(nw, ip, w);
320	0	}
321	0	nc = ip[1];
322	0	if (n > (nc << 2)) {
323	0	nc = n >> 2;
324	0	makect(nc, ip, w + nw);
325	0	}
326	0	if (isgn >= 0) {
327	0	if (n > 4) {
328	0	cftfsub(n, a, ip, nw, w);
329	0	rftfsub(n, a, nc, w + nw);
330	0	} else if (n == 4) {
331	0	cftfsub(n, a, ip, nw, w);
332	0	}
333	0	xi = a[0] - a[1];
334	0	a[0] += a[1];
335	0	a[1] = xi;
336	0	} else {
337	0	a[1] = 0.5 * (a[0] - a[1]);
338	0	a[0] -= a[1];
339	0	if (n > 4) {
340	0	rftbsub(n, a, nc, w + nw);
341	0	cftbsub(n, a, ip, nw, w);
342	0	} else if (n == 4) {
343	0	cftbsub(n, a, ip, nw, w);
344	0	}
345	0	}
346	0	}
347
348
349		void ddct(int n, int isgn, double a, int ip, double *w)
350	0	{
351	0	void makewt(int nw, int ip, double w);
352	0	void makect(int nc, int ip, double c);
353	0	void cftfsub(int n, double a, int ip, int nw, double *w);
354	0	void cftbsub(int n, double a, int ip, int nw, double *w);
355	0	void rftfsub(int n, double a, int nc, double c);
356	0	void rftbsub(int n, double a, int nc, double c);
357	0	void dctsub(int n, double a, int nc, double c);
358	0	int j, nw, nc;
359	0	double xr;
360
361	0	nw = ip[0];
362	0	if (n > (nw << 2)) {
363	0	nw = n >> 2;
364	0	makewt(nw, ip, w);
365	0	}
366	0	nc = ip[1];
367	0	if (n > nc) {
368	0	nc = n;
369	0	makect(nc, ip, w + nw);
370	0	}
371	0	if (isgn < 0) {
372	0	xr = a[n - 1];
373	0	for (j = n - 2; j >= 2; j -= 2) {
374	0	a[j + 1] = a[j] - a[j - 1];
375	0	a[j] += a[j - 1];
376	0	}
377	0	a[1] = a[0] - xr;
378	0	a[0] += xr;
379	0	if (n > 4) {
380	0	rftbsub(n, a, nc, w + nw);
381	0	cftbsub(n, a, ip, nw, w);
382	0	} else if (n == 4) {
383	0	cftbsub(n, a, ip, nw, w);
384	0	}
385	0	}
386	0	dctsub(n, a, nc, w + nw);
387	0	if (isgn >= 0) {
388	0	if (n > 4) {
389	0	cftfsub(n, a, ip, nw, w);
390	0	rftfsub(n, a, nc, w + nw);
391	0	} else if (n == 4) {
392	0	cftfsub(n, a, ip, nw, w);
393	0	}
394	0	xr = a[0] - a[1];
395	0	a[0] += a[1];
396	0	for (j = 2; j < n; j += 2) {
397	0	a[j - 1] = a[j] - a[j + 1];
398	0	a[j] += a[j + 1];
399	0	}
400	0	a[n - 1] = xr;
401	0	}
402	0	}
403
404
405		void ddst(int n, int isgn, double a, int ip, double *w)
406	0	{
407	0	void makewt(int nw, int ip, double w);
408	0	void makect(int nc, int ip, double c);
409	0	void cftfsub(int n, double a, int ip, int nw, double *w);
410	0	void cftbsub(int n, double a, int ip, int nw, double *w);
411	0	void rftfsub(int n, double a, int nc, double c);
412	0	void rftbsub(int n, double a, int nc, double c);
413	0	void dstsub(int n, double a, int nc, double c);
414	0	int j, nw, nc;
415	0	double xr;
416
417	0	nw = ip[0];
418	0	if (n > (nw << 2)) {
419	0	nw = n >> 2;
420	0	makewt(nw, ip, w);
421	0	}
422	0	nc = ip[1];
423	0	if (n > nc) {
424	0	nc = n;
425	0	makect(nc, ip, w + nw);
426	0	}
427	0	if (isgn < 0) {
428	0	xr = a[n - 1];
429	0	for (j = n - 2; j >= 2; j -= 2) {
430	0	a[j + 1] = -a[j] - a[j - 1];
431	0	a[j] -= a[j - 1];
432	0	}
433	0	a[1] = a[0] + xr;
434	0	a[0] -= xr;
435	0	if (n > 4) {
436	0	rftbsub(n, a, nc, w + nw);
437	0	cftbsub(n, a, ip, nw, w);
438	0	} else if (n == 4) {
439	0	cftbsub(n, a, ip, nw, w);
440	0	}
441	0	}
442	0	dstsub(n, a, nc, w + nw);
443	0	if (isgn >= 0) {
444	0	if (n > 4) {
445	0	cftfsub(n, a, ip, nw, w);
446	0	rftfsub(n, a, nc, w + nw);
447	0	} else if (n == 4) {
448	0	cftfsub(n, a, ip, nw, w);
449	0	}
450	0	xr = a[0] - a[1];
451	0	a[0] += a[1];
452	0	for (j = 2; j < n; j += 2) {
453	0	a[j - 1] = -a[j] - a[j + 1];
454	0	a[j] -= a[j + 1];
455	0	}
456	0	a[n - 1] = -xr;
457	0	}
458	0	}
459
460
461		void dfct(int n, double a, double t, int ip, double w)
462	0	{
463	0	void makewt(int nw, int ip, double w);
464	0	void makect(int nc, int ip, double c);
465	0	void cftfsub(int n, double a, int ip, int nw, double *w);
466	0	void rftfsub(int n, double a, int nc, double c);
467	0	void dctsub(int n, double a, int nc, double c);
468	0	int j, k, l, m, mh, nw, nc;
469	0	double xr, xi, yr, yi;
470
471	0	nw = ip[0];
472	0	if (n > (nw << 3)) {
473	0	nw = n >> 3;
474	0	makewt(nw, ip, w);
475	0	}
476	0	nc = ip[1];
477	0	if (n > (nc << 1)) {
478	0	nc = n >> 1;
479	0	makect(nc, ip, w + nw);
480	0	}
481	0	m = n >> 1;
482	0	yi = a[m];
483	0	xi = a[0] + a[n];
484	0	a[0] -= a[n];
485	0	t[0] = xi - yi;
486	0	t[m] = xi + yi;
487	0	if (n > 2) {
488	0	mh = m >> 1;
489	0	for (j = 1; j < mh; j++) {
490	0	k = m - j;
491	0	xr = a[j] - a[n - j];
492	0	xi = a[j] + a[n - j];
493	0	yr = a[k] - a[n - k];
494	0	yi = a[k] + a[n - k];
495	0	a[j] = xr;
496	0	a[k] = yr;
497	0	t[j] = xi - yi;
498	0	t[k] = xi + yi;
499	0	}
500	0	t[mh] = a[mh] + a[n - mh];
501	0	a[mh] -= a[n - mh];
502	0	dctsub(m, a, nc, w + nw);
503	0	if (m > 4) {
504	0	cftfsub(m, a, ip, nw, w);
505	0	rftfsub(m, a, nc, w + nw);
506	0	} else if (m == 4) {
507	0	cftfsub(m, a, ip, nw, w);
508	0	}
509	0	a[n - 1] = a[0] - a[1];
510	0	a[1] = a[0] + a[1];
511	0	for (j = m - 2; j >= 2; j -= 2) {
512	0	a[2 * j + 1] = a[j] + a[j + 1];
513	0	a[2 * j - 1] = a[j] - a[j + 1];
514	0	}
515	0	l = 2;
516	0	m = mh;
517	0	while (m >= 2) {
518	0	dctsub(m, t, nc, w + nw);
519	0	if (m > 4) {
520	0	cftfsub(m, t, ip, nw, w);
521	0	rftfsub(m, t, nc, w + nw);
522	0	} else if (m == 4) {
523	0	cftfsub(m, t, ip, nw, w);
524	0	}
525	0	a[n - l] = t[0] - t[1];
526	0	a[l] = t[0] + t[1];
527	0	k = 0;
528	0	for (j = 2; j < m; j += 2) {
529	0	k += l << 2;
530	0	a[k - l] = t[j] - t[j + 1];
531	0	a[k + l] = t[j] + t[j + 1];
532	0	}
533	0	l <<= 1;
534	0	mh = m >> 1;
535	0	for (j = 0; j < mh; j++) {
536	0	k = m - j;
537	0	t[j] = t[m + k] - t[m + j];
538	0	t[k] = t[m + k] + t[m + j];
539	0	}
540	0	t[mh] = t[m + mh];
541	0	m = mh;
542	0	}
543	0	a[l] = t[0];
544	0	a[n] = t[2] - t[1];
545	0	a[0] = t[2] + t[1];
546	0	} else {
547	0	a[1] = a[0];
548	0	a[2] = t[0];
549	0	a[0] = t[1];
550	0	}
551	0	}
552
553
554		void dfst(int n, double a, double t, int ip, double w)
555	0	{
556	0	void makewt(int nw, int ip, double w);
557	0	void makect(int nc, int ip, double c);
558	0	void cftfsub(int n, double a, int ip, int nw, double *w);
559	0	void rftfsub(int n, double a, int nc, double c);
560	0	void dstsub(int n, double a, int nc, double c);
561	0	int j, k, l, m, mh, nw, nc;
562	0	double xr, xi, yr, yi;
563
564	0	nw = ip[0];
565	0	if (n > (nw << 3)) {
566	0	nw = n >> 3;
567	0	makewt(nw, ip, w);
568	0	}
569	0	nc = ip[1];
570	0	if (n > (nc << 1)) {
571	0	nc = n >> 1;
572	0	makect(nc, ip, w + nw);
573	0	}
574	0	if (n > 2) {
575	0	m = n >> 1;
576	0	mh = m >> 1;
577	0	for (j = 1; j < mh; j++) {
578	0	k = m - j;
579	0	xr = a[j] + a[n - j];
580	0	xi = a[j] - a[n - j];
581	0	yr = a[k] + a[n - k];
582	0	yi = a[k] - a[n - k];
583	0	a[j] = xr;
584	0	a[k] = yr;
585	0	t[j] = xi + yi;
586	0	t[k] = xi - yi;
587	0	}
588	0	t[0] = a[mh] - a[n - mh];
589	0	a[mh] += a[n - mh];
590	0	a[0] = a[m];
591	0	dstsub(m, a, nc, w + nw);
592	0	if (m > 4) {
593	0	cftfsub(m, a, ip, nw, w);
594	0	rftfsub(m, a, nc, w + nw);
595	0	} else if (m == 4) {
596	0	cftfsub(m, a, ip, nw, w);
597	0	}
598	0	a[n - 1] = a[1] - a[0];
599	0	a[1] = a[0] + a[1];
600	0	for (j = m - 2; j >= 2; j -= 2) {
601	0	a[2 * j + 1] = a[j] - a[j + 1];
602	0	a[2 * j - 1] = -a[j] - a[j + 1];
603	0	}
604	0	l = 2;
605	0	m = mh;
606	0	while (m >= 2) {
607	0	dstsub(m, t, nc, w + nw);
608	0	if (m > 4) {
609	0	cftfsub(m, t, ip, nw, w);
610	0	rftfsub(m, t, nc, w + nw);
611	0	} else if (m == 4) {
612	0	cftfsub(m, t, ip, nw, w);
613	0	}
614	0	a[n - l] = t[1] - t[0];
615	0	a[l] = t[0] + t[1];
616	0	k = 0;
617	0	for (j = 2; j < m; j += 2) {
618	0	k += l << 2;
619	0	a[k - l] = -t[j] - t[j + 1];
620	0	a[k + l] = t[j] - t[j + 1];
621	0	}
622	0	l <<= 1;
623	0	mh = m >> 1;
624	0	for (j = 1; j < mh; j++) {
625	0	k = m - j;
626	0	t[j] = t[m + k] + t[m + j];
627	0	t[k] = t[m + k] - t[m + j];
628	0	}
629	0	t[0] = t[m + mh];
630	0	m = mh;
631	0	}
632	0	a[l] = t[0];
633	0	}
634	0	a[0] = 0;
635	0	}
636
637
638		/* -------- initializing routines -------- */
639
640
641		#include <math.h>
642
643		void makewt(int nw, int ip, double w)
644	0	{
645	0	void makeipt(int nw, int *ip);
646	0	int j, nwh, nw0, nw1;
647	0	double delta, wn4r, wk1r, wk1i, wk3r, wk3i;
648
649	0	ip[0] = nw;
650	0	ip[1] = 1;
651	0	if (nw > 2) {
652	0	nwh = nw >> 1;
653	0	delta = atan(1.0) / nwh;
654	0	wn4r = cos(delta * nwh);
655	0	w[0] = 1;
656	0	w[1] = wn4r;
657	0	if (nwh == 4) {
658	0	w[2] = cos(delta * 2);
659	0	w[3] = sin(delta * 2);
660	0	} else if (nwh > 4) {
661	0	makeipt(nw, ip);
662	0	w[2] = 0.5 / cos(delta * 2);
663	0	w[3] = 0.5 / cos(delta * 6);
664	0	for (j = 4; j < nwh; j += 4) {
665	0	w[j] = cos(delta * j);
666	0	w[j + 1] = sin(delta * j);
667	0	w[j + 2] = cos(3 * delta * j);
668	0	w[j + 3] = -sin(3 * delta * j);
669	0	}
670	0	}
671	0	nw0 = 0;
672	0	while (nwh > 2) {
673	0	nw1 = nw0 + nwh;
674	0	nwh >>= 1;
675	0	w[nw1] = 1;
676	0	w[nw1 + 1] = wn4r;
677	0	if (nwh == 4) {
678	0	wk1r = w[nw0 + 4];
679	0	wk1i = w[nw0 + 5];
680	0	w[nw1 + 2] = wk1r;
681	0	w[nw1 + 3] = wk1i;
682	0	} else if (nwh > 4) {
683	0	wk1r = w[nw0 + 4];
684	0	wk3r = w[nw0 + 6];
685	0	w[nw1 + 2] = 0.5 / wk1r;
686	0	w[nw1 + 3] = 0.5 / wk3r;
687	0	for (j = 4; j < nwh; j += 4) {
688	0	wk1r = w[nw0 + 2 * j];
689	0	wk1i = w[nw0 + 2 * j + 1];
690	0	wk3r = w[nw0 + 2 * j + 2];
691	0	wk3i = w[nw0 + 2 * j + 3];
692	0	w[nw1 + j] = wk1r;
693	0	w[nw1 + j + 1] = wk1i;
694	0	w[nw1 + j + 2] = wk3r;
695	0	w[nw1 + j + 3] = wk3i;
696	0	}
697	0	}
698	0	nw0 = nw1;
699	0	}
700	0	}
701	0	}
702
703
704		void makeipt(int nw, int *ip)
705	0	{
706	0	int j, l, m, m2, p, q;
707
708	0	ip[2] = 0;
709	0	ip[3] = 16;
710	0	m = 2;
711	0	for (l = nw; l > 32; l >>= 2) {
712	0	m2 = m << 1;
713	0	q = m2 << 3;
714	0	for (j = m; j < m2; j++) {
715	0	p = ip[j] << 2;
716	0	ip[m + j] = p;
717	0	ip[m2 + j] = p + q;
718	0	}
719	0	m = m2;
720	0	}
721	0	}
722
723
724		void makect(int nc, int ip, double c)
725	0	{
726	0	int j, nch;
727	0	double delta;
728
729	0	ip[1] = nc;
730	0	if (nc > 1) {
731	0	nch = nc >> 1;
732	0	delta = atan(1.0) / nch;
733	0	c[0] = cos(delta * nch);
734	0	c[nch] = 0.5 * c[0];
735	0	for (j = 1; j < nch; j++) {
736	0	c[j] = 0.5 * cos(delta * j);
737	0	c[nc - j] = 0.5 * sin(delta * j);
738	0	}
739	0	}
740	0	}
741
742
743		/* -------- child routines -------- */
744
745
746		#ifdef USE_CDFT_PTHREADS
747		#define USE_CDFT_THREADS
748		#ifndef CDFT_THREADS_BEGIN_N
749		#define CDFT_THREADS_BEGIN_N 8192
750		#endif
751		#ifndef CDFT_4THREADS_BEGIN_N
752		#define CDFT_4THREADS_BEGIN_N 65536
753		#endif
754		#include <pthread.h>
755		#include <stdio.h>
756		#include <stdlib.h>
757		#define cdft_thread_t pthread_t
758		#define cdft_thread_create(thp,func,argp) { \
759		if (pthread_create(thp, NULL, func, (void *) argp) != 0) { \
760		fprintf(stderr, "cdft thread error\n"); \
761		exit(1); \
762		} \
763		}
764		#define cdft_thread_wait(th) { \
765		if (pthread_join(th, NULL) != 0) { \
766		fprintf(stderr, "cdft thread error\n"); \
767		exit(1); \
768		} \
769		}
770		#endif /* USE_CDFT_PTHREADS */
771
772
773		#ifdef USE_CDFT_WINTHREADS
774		#define USE_CDFT_THREADS
775		#ifndef CDFT_THREADS_BEGIN_N
776		#define CDFT_THREADS_BEGIN_N 32768
777		#endif
778		#ifndef CDFT_4THREADS_BEGIN_N
779		#define CDFT_4THREADS_BEGIN_N 524288
780		#endif
781		#include <windows.h>
782		#include <stdio.h>
783		#include <stdlib.h>
784		#define cdft_thread_t HANDLE
785		#define cdft_thread_create(thp,func,argp) { \
786		DWORD thid; \
787		*(thp) = CreateThread(NULL, 0, (LPTHREAD_START_ROUTINE) func, (LPVOID) argp, 0, &thid); \
788		if (*(thp) == 0) { \
789		fprintf(stderr, "cdft thread error\n"); \
790		exit(1); \
791		} \
792		}
793		#define cdft_thread_wait(th) { \
794		WaitForSingleObject(th, INFINITE); \
795		CloseHandle(th); \
796		}
797		#endif /* USE_CDFT_WINTHREADS */
798
799
800		void cftfsub(int n, double a, int ip, int nw, double *w)
801	0	{
802	0	void bitrv2(int n, int ip, double a);
803	0	void bitrv216(double *a);
804	0	void bitrv208(double *a);
805	0	void cftf1st(int n, double a, double w);
806	0	void cftrec4(int n, double a, int nw, double w);
807	0	void cftleaf(int n, int isplt, double a, int nw, double w);
808	0	void cftfx41(int n, double a, int nw, double w);
809	0	void cftf161(double a, double w);
810	0	void cftf081(double a, double w);
811	0	void cftf040(double *a);
812	0	void cftx020(double *a);
813		#ifdef USE_CDFT_THREADS
814		void cftrec4_th(int n, double a, int nw, double w);
815		#endif /* USE_CDFT_THREADS */
816
817	0	if (n > 8) {
818	0	if (n > 32) {
819	0	cftf1st(n, a, &w[nw - (n >> 2)]);
820		#ifdef USE_CDFT_THREADS
821		if (n > CDFT_THREADS_BEGIN_N) {
822		cftrec4_th(n, a, nw, w);
823		} else
824		#endif /* USE_CDFT_THREADS */
825	0	if (n > 512) {
826	0	cftrec4(n, a, nw, w);
827	0	} else if (n > 128) {
828	0	cftleaf(n, 1, a, nw, w);
829	0	} else {
830	0	cftfx41(n, a, nw, w);
831	0	}
832	0	bitrv2(n, ip, a);
833	0	} else if (n == 32) {
834	0	cftf161(a, &w[nw - 8]);
835	0	bitrv216(a);
836	0	} else {
837	0	cftf081(a, w);
838	0	bitrv208(a);
839	0	}
840	0	} else if (n == 8) {
841	0	cftf040(a);
842	0	} else if (n == 4) {
843	0	cftx020(a);
844	0	}
845	0	}
846
847
848		void cftbsub(int n, double a, int ip, int nw, double *w)
849	0	{
850	0	void bitrv2conj(int n, int ip, double a);
851	0	void bitrv216neg(double *a);
852	0	void bitrv208neg(double *a);
853	0	void cftb1st(int n, double a, double w);
854	0	void cftrec4(int n, double a, int nw, double w);
855	0	void cftleaf(int n, int isplt, double a, int nw, double w);
856	0	void cftfx41(int n, double a, int nw, double w);
857	0	void cftf161(double a, double w);
858	0	void cftf081(double a, double w);
859	0	void cftb040(double *a);
860	0	void cftx020(double *a);
861		#ifdef USE_CDFT_THREADS
862		void cftrec4_th(int n, double a, int nw, double w);
863		#endif /* USE_CDFT_THREADS */
864
865	0	if (n > 8) {
866	0	if (n > 32) {
867	0	cftb1st(n, a, &w[nw - (n >> 2)]);
868		#ifdef USE_CDFT_THREADS
869		if (n > CDFT_THREADS_BEGIN_N) {
870		cftrec4_th(n, a, nw, w);
871		} else
872		#endif /* USE_CDFT_THREADS */
873	0	if (n > 512) {
874	0	cftrec4(n, a, nw, w);
875	0	} else if (n > 128) {
876	0	cftleaf(n, 1, a, nw, w);
877	0	} else {
878	0	cftfx41(n, a, nw, w);
879	0	}
880	0	bitrv2conj(n, ip, a);
881	0	} else if (n == 32) {
882	0	cftf161(a, &w[nw - 8]);
883	0	bitrv216neg(a);
884	0	} else {
885	0	cftf081(a, w);
886	0	bitrv208neg(a);
887	0	}
888	0	} else if (n == 8) {
889	0	cftb040(a);
890	0	} else if (n == 4) {
891	0	cftx020(a);
892	0	}
893	0	}
894
895
896		void bitrv2(int n, int ip, double a)
897	0	{
898	0	int j, j1, k, k1, l, m, nh, nm;
899	0	double xr, xi, yr, yi;
900
901	0	m = 1;
902	0	for (l = n >> 2; l > 8; l >>= 2) {
903	0	m <<= 1;
904	0	}
905	0	nh = n >> 1;
906	0	nm = 4 * m;
907	0	if (l == 8) {
908	0	for (k = 0; k < m; k++) {
909	0	for (j = 0; j < k; j++) {
910	0	j1 = 4 * j + 2 * ip[m + k];
911	0	k1 = 4 * k + 2 * ip[m + j];
912	0	xr = a[j1];
913	0	xi = a[j1 + 1];
914	0	yr = a[k1];
915	0	yi = a[k1 + 1];
916	0	a[j1] = yr;
917	0	a[j1 + 1] = yi;
918	0	a[k1] = xr;
919	0	a[k1 + 1] = xi;
920	0	j1 += nm;
921	0	k1 += 2 * nm;
922	0	xr = a[j1];
923	0	xi = a[j1 + 1];
924	0	yr = a[k1];
925	0	yi = a[k1 + 1];
926	0	a[j1] = yr;
927	0	a[j1 + 1] = yi;
928	0	a[k1] = xr;
929	0	a[k1 + 1] = xi;
930	0	j1 += nm;
931	0	k1 -= nm;
932	0	xr = a[j1];
933	0	xi = a[j1 + 1];
934	0	yr = a[k1];
935	0	yi = a[k1 + 1];
936	0	a[j1] = yr;
937	0	a[j1 + 1] = yi;
938	0	a[k1] = xr;
939	0	a[k1 + 1] = xi;
940	0	j1 += nm;
941	0	k1 += 2 * nm;
942	0	xr = a[j1];
943	0	xi = a[j1 + 1];
944	0	yr = a[k1];
945	0	yi = a[k1 + 1];
946	0	a[j1] = yr;
947	0	a[j1 + 1] = yi;
948	0	a[k1] = xr;
949	0	a[k1 + 1] = xi;
950	0	j1 += nh;
951	0	k1 += 2;
952	0	xr = a[j1];
953	0	xi = a[j1 + 1];
954	0	yr = a[k1];
955	0	yi = a[k1 + 1];
956	0	a[j1] = yr;
957	0	a[j1 + 1] = yi;
958	0	a[k1] = xr;
959	0	a[k1 + 1] = xi;
960	0	j1 -= nm;
961	0	k1 -= 2 * nm;
962	0	xr = a[j1];
963	0	xi = a[j1 + 1];
964	0	yr = a[k1];
965	0	yi = a[k1 + 1];
966	0	a[j1] = yr;
967	0	a[j1 + 1] = yi;
968	0	a[k1] = xr;
969	0	a[k1 + 1] = xi;
970	0	j1 -= nm;
971	0	k1 += nm;
972	0	xr = a[j1];
973	0	xi = a[j1 + 1];
974	0	yr = a[k1];
975	0	yi = a[k1 + 1];
976	0	a[j1] = yr;
977	0	a[j1 + 1] = yi;
978	0	a[k1] = xr;
979	0	a[k1 + 1] = xi;
980	0	j1 -= nm;
981	0	k1 -= 2 * nm;
982	0	xr = a[j1];
983	0	xi = a[j1 + 1];
984	0	yr = a[k1];
985	0	yi = a[k1 + 1];
986	0	a[j1] = yr;
987	0	a[j1 + 1] = yi;
988	0	a[k1] = xr;
989	0	a[k1 + 1] = xi;
990	0	j1 += 2;
991	0	k1 += nh;
992	0	xr = a[j1];
993	0	xi = a[j1 + 1];
994	0	yr = a[k1];
995	0	yi = a[k1 + 1];
996	0	a[j1] = yr;
997	0	a[j1 + 1] = yi;
998	0	a[k1] = xr;
999	0	a[k1 + 1] = xi;
1000	0	j1 += nm;
1001	0	k1 += 2 * nm;
1002	0	xr = a[j1];
1003	0	xi = a[j1 + 1];
1004	0	yr = a[k1];
1005	0	yi = a[k1 + 1];
1006	0	a[j1] = yr;
1007	0	a[j1 + 1] = yi;
1008	0	a[k1] = xr;
1009	0	a[k1 + 1] = xi;
1010	0	j1 += nm;
1011	0	k1 -= nm;
1012	0	xr = a[j1];
1013	0	xi = a[j1 + 1];
1014	0	yr = a[k1];
1015	0	yi = a[k1 + 1];
1016	0	a[j1] = yr;
1017	0	a[j1 + 1] = yi;
1018	0	a[k1] = xr;
1019	0	a[k1 + 1] = xi;
1020	0	j1 += nm;
1021	0	k1 += 2 * nm;
1022	0	xr = a[j1];
1023	0	xi = a[j1 + 1];
1024	0	yr = a[k1];
1025	0	yi = a[k1 + 1];
1026	0	a[j1] = yr;
1027	0	a[j1 + 1] = yi;
1028	0	a[k1] = xr;
1029	0	a[k1 + 1] = xi;
1030	0	j1 -= nh;
1031	0	k1 -= 2;
1032	0	xr = a[j1];
1033	0	xi = a[j1 + 1];
1034	0	yr = a[k1];
1035	0	yi = a[k1 + 1];
1036	0	a[j1] = yr;
1037	0	a[j1 + 1] = yi;
1038	0	a[k1] = xr;
1039	0	a[k1 + 1] = xi;
1040	0	j1 -= nm;
1041	0	k1 -= 2 * nm;
1042	0	xr = a[j1];
1043	0	xi = a[j1 + 1];
1044	0	yr = a[k1];
1045	0	yi = a[k1 + 1];
1046	0	a[j1] = yr;
1047	0	a[j1 + 1] = yi;
1048	0	a[k1] = xr;
1049	0	a[k1 + 1] = xi;
1050	0	j1 -= nm;
1051	0	k1 += nm;
1052	0	xr = a[j1];
1053	0	xi = a[j1 + 1];
1054	0	yr = a[k1];
1055	0	yi = a[k1 + 1];
1056	0	a[j1] = yr;
1057	0	a[j1 + 1] = yi;
1058	0	a[k1] = xr;
1059	0	a[k1 + 1] = xi;
1060	0	j1 -= nm;
1061	0	k1 -= 2 * nm;
1062	0	xr = a[j1];
1063	0	xi = a[j1 + 1];
1064	0	yr = a[k1];
1065	0	yi = a[k1 + 1];
1066	0	a[j1] = yr;
1067	0	a[j1 + 1] = yi;
1068	0	a[k1] = xr;
1069	0	a[k1 + 1] = xi;
1070	0	}
1071	0	k1 = 4 * k + 2 * ip[m + k];
1072	0	j1 = k1 + 2;
1073	0	k1 += nh;
1074	0	xr = a[j1];
1075	0	xi = a[j1 + 1];
1076	0	yr = a[k1];
1077	0	yi = a[k1 + 1];
1078	0	a[j1] = yr;
1079	0	a[j1 + 1] = yi;
1080	0	a[k1] = xr;
1081	0	a[k1 + 1] = xi;
1082	0	j1 += nm;
1083	0	k1 += 2 * nm;
1084	0	xr = a[j1];
1085	0	xi = a[j1 + 1];
1086	0	yr = a[k1];
1087	0	yi = a[k1 + 1];
1088	0	a[j1] = yr;
1089	0	a[j1 + 1] = yi;
1090	0	a[k1] = xr;
1091	0	a[k1 + 1] = xi;
1092	0	j1 += nm;
1093	0	k1 -= nm;
1094	0	xr = a[j1];
1095	0	xi = a[j1 + 1];
1096	0	yr = a[k1];
1097	0	yi = a[k1 + 1];
1098	0	a[j1] = yr;
1099	0	a[j1 + 1] = yi;
1100	0	a[k1] = xr;
1101	0	a[k1 + 1] = xi;
1102	0	j1 -= 2;
1103	0	k1 -= nh;
1104	0	xr = a[j1];
1105	0	xi = a[j1 + 1];
1106	0	yr = a[k1];
1107	0	yi = a[k1 + 1];
1108	0	a[j1] = yr;
1109	0	a[j1 + 1] = yi;
1110	0	a[k1] = xr;
1111	0	a[k1 + 1] = xi;
1112	0	j1 += nh + 2;
1113	0	k1 += nh + 2;
1114	0	xr = a[j1];
1115	0	xi = a[j1 + 1];
1116	0	yr = a[k1];
1117	0	yi = a[k1 + 1];
1118	0	a[j1] = yr;
1119	0	a[j1 + 1] = yi;
1120	0	a[k1] = xr;
1121	0	a[k1 + 1] = xi;
1122	0	j1 -= nh - nm;
1123	0	k1 += 2 * nm - 2;
1124	0	xr = a[j1];
1125	0	xi = a[j1 + 1];
1126	0	yr = a[k1];
1127	0	yi = a[k1 + 1];
1128	0	a[j1] = yr;
1129	0	a[j1 + 1] = yi;
1130	0	a[k1] = xr;
1131	0	a[k1 + 1] = xi;
1132	0	}
1133	0	} else {
1134	0	for (k = 0; k < m; k++) {
1135	0	for (j = 0; j < k; j++) {
1136	0	j1 = 4 * j + ip[m + k];
1137	0	k1 = 4 * k + ip[m + j];
1138	0	xr = a[j1];
1139	0	xi = a[j1 + 1];
1140	0	yr = a[k1];
1141	0	yi = a[k1 + 1];
1142	0	a[j1] = yr;
1143	0	a[j1 + 1] = yi;
1144	0	a[k1] = xr;
1145	0	a[k1 + 1] = xi;
1146	0	j1 += nm;
1147	0	k1 += nm;
1148	0	xr = a[j1];
1149	0	xi = a[j1 + 1];
1150	0	yr = a[k1];
1151	0	yi = a[k1 + 1];
1152	0	a[j1] = yr;
1153	0	a[j1 + 1] = yi;
1154	0	a[k1] = xr;
1155	0	a[k1 + 1] = xi;
1156	0	j1 += nh;
1157	0	k1 += 2;
1158	0	xr = a[j1];
1159	0	xi = a[j1 + 1];
1160	0	yr = a[k1];
1161	0	yi = a[k1 + 1];
1162	0	a[j1] = yr;
1163	0	a[j1 + 1] = yi;
1164	0	a[k1] = xr;
1165	0	a[k1 + 1] = xi;
1166	0	j1 -= nm;
1167	0	k1 -= nm;
1168	0	xr = a[j1];
1169	0	xi = a[j1 + 1];
1170	0	yr = a[k1];
1171	0	yi = a[k1 + 1];
1172	0	a[j1] = yr;
1173	0	a[j1 + 1] = yi;
1174	0	a[k1] = xr;
1175	0	a[k1 + 1] = xi;
1176	0	j1 += 2;
1177	0	k1 += nh;
1178	0	xr = a[j1];
1179	0	xi = a[j1 + 1];
1180	0	yr = a[k1];
1181	0	yi = a[k1 + 1];
1182	0	a[j1] = yr;
1183	0	a[j1 + 1] = yi;
1184	0	a[k1] = xr;
1185	0	a[k1 + 1] = xi;
1186	0	j1 += nm;
1187	0	k1 += nm;
1188	0	xr = a[j1];
1189	0	xi = a[j1 + 1];
1190	0	yr = a[k1];
1191	0	yi = a[k1 + 1];
1192	0	a[j1] = yr;
1193	0	a[j1 + 1] = yi;
1194	0	a[k1] = xr;
1195	0	a[k1 + 1] = xi;
1196	0	j1 -= nh;
1197	0	k1 -= 2;
1198	0	xr = a[j1];
1199	0	xi = a[j1 + 1];
1200	0	yr = a[k1];
1201	0	yi = a[k1 + 1];
1202	0	a[j1] = yr;
1203	0	a[j1 + 1] = yi;
1204	0	a[k1] = xr;
1205	0	a[k1 + 1] = xi;
1206	0	j1 -= nm;
1207	0	k1 -= nm;
1208	0	xr = a[j1];
1209	0	xi = a[j1 + 1];
1210	0	yr = a[k1];
1211	0	yi = a[k1 + 1];
1212	0	a[j1] = yr;
1213	0	a[j1 + 1] = yi;
1214	0	a[k1] = xr;
1215	0	a[k1 + 1] = xi;
1216	0	}
1217	0	k1 = 4 * k + ip[m + k];
1218	0	j1 = k1 + 2;
1219	0	k1 += nh;
1220	0	xr = a[j1];
1221	0	xi = a[j1 + 1];
1222	0	yr = a[k1];
1223	0	yi = a[k1 + 1];
1224	0	a[j1] = yr;
1225	0	a[j1 + 1] = yi;
1226	0	a[k1] = xr;
1227	0	a[k1 + 1] = xi;
1228	0	j1 += nm;
1229	0	k1 += nm;
1230	0	xr = a[j1];
1231	0	xi = a[j1 + 1];
1232	0	yr = a[k1];
1233	0	yi = a[k1 + 1];
1234	0	a[j1] = yr;
1235	0	a[j1 + 1] = yi;
1236	0	a[k1] = xr;
1237	0	a[k1 + 1] = xi;
1238	0	}
1239	0	}
1240	0	}
1241
1242
1243		void bitrv2conj(int n, int ip, double a)
1244	0	{
1245	0	int j, j1, k, k1, l, m, nh, nm;
1246	0	double xr, xi, yr, yi;
1247
1248	0	m = 1;
1249	0	for (l = n >> 2; l > 8; l >>= 2) {
1250	0	m <<= 1;
1251	0	}
1252	0	nh = n >> 1;
1253	0	nm = 4 * m;
1254	0	if (l == 8) {
1255	0	for (k = 0; k < m; k++) {
1256	0	for (j = 0; j < k; j++) {
1257	0	j1 = 4 * j + 2 * ip[m + k];
1258	0	k1 = 4 * k + 2 * ip[m + j];
1259	0	xr = a[j1];
1260	0	xi = -a[j1 + 1];
1261	0	yr = a[k1];
1262	0	yi = -a[k1 + 1];
1263	0	a[j1] = yr;
1264	0	a[j1 + 1] = yi;
1265	0	a[k1] = xr;
1266	0	a[k1 + 1] = xi;
1267	0	j1 += nm;
1268	0	k1 += 2 * nm;
1269	0	xr = a[j1];
1270	0	xi = -a[j1 + 1];
1271	0	yr = a[k1];
1272	0	yi = -a[k1 + 1];
1273	0	a[j1] = yr;
1274	0	a[j1 + 1] = yi;
1275	0	a[k1] = xr;
1276	0	a[k1 + 1] = xi;
1277	0	j1 += nm;
1278	0	k1 -= nm;
1279	0	xr = a[j1];
1280	0	xi = -a[j1 + 1];
1281	0	yr = a[k1];
1282	0	yi = -a[k1 + 1];
1283	0	a[j1] = yr;
1284	0	a[j1 + 1] = yi;
1285	0	a[k1] = xr;
1286	0	a[k1 + 1] = xi;
1287	0	j1 += nm;
1288	0	k1 += 2 * nm;
1289	0	xr = a[j1];
1290	0	xi = -a[j1 + 1];
1291	0	yr = a[k1];
1292	0	yi = -a[k1 + 1];
1293	0	a[j1] = yr;
1294	0	a[j1 + 1] = yi;
1295	0	a[k1] = xr;
1296	0	a[k1 + 1] = xi;
1297	0	j1 += nh;
1298	0	k1 += 2;
1299	0	xr = a[j1];
1300	0	xi = -a[j1 + 1];
1301	0	yr = a[k1];
1302	0	yi = -a[k1 + 1];
1303	0	a[j1] = yr;
1304	0	a[j1 + 1] = yi;
1305	0	a[k1] = xr;
1306	0	a[k1 + 1] = xi;
1307	0	j1 -= nm;
1308	0	k1 -= 2 * nm;
1309	0	xr = a[j1];
1310	0	xi = -a[j1 + 1];
1311	0	yr = a[k1];
1312	0	yi = -a[k1 + 1];
1313	0	a[j1] = yr;
1314	0	a[j1 + 1] = yi;
1315	0	a[k1] = xr;
1316	0	a[k1 + 1] = xi;
1317	0	j1 -= nm;
1318	0	k1 += nm;
1319	0	xr = a[j1];
1320	0	xi = -a[j1 + 1];
1321	0	yr = a[k1];
1322	0	yi = -a[k1 + 1];
1323	0	a[j1] = yr;
1324	0	a[j1 + 1] = yi;
1325	0	a[k1] = xr;
1326	0	a[k1 + 1] = xi;
1327	0	j1 -= nm;
1328	0	k1 -= 2 * nm;
1329	0	xr = a[j1];
1330	0	xi = -a[j1 + 1];
1331	0	yr = a[k1];
1332	0	yi = -a[k1 + 1];
1333	0	a[j1] = yr;
1334	0	a[j1 + 1] = yi;
1335	0	a[k1] = xr;
1336	0	a[k1 + 1] = xi;
1337	0	j1 += 2;
1338	0	k1 += nh;
1339	0	xr = a[j1];
1340	0	xi = -a[j1 + 1];
1341	0	yr = a[k1];
1342	0	yi = -a[k1 + 1];
1343	0	a[j1] = yr;
1344	0	a[j1 + 1] = yi;
1345	0	a[k1] = xr;
1346	0	a[k1 + 1] = xi;
1347	0	j1 += nm;
1348	0	k1 += 2 * nm;
1349	0	xr = a[j1];
1350	0	xi = -a[j1 + 1];
1351	0	yr = a[k1];
1352	0	yi = -a[k1 + 1];
1353	0	a[j1] = yr;
1354	0	a[j1 + 1] = yi;
1355	0	a[k1] = xr;
1356	0	a[k1 + 1] = xi;
1357	0	j1 += nm;
1358	0	k1 -= nm;
1359	0	xr = a[j1];
1360	0	xi = -a[j1 + 1];
1361	0	yr = a[k1];
1362	0	yi = -a[k1 + 1];
1363	0	a[j1] = yr;
1364	0	a[j1 + 1] = yi;
1365	0	a[k1] = xr;
1366	0	a[k1 + 1] = xi;
1367	0	j1 += nm;
1368	0	k1 += 2 * nm;
1369	0	xr = a[j1];
1370	0	xi = -a[j1 + 1];
1371	0	yr = a[k1];
1372	0	yi = -a[k1 + 1];
1373	0	a[j1] = yr;
1374	0	a[j1 + 1] = yi;
1375	0	a[k1] = xr;
1376	0	a[k1 + 1] = xi;
1377	0	j1 -= nh;
1378	0	k1 -= 2;
1379	0	xr = a[j1];
1380	0	xi = -a[j1 + 1];
1381	0	yr = a[k1];
1382	0	yi = -a[k1 + 1];
1383	0	a[j1] = yr;
1384	0	a[j1 + 1] = yi;
1385	0	a[k1] = xr;
1386	0	a[k1 + 1] = xi;
1387	0	j1 -= nm;
1388	0	k1 -= 2 * nm;
1389	0	xr = a[j1];
1390	0	xi = -a[j1 + 1];
1391	0	yr = a[k1];
1392	0	yi = -a[k1 + 1];
1393	0	a[j1] = yr;
1394	0	a[j1 + 1] = yi;
1395	0	a[k1] = xr;
1396	0	a[k1 + 1] = xi;
1397	0	j1 -= nm;
1398	0	k1 += nm;
1399	0	xr = a[j1];
1400	0	xi = -a[j1 + 1];
1401	0	yr = a[k1];
1402	0	yi = -a[k1 + 1];
1403	0	a[j1] = yr;
1404	0	a[j1 + 1] = yi;
1405	0	a[k1] = xr;
1406	0	a[k1 + 1] = xi;
1407	0	j1 -= nm;
1408	0	k1 -= 2 * nm;
1409	0	xr = a[j1];
1410	0	xi = -a[j1 + 1];
1411	0	yr = a[k1];
1412	0	yi = -a[k1 + 1];
1413	0	a[j1] = yr;
1414	0	a[j1 + 1] = yi;
1415	0	a[k1] = xr;
1416	0	a[k1 + 1] = xi;
1417	0	}
1418	0	k1 = 4 * k + 2 * ip[m + k];
1419	0	j1 = k1 + 2;
1420	0	k1 += nh;
1421	0	a[j1 - 1] = -a[j1 - 1];
1422	0	xr = a[j1];
1423	0	xi = -a[j1 + 1];
1424	0	yr = a[k1];
1425	0	yi = -a[k1 + 1];
1426	0	a[j1] = yr;
1427	0	a[j1 + 1] = yi;
1428	0	a[k1] = xr;
1429	0	a[k1 + 1] = xi;
1430	0	a[k1 + 3] = -a[k1 + 3];
1431	0	j1 += nm;
1432	0	k1 += 2 * nm;
1433	0	xr = a[j1];
1434	0	xi = -a[j1 + 1];
1435	0	yr = a[k1];
1436	0	yi = -a[k1 + 1];
1437	0	a[j1] = yr;
1438	0	a[j1 + 1] = yi;
1439	0	a[k1] = xr;
1440	0	a[k1 + 1] = xi;
1441	0	j1 += nm;
1442	0	k1 -= nm;
1443	0	xr = a[j1];
1444	0	xi = -a[j1 + 1];
1445	0	yr = a[k1];
1446	0	yi = -a[k1 + 1];
1447	0	a[j1] = yr;
1448	0	a[j1 + 1] = yi;
1449	0	a[k1] = xr;
1450	0	a[k1 + 1] = xi;
1451	0	j1 -= 2;
1452	0	k1 -= nh;
1453	0	xr = a[j1];
1454	0	xi = -a[j1 + 1];
1455	0	yr = a[k1];
1456	0	yi = -a[k1 + 1];
1457	0	a[j1] = yr;
1458	0	a[j1 + 1] = yi;
1459	0	a[k1] = xr;
1460	0	a[k1 + 1] = xi;
1461	0	j1 += nh + 2;
1462	0	k1 += nh + 2;
1463	0	xr = a[j1];
1464	0	xi = -a[j1 + 1];
1465	0	yr = a[k1];
1466	0	yi = -a[k1 + 1];
1467	0	a[j1] = yr;
1468	0	a[j1 + 1] = yi;
1469	0	a[k1] = xr;
1470	0	a[k1 + 1] = xi;
1471	0	j1 -= nh - nm;
1472	0	k1 += 2 * nm - 2;
1473	0	a[j1 - 1] = -a[j1 - 1];
1474	0	xr = a[j1];
1475	0	xi = -a[j1 + 1];
1476	0	yr = a[k1];
1477	0	yi = -a[k1 + 1];
1478	0	a[j1] = yr;
1479	0	a[j1 + 1] = yi;
1480	0	a[k1] = xr;
1481	0	a[k1 + 1] = xi;
1482	0	a[k1 + 3] = -a[k1 + 3];
1483	0	}
1484	0	} else {
1485	0	for (k = 0; k < m; k++) {
1486	0	for (j = 0; j < k; j++) {
1487	0	j1 = 4 * j + ip[m + k];
1488	0	k1 = 4 * k + ip[m + j];
1489	0	xr = a[j1];
1490	0	xi = -a[j1 + 1];
1491	0	yr = a[k1];
1492	0	yi = -a[k1 + 1];
1493	0	a[j1] = yr;
1494	0	a[j1 + 1] = yi;
1495	0	a[k1] = xr;
1496	0	a[k1 + 1] = xi;
1497	0	j1 += nm;
1498	0	k1 += nm;
1499	0	xr = a[j1];
1500	0	xi = -a[j1 + 1];
1501	0	yr = a[k1];
1502	0	yi = -a[k1 + 1];
1503	0	a[j1] = yr;
1504	0	a[j1 + 1] = yi;
1505	0	a[k1] = xr;
1506	0	a[k1 + 1] = xi;
1507	0	j1 += nh;
1508	0	k1 += 2;
1509	0	xr = a[j1];
1510	0	xi = -a[j1 + 1];
1511	0	yr = a[k1];
1512	0	yi = -a[k1 + 1];
1513	0	a[j1] = yr;
1514	0	a[j1 + 1] = yi;
1515	0	a[k1] = xr;
1516	0	a[k1 + 1] = xi;
1517	0	j1 -= nm;
1518	0	k1 -= nm;
1519	0	xr = a[j1];
1520	0	xi = -a[j1 + 1];
1521	0	yr = a[k1];
1522	0	yi = -a[k1 + 1];
1523	0	a[j1] = yr;
1524	0	a[j1 + 1] = yi;
1525	0	a[k1] = xr;
1526	0	a[k1 + 1] = xi;
1527	0	j1 += 2;
1528	0	k1 += nh;
1529	0	xr = a[j1];
1530	0	xi = -a[j1 + 1];
1531	0	yr = a[k1];
1532	0	yi = -a[k1 + 1];
1533	0	a[j1] = yr;
1534	0	a[j1 + 1] = yi;
1535	0	a[k1] = xr;
1536	0	a[k1 + 1] = xi;
1537	0	j1 += nm;
1538	0	k1 += nm;
1539	0	xr = a[j1];
1540	0	xi = -a[j1 + 1];
1541	0	yr = a[k1];
1542	0	yi = -a[k1 + 1];
1543	0	a[j1] = yr;
1544	0	a[j1 + 1] = yi;
1545	0	a[k1] = xr;
1546	0	a[k1 + 1] = xi;
1547	0	j1 -= nh;
1548	0	k1 -= 2;
1549	0	xr = a[j1];
1550	0	xi = -a[j1 + 1];
1551	0	yr = a[k1];
1552	0	yi = -a[k1 + 1];
1553	0	a[j1] = yr;
1554	0	a[j1 + 1] = yi;
1555	0	a[k1] = xr;
1556	0	a[k1 + 1] = xi;
1557	0	j1 -= nm;
1558	0	k1 -= nm;
1559	0	xr = a[j1];
1560	0	xi = -a[j1 + 1];
1561	0	yr = a[k1];
1562	0	yi = -a[k1 + 1];
1563	0	a[j1] = yr;
1564	0	a[j1 + 1] = yi;
1565	0	a[k1] = xr;
1566	0	a[k1 + 1] = xi;
1567	0	}
1568	0	k1 = 4 * k + ip[m + k];
1569	0	j1 = k1 + 2;
1570	0	k1 += nh;
1571	0	a[j1 - 1] = -a[j1 - 1];
1572	0	xr = a[j1];
1573	0	xi = -a[j1 + 1];
1574	0	yr = a[k1];
1575	0	yi = -a[k1 + 1];
1576	0	a[j1] = yr;
1577	0	a[j1 + 1] = yi;
1578	0	a[k1] = xr;
1579	0	a[k1 + 1] = xi;
1580	0	a[k1 + 3] = -a[k1 + 3];
1581	0	j1 += nm;
1582	0	k1 += nm;
1583	0	a[j1 - 1] = -a[j1 - 1];
1584	0	xr = a[j1];
1585	0	xi = -a[j1 + 1];
1586	0	yr = a[k1];
1587	0	yi = -a[k1 + 1];
1588	0	a[j1] = yr;
1589	0	a[j1 + 1] = yi;
1590	0	a[k1] = xr;
1591	0	a[k1 + 1] = xi;
1592	0	a[k1 + 3] = -a[k1 + 3];
1593	0	}
1594	0	}
1595	0	}
1596
1597
1598		void bitrv216(double *a)
1599	0	{
1600	0	double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1601	0	x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i,
1602	0	x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i;
1603
1604	0	x1r = a[2];
1605	0	x1i = a[3];
1606	0	x2r = a[4];
1607	0	x2i = a[5];
1608	0	x3r = a[6];
1609	0	x3i = a[7];
1610	0	x4r = a[8];
1611	0	x4i = a[9];
1612	0	x5r = a[10];
1613	0	x5i = a[11];
1614	0	x7r = a[14];
1615	0	x7i = a[15];
1616	0	x8r = a[16];
1617	0	x8i = a[17];
1618	0	x10r = a[20];
1619	0	x10i = a[21];
1620	0	x11r = a[22];
1621	0	x11i = a[23];
1622	0	x12r = a[24];
1623	0	x12i = a[25];
1624	0	x13r = a[26];
1625	0	x13i = a[27];
1626	0	x14r = a[28];
1627	0	x14i = a[29];
1628	0	a[2] = x8r;
1629	0	a[3] = x8i;
1630	0	a[4] = x4r;
1631	0	a[5] = x4i;
1632	0	a[6] = x12r;
1633	0	a[7] = x12i;
1634	0	a[8] = x2r;
1635	0	a[9] = x2i;
1636	0	a[10] = x10r;
1637	0	a[11] = x10i;
1638	0	a[14] = x14r;
1639	0	a[15] = x14i;
1640	0	a[16] = x1r;
1641	0	a[17] = x1i;
1642	0	a[20] = x5r;
1643	0	a[21] = x5i;
1644	0	a[22] = x13r;
1645	0	a[23] = x13i;
1646	0	a[24] = x3r;
1647	0	a[25] = x3i;
1648	0	a[26] = x11r;
1649	0	a[27] = x11i;
1650	0	a[28] = x7r;
1651	0	a[29] = x7i;
1652	0	}
1653
1654
1655		void bitrv216neg(double *a)
1656	0	{
1657	0	double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1658	0	x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i,
1659	0	x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i,
1660	0	x13r, x13i, x14r, x14i, x15r, x15i;
1661
1662	0	x1r = a[2];
1663	0	x1i = a[3];
1664	0	x2r = a[4];
1665	0	x2i = a[5];
1666	0	x3r = a[6];
1667	0	x3i = a[7];
1668	0	x4r = a[8];
1669	0	x4i = a[9];
1670	0	x5r = a[10];
1671	0	x5i = a[11];
1672	0	x6r = a[12];
1673	0	x6i = a[13];
1674	0	x7r = a[14];
1675	0	x7i = a[15];
1676	0	x8r = a[16];
1677	0	x8i = a[17];
1678	0	x9r = a[18];
1679	0	x9i = a[19];
1680	0	x10r = a[20];
1681	0	x10i = a[21];
1682	0	x11r = a[22];
1683	0	x11i = a[23];
1684	0	x12r = a[24];
1685	0	x12i = a[25];
1686	0	x13r = a[26];
1687	0	x13i = a[27];
1688	0	x14r = a[28];
1689	0	x14i = a[29];
1690	0	x15r = a[30];
1691	0	x15i = a[31];
1692	0	a[2] = x15r;
1693	0	a[3] = x15i;
1694	0	a[4] = x7r;
1695	0	a[5] = x7i;
1696	0	a[6] = x11r;
1697	0	a[7] = x11i;
1698	0	a[8] = x3r;
1699	0	a[9] = x3i;
1700	0	a[10] = x13r;
1701	0	a[11] = x13i;
1702	0	a[12] = x5r;
1703	0	a[13] = x5i;
1704	0	a[14] = x9r;
1705	0	a[15] = x9i;
1706	0	a[16] = x1r;
1707	0	a[17] = x1i;
1708	0	a[18] = x14r;
1709	0	a[19] = x14i;
1710	0	a[20] = x6r;
1711	0	a[21] = x6i;
1712	0	a[22] = x10r;
1713	0	a[23] = x10i;
1714	0	a[24] = x2r;
1715	0	a[25] = x2i;
1716	0	a[26] = x12r;
1717	0	a[27] = x12i;
1718	0	a[28] = x4r;
1719	0	a[29] = x4i;
1720	0	a[30] = x8r;
1721	0	a[31] = x8i;
1722	0	}
1723
1724
1725		void bitrv208(double *a)
1726	0	{
1727	0	double x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
1728
1729	0	x1r = a[2];
1730	0	x1i = a[3];
1731	0	x3r = a[6];
1732	0	x3i = a[7];
1733	0	x4r = a[8];
1734	0	x4i = a[9];
1735	0	x6r = a[12];
1736	0	x6i = a[13];
1737	0	a[2] = x4r;
1738	0	a[3] = x4i;
1739	0	a[6] = x6r;
1740	0	a[7] = x6i;
1741	0	a[8] = x1r;
1742	0	a[9] = x1i;
1743	0	a[12] = x3r;
1744	0	a[13] = x3i;
1745	0	}
1746
1747
1748		void bitrv208neg(double *a)
1749	0	{
1750	0	double x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i,
1751	0	x5r, x5i, x6r, x6i, x7r, x7i;
1752
1753	0	x1r = a[2];
1754	0	x1i = a[3];
1755	0	x2r = a[4];
1756	0	x2i = a[5];
1757	0	x3r = a[6];
1758	0	x3i = a[7];
1759	0	x4r = a[8];
1760	0	x4i = a[9];
1761	0	x5r = a[10];
1762	0	x5i = a[11];
1763	0	x6r = a[12];
1764	0	x6i = a[13];
1765	0	x7r = a[14];
1766	0	x7i = a[15];
1767	0	a[2] = x7r;
1768	0	a[3] = x7i;
1769	0	a[4] = x3r;
1770	0	a[5] = x3i;
1771	0	a[6] = x5r;
1772	0	a[7] = x5i;
1773	0	a[8] = x1r;
1774	0	a[9] = x1i;
1775	0	a[10] = x6r;
1776	0	a[11] = x6i;
1777	0	a[12] = x2r;
1778	0	a[13] = x2i;
1779	0	a[14] = x4r;
1780	0	a[15] = x4i;
1781	0	}
1782
1783
1784		void cftf1st(int n, double a, double w)
1785	0	{
1786	0	int j, j0, j1, j2, j3, k, m, mh;
1787	0	double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
1788	0	wd1r, wd1i, wd3r, wd3i;
1789	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1790	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1791
1792	0	mh = n >> 3;
1793	0	m = 2 * mh;
1794	0	j1 = m;
1795	0	j2 = j1 + m;
1796	0	j3 = j2 + m;
1797	0	x0r = a[0] + a[j2];
1798	0	x0i = a[1] + a[j2 + 1];
1799	0	x1r = a[0] - a[j2];
1800	0	x1i = a[1] - a[j2 + 1];
1801	0	x2r = a[j1] + a[j3];
1802	0	x2i = a[j1 + 1] + a[j3 + 1];
1803	0	x3r = a[j1] - a[j3];
1804	0	x3i = a[j1 + 1] - a[j3 + 1];
1805	0	a[0] = x0r + x2r;
1806	0	a[1] = x0i + x2i;
1807	0	a[j1] = x0r - x2r;
1808	0	a[j1 + 1] = x0i - x2i;
1809	0	a[j2] = x1r - x3i;
1810	0	a[j2 + 1] = x1i + x3r;
1811	0	a[j3] = x1r + x3i;
1812	0	a[j3 + 1] = x1i - x3r;
1813	0	wn4r = w[1];
1814	0	csc1 = w[2];
1815	0	csc3 = w[3];
1816	0	wd1r = 1;
1817	0	wd1i = 0;
1818	0	wd3r = 1;
1819	0	wd3i = 0;
1820	0	k = 0;
1821	0	for (j = 2; j < mh - 2; j += 4) {
1822	0	k += 4;
1823	0	wk1r = csc1 * (wd1r + w[k]);
1824	0	wk1i = csc1 * (wd1i + w[k + 1]);
1825	0	wk3r = csc3 * (wd3r + w[k + 2]);
1826	0	wk3i = csc3 * (wd3i + w[k + 3]);
1827	0	wd1r = w[k];
1828	0	wd1i = w[k + 1];
1829	0	wd3r = w[k + 2];
1830	0	wd3i = w[k + 3];
1831	0	j1 = j + m;
1832	0	j2 = j1 + m;
1833	0	j3 = j2 + m;
1834	0	x0r = a[j] + a[j2];
1835	0	x0i = a[j + 1] + a[j2 + 1];
1836	0	x1r = a[j] - a[j2];
1837	0	x1i = a[j + 1] - a[j2 + 1];
1838	0	y0r = a[j + 2] + a[j2 + 2];
1839	0	y0i = a[j + 3] + a[j2 + 3];
1840	0	y1r = a[j + 2] - a[j2 + 2];
1841	0	y1i = a[j + 3] - a[j2 + 3];
1842	0	x2r = a[j1] + a[j3];
1843	0	x2i = a[j1 + 1] + a[j3 + 1];
1844	0	x3r = a[j1] - a[j3];
1845	0	x3i = a[j1 + 1] - a[j3 + 1];
1846	0	y2r = a[j1 + 2] + a[j3 + 2];
1847	0	y2i = a[j1 + 3] + a[j3 + 3];
1848	0	y3r = a[j1 + 2] - a[j3 + 2];
1849	0	y3i = a[j1 + 3] - a[j3 + 3];
1850	0	a[j] = x0r + x2r;
1851	0	a[j + 1] = x0i + x2i;
1852	0	a[j + 2] = y0r + y2r;
1853	0	a[j + 3] = y0i + y2i;
1854	0	a[j1] = x0r - x2r;
1855	0	a[j1 + 1] = x0i - x2i;
1856	0	a[j1 + 2] = y0r - y2r;
1857	0	a[j1 + 3] = y0i - y2i;
1858	0	x0r = x1r - x3i;
1859	0	x0i = x1i + x3r;
1860	0	a[j2] = wk1r * x0r - wk1i * x0i;
1861	0	a[j2 + 1] = wk1r * x0i + wk1i * x0r;
1862	0	x0r = y1r - y3i;
1863	0	x0i = y1i + y3r;
1864	0	a[j2 + 2] = wd1r * x0r - wd1i * x0i;
1865	0	a[j2 + 3] = wd1r * x0i + wd1i * x0r;
1866	0	x0r = x1r + x3i;
1867	0	x0i = x1i - x3r;
1868	0	a[j3] = wk3r * x0r + wk3i * x0i;
1869	0	a[j3 + 1] = wk3r * x0i - wk3i * x0r;
1870	0	x0r = y1r + y3i;
1871	0	x0i = y1i - y3r;
1872	0	a[j3 + 2] = wd3r * x0r + wd3i * x0i;
1873	0	a[j3 + 3] = wd3r * x0i - wd3i * x0r;
1874	0	j0 = m - j;
1875	0	j1 = j0 + m;
1876	0	j2 = j1 + m;
1877	0	j3 = j2 + m;
1878	0	x0r = a[j0] + a[j2];
1879	0	x0i = a[j0 + 1] + a[j2 + 1];
1880	0	x1r = a[j0] - a[j2];
1881	0	x1i = a[j0 + 1] - a[j2 + 1];
1882	0	y0r = a[j0 - 2] + a[j2 - 2];
1883	0	y0i = a[j0 - 1] + a[j2 - 1];
1884	0	y1r = a[j0 - 2] - a[j2 - 2];
1885	0	y1i = a[j0 - 1] - a[j2 - 1];
1886	0	x2r = a[j1] + a[j3];
1887	0	x2i = a[j1 + 1] + a[j3 + 1];
1888	0	x3r = a[j1] - a[j3];
1889	0	x3i = a[j1 + 1] - a[j3 + 1];
1890	0	y2r = a[j1 - 2] + a[j3 - 2];
1891	0	y2i = a[j1 - 1] + a[j3 - 1];
1892	0	y3r = a[j1 - 2] - a[j3 - 2];
1893	0	y3i = a[j1 - 1] - a[j3 - 1];
1894	0	a[j0] = x0r + x2r;
1895	0	a[j0 + 1] = x0i + x2i;
1896	0	a[j0 - 2] = y0r + y2r;
1897	0	a[j0 - 1] = y0i + y2i;
1898	0	a[j1] = x0r - x2r;
1899	0	a[j1 + 1] = x0i - x2i;
1900	0	a[j1 - 2] = y0r - y2r;
1901	0	a[j1 - 1] = y0i - y2i;
1902	0	x0r = x1r - x3i;
1903	0	x0i = x1i + x3r;
1904	0	a[j2] = wk1i * x0r - wk1r * x0i;
1905	0	a[j2 + 1] = wk1i * x0i + wk1r * x0r;
1906	0	x0r = y1r - y3i;
1907	0	x0i = y1i + y3r;
1908	0	a[j2 - 2] = wd1i * x0r - wd1r * x0i;
1909	0	a[j2 - 1] = wd1i * x0i + wd1r * x0r;
1910	0	x0r = x1r + x3i;
1911	0	x0i = x1i - x3r;
1912	0	a[j3] = wk3i * x0r + wk3r * x0i;
1913	0	a[j3 + 1] = wk3i * x0i - wk3r * x0r;
1914	0	x0r = y1r + y3i;
1915	0	x0i = y1i - y3r;
1916	0	a[j3 - 2] = wd3i * x0r + wd3r * x0i;
1917	0	a[j3 - 1] = wd3i * x0i - wd3r * x0r;
1918	0	}
1919	0	wk1r = csc1 * (wd1r + wn4r);
1920	0	wk1i = csc1 * (wd1i + wn4r);
1921	0	wk3r = csc3 * (wd3r - wn4r);
1922	0	wk3i = csc3 * (wd3i - wn4r);
1923	0	j0 = mh;
1924	0	j1 = j0 + m;
1925	0	j2 = j1 + m;
1926	0	j3 = j2 + m;
1927	0	x0r = a[j0 - 2] + a[j2 - 2];
1928	0	x0i = a[j0 - 1] + a[j2 - 1];
1929	0	x1r = a[j0 - 2] - a[j2 - 2];
1930	0	x1i = a[j0 - 1] - a[j2 - 1];
1931	0	x2r = a[j1 - 2] + a[j3 - 2];
1932	0	x2i = a[j1 - 1] + a[j3 - 1];
1933	0	x3r = a[j1 - 2] - a[j3 - 2];
1934	0	x3i = a[j1 - 1] - a[j3 - 1];
1935	0	a[j0 - 2] = x0r + x2r;
1936	0	a[j0 - 1] = x0i + x2i;
1937	0	a[j1 - 2] = x0r - x2r;
1938	0	a[j1 - 1] = x0i - x2i;
1939	0	x0r = x1r - x3i;
1940	0	x0i = x1i + x3r;
1941	0	a[j2 - 2] = wk1r * x0r - wk1i * x0i;
1942	0	a[j2 - 1] = wk1r * x0i + wk1i * x0r;
1943	0	x0r = x1r + x3i;
1944	0	x0i = x1i - x3r;
1945	0	a[j3 - 2] = wk3r * x0r + wk3i * x0i;
1946	0	a[j3 - 1] = wk3r * x0i - wk3i * x0r;
1947	0	x0r = a[j0] + a[j2];
1948	0	x0i = a[j0 + 1] + a[j2 + 1];
1949	0	x1r = a[j0] - a[j2];
1950	0	x1i = a[j0 + 1] - a[j2 + 1];
1951	0	x2r = a[j1] + a[j3];
1952	0	x2i = a[j1 + 1] + a[j3 + 1];
1953	0	x3r = a[j1] - a[j3];
1954	0	x3i = a[j1 + 1] - a[j3 + 1];
1955	0	a[j0] = x0r + x2r;
1956	0	a[j0 + 1] = x0i + x2i;
1957	0	a[j1] = x0r - x2r;
1958	0	a[j1 + 1] = x0i - x2i;
1959	0	x0r = x1r - x3i;
1960	0	x0i = x1i + x3r;
1961	0	a[j2] = wn4r * (x0r - x0i);
1962	0	a[j2 + 1] = wn4r * (x0i + x0r);
1963	0	x0r = x1r + x3i;
1964	0	x0i = x1i - x3r;
1965	0	a[j3] = -wn4r * (x0r + x0i);
1966	0	a[j3 + 1] = -wn4r * (x0i - x0r);
1967	0	x0r = a[j0 + 2] + a[j2 + 2];
1968	0	x0i = a[j0 + 3] + a[j2 + 3];
1969	0	x1r = a[j0 + 2] - a[j2 + 2];
1970	0	x1i = a[j0 + 3] - a[j2 + 3];
1971	0	x2r = a[j1 + 2] + a[j3 + 2];
1972	0	x2i = a[j1 + 3] + a[j3 + 3];
1973	0	x3r = a[j1 + 2] - a[j3 + 2];
1974	0	x3i = a[j1 + 3] - a[j3 + 3];
1975	0	a[j0 + 2] = x0r + x2r;
1976	0	a[j0 + 3] = x0i + x2i;
1977	0	a[j1 + 2] = x0r - x2r;
1978	0	a[j1 + 3] = x0i - x2i;
1979	0	x0r = x1r - x3i;
1980	0	x0i = x1i + x3r;
1981	0	a[j2 + 2] = wk1i * x0r - wk1r * x0i;
1982	0	a[j2 + 3] = wk1i * x0i + wk1r * x0r;
1983	0	x0r = x1r + x3i;
1984	0	x0i = x1i - x3r;
1985	0	a[j3 + 2] = wk3i * x0r + wk3r * x0i;
1986	0	a[j3 + 3] = wk3i * x0i - wk3r * x0r;
1987	0	}
1988
1989
1990		void cftb1st(int n, double a, double w)
1991	0	{
1992	0	int j, j0, j1, j2, j3, k, m, mh;
1993	0	double wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i,
1994	0	wd1r, wd1i, wd3r, wd3i;
1995	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
1996	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i;
1997
1998	0	mh = n >> 3;
1999	0	m = 2 * mh;
2000	0	j1 = m;
2001	0	j2 = j1 + m;
2002	0	j3 = j2 + m;
2003	0	x0r = a[0] + a[j2];
2004	0	x0i = -a[1] - a[j2 + 1];
2005	0	x1r = a[0] - a[j2];
2006	0	x1i = -a[1] + a[j2 + 1];
2007	0	x2r = a[j1] + a[j3];
2008	0	x2i = a[j1 + 1] + a[j3 + 1];
2009	0	x3r = a[j1] - a[j3];
2010	0	x3i = a[j1 + 1] - a[j3 + 1];
2011	0	a[0] = x0r + x2r;
2012	0	a[1] = x0i - x2i;
2013	0	a[j1] = x0r - x2r;
2014	0	a[j1 + 1] = x0i + x2i;
2015	0	a[j2] = x1r + x3i;
2016	0	a[j2 + 1] = x1i + x3r;
2017	0	a[j3] = x1r - x3i;
2018	0	a[j3 + 1] = x1i - x3r;
2019	0	wn4r = w[1];
2020	0	csc1 = w[2];
2021	0	csc3 = w[3];
2022	0	wd1r = 1;
2023	0	wd1i = 0;
2024	0	wd3r = 1;
2025	0	wd3i = 0;
2026	0	k = 0;
2027	0	for (j = 2; j < mh - 2; j += 4) {
2028	0	k += 4;
2029	0	wk1r = csc1 * (wd1r + w[k]);
2030	0	wk1i = csc1 * (wd1i + w[k + 1]);
2031	0	wk3r = csc3 * (wd3r + w[k + 2]);
2032	0	wk3i = csc3 * (wd3i + w[k + 3]);
2033	0	wd1r = w[k];
2034	0	wd1i = w[k + 1];
2035	0	wd3r = w[k + 2];
2036	0	wd3i = w[k + 3];
2037	0	j1 = j + m;
2038	0	j2 = j1 + m;
2039	0	j3 = j2 + m;
2040	0	x0r = a[j] + a[j2];
2041	0	x0i = -a[j + 1] - a[j2 + 1];
2042	0	x1r = a[j] - a[j2];
2043	0	x1i = -a[j + 1] + a[j2 + 1];
2044	0	y0r = a[j + 2] + a[j2 + 2];
2045	0	y0i = -a[j + 3] - a[j2 + 3];
2046	0	y1r = a[j + 2] - a[j2 + 2];
2047	0	y1i = -a[j + 3] + a[j2 + 3];
2048	0	x2r = a[j1] + a[j3];
2049	0	x2i = a[j1 + 1] + a[j3 + 1];
2050	0	x3r = a[j1] - a[j3];
2051	0	x3i = a[j1 + 1] - a[j3 + 1];
2052	0	y2r = a[j1 + 2] + a[j3 + 2];
2053	0	y2i = a[j1 + 3] + a[j3 + 3];
2054	0	y3r = a[j1 + 2] - a[j3 + 2];
2055	0	y3i = a[j1 + 3] - a[j3 + 3];
2056	0	a[j] = x0r + x2r;
2057	0	a[j + 1] = x0i - x2i;
2058	0	a[j + 2] = y0r + y2r;
2059	0	a[j + 3] = y0i - y2i;
2060	0	a[j1] = x0r - x2r;
2061	0	a[j1 + 1] = x0i + x2i;
2062	0	a[j1 + 2] = y0r - y2r;
2063	0	a[j1 + 3] = y0i + y2i;
2064	0	x0r = x1r + x3i;
2065	0	x0i = x1i + x3r;
2066	0	a[j2] = wk1r * x0r - wk1i * x0i;
2067	0	a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2068	0	x0r = y1r + y3i;
2069	0	x0i = y1i + y3r;
2070	0	a[j2 + 2] = wd1r * x0r - wd1i * x0i;
2071	0	a[j2 + 3] = wd1r * x0i + wd1i * x0r;
2072	0	x0r = x1r - x3i;
2073	0	x0i = x1i - x3r;
2074	0	a[j3] = wk3r * x0r + wk3i * x0i;
2075	0	a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2076	0	x0r = y1r - y3i;
2077	0	x0i = y1i - y3r;
2078	0	a[j3 + 2] = wd3r * x0r + wd3i * x0i;
2079	0	a[j3 + 3] = wd3r * x0i - wd3i * x0r;
2080	0	j0 = m - j;
2081	0	j1 = j0 + m;
2082	0	j2 = j1 + m;
2083	0	j3 = j2 + m;
2084	0	x0r = a[j0] + a[j2];
2085	0	x0i = -a[j0 + 1] - a[j2 + 1];
2086	0	x1r = a[j0] - a[j2];
2087	0	x1i = -a[j0 + 1] + a[j2 + 1];
2088	0	y0r = a[j0 - 2] + a[j2 - 2];
2089	0	y0i = -a[j0 - 1] - a[j2 - 1];
2090	0	y1r = a[j0 - 2] - a[j2 - 2];
2091	0	y1i = -a[j0 - 1] + a[j2 - 1];
2092	0	x2r = a[j1] + a[j3];
2093	0	x2i = a[j1 + 1] + a[j3 + 1];
2094	0	x3r = a[j1] - a[j3];
2095	0	x3i = a[j1 + 1] - a[j3 + 1];
2096	0	y2r = a[j1 - 2] + a[j3 - 2];
2097	0	y2i = a[j1 - 1] + a[j3 - 1];
2098	0	y3r = a[j1 - 2] - a[j3 - 2];
2099	0	y3i = a[j1 - 1] - a[j3 - 1];
2100	0	a[j0] = x0r + x2r;
2101	0	a[j0 + 1] = x0i - x2i;
2102	0	a[j0 - 2] = y0r + y2r;
2103	0	a[j0 - 1] = y0i - y2i;
2104	0	a[j1] = x0r - x2r;
2105	0	a[j1 + 1] = x0i + x2i;
2106	0	a[j1 - 2] = y0r - y2r;
2107	0	a[j1 - 1] = y0i + y2i;
2108	0	x0r = x1r + x3i;
2109	0	x0i = x1i + x3r;
2110	0	a[j2] = wk1i * x0r - wk1r * x0i;
2111	0	a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2112	0	x0r = y1r + y3i;
2113	0	x0i = y1i + y3r;
2114	0	a[j2 - 2] = wd1i * x0r - wd1r * x0i;
2115	0	a[j2 - 1] = wd1i * x0i + wd1r * x0r;
2116	0	x0r = x1r - x3i;
2117	0	x0i = x1i - x3r;
2118	0	a[j3] = wk3i * x0r + wk3r * x0i;
2119	0	a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2120	0	x0r = y1r - y3i;
2121	0	x0i = y1i - y3r;
2122	0	a[j3 - 2] = wd3i * x0r + wd3r * x0i;
2123	0	a[j3 - 1] = wd3i * x0i - wd3r * x0r;
2124	0	}
2125	0	wk1r = csc1 * (wd1r + wn4r);
2126	0	wk1i = csc1 * (wd1i + wn4r);
2127	0	wk3r = csc3 * (wd3r - wn4r);
2128	0	wk3i = csc3 * (wd3i - wn4r);
2129	0	j0 = mh;
2130	0	j1 = j0 + m;
2131	0	j2 = j1 + m;
2132	0	j3 = j2 + m;
2133	0	x0r = a[j0 - 2] + a[j2 - 2];
2134	0	x0i = -a[j0 - 1] - a[j2 - 1];
2135	0	x1r = a[j0 - 2] - a[j2 - 2];
2136	0	x1i = -a[j0 - 1] + a[j2 - 1];
2137	0	x2r = a[j1 - 2] + a[j3 - 2];
2138	0	x2i = a[j1 - 1] + a[j3 - 1];
2139	0	x3r = a[j1 - 2] - a[j3 - 2];
2140	0	x3i = a[j1 - 1] - a[j3 - 1];
2141	0	a[j0 - 2] = x0r + x2r;
2142	0	a[j0 - 1] = x0i - x2i;
2143	0	a[j1 - 2] = x0r - x2r;
2144	0	a[j1 - 1] = x0i + x2i;
2145	0	x0r = x1r + x3i;
2146	0	x0i = x1i + x3r;
2147	0	a[j2 - 2] = wk1r * x0r - wk1i * x0i;
2148	0	a[j2 - 1] = wk1r * x0i + wk1i * x0r;
2149	0	x0r = x1r - x3i;
2150	0	x0i = x1i - x3r;
2151	0	a[j3 - 2] = wk3r * x0r + wk3i * x0i;
2152	0	a[j3 - 1] = wk3r * x0i - wk3i * x0r;
2153	0	x0r = a[j0] + a[j2];
2154	0	x0i = -a[j0 + 1] - a[j2 + 1];
2155	0	x1r = a[j0] - a[j2];
2156	0	x1i = -a[j0 + 1] + a[j2 + 1];
2157	0	x2r = a[j1] + a[j3];
2158	0	x2i = a[j1 + 1] + a[j3 + 1];
2159	0	x3r = a[j1] - a[j3];
2160	0	x3i = a[j1 + 1] - a[j3 + 1];
2161	0	a[j0] = x0r + x2r;
2162	0	a[j0 + 1] = x0i - x2i;
2163	0	a[j1] = x0r - x2r;
2164	0	a[j1 + 1] = x0i + x2i;
2165	0	x0r = x1r + x3i;
2166	0	x0i = x1i + x3r;
2167	0	a[j2] = wn4r * (x0r - x0i);
2168	0	a[j2 + 1] = wn4r * (x0i + x0r);
2169	0	x0r = x1r - x3i;
2170	0	x0i = x1i - x3r;
2171	0	a[j3] = -wn4r * (x0r + x0i);
2172	0	a[j3 + 1] = -wn4r * (x0i - x0r);
2173	0	x0r = a[j0 + 2] + a[j2 + 2];
2174	0	x0i = -a[j0 + 3] - a[j2 + 3];
2175	0	x1r = a[j0 + 2] - a[j2 + 2];
2176	0	x1i = -a[j0 + 3] + a[j2 + 3];
2177	0	x2r = a[j1 + 2] + a[j3 + 2];
2178	0	x2i = a[j1 + 3] + a[j3 + 3];
2179	0	x3r = a[j1 + 2] - a[j3 + 2];
2180	0	x3i = a[j1 + 3] - a[j3 + 3];
2181	0	a[j0 + 2] = x0r + x2r;
2182	0	a[j0 + 3] = x0i - x2i;
2183	0	a[j1 + 2] = x0r - x2r;
2184	0	a[j1 + 3] = x0i + x2i;
2185	0	x0r = x1r + x3i;
2186	0	x0i = x1i + x3r;
2187	0	a[j2 + 2] = wk1i * x0r - wk1r * x0i;
2188	0	a[j2 + 3] = wk1i * x0i + wk1r * x0r;
2189	0	x0r = x1r - x3i;
2190	0	x0i = x1i - x3r;
2191	0	a[j3 + 2] = wk3i * x0r + wk3r * x0i;
2192	0	a[j3 + 3] = wk3i * x0i - wk3r * x0r;
2193	0	}
2194
2195
2196		#ifdef USE_CDFT_THREADS
2197		struct cdft_arg_st {
2198		int n0;
2199		int n;
2200		double *a;
2201		int nw;
2202		double *w;
2203		};
2204		typedef struct cdft_arg_st cdft_arg_t;
2205
2206
2207		void cftrec4_th(int n, double a, int nw, double w)
2208		{
2209		void cftrec1_th(void p);
2210		void cftrec2_th(void p);
2211		int i, idiv4, m, nthread;
2212		cdft_thread_t th[4];
2213		cdft_arg_t ag[4];
2214
2215		nthread = 2;
2216		idiv4 = 0;
2217		m = n >> 1;
2218		if (n > CDFT_4THREADS_BEGIN_N) {
2219		nthread = 4;
2220		idiv4 = 1;
2221		m >>= 1;
2222		}
2223		for (i = 0; i < nthread; i++) {
2224		ag[i].n0 = n;
2225		ag[i].n = m;
2226		ag[i].a = &a[i * m];
2227		ag[i].nw = nw;
2228		ag[i].w = w;
2229		if (i != idiv4) {
2230		cdft_thread_create(&th[i], cftrec1_th, &ag[i]);
2231		} else {
2232		cdft_thread_create(&th[i], cftrec2_th, &ag[i]);
2233		}
2234		}
2235		for (i = 0; i < nthread; i++) {
2236		cdft_thread_wait(th[i]);
2237		}
2238		}
2239
2240
2241		void cftrec1_th(void p)
2242		{
2243		int cfttree(int n, int j, int k, double a, int nw, double w);
2244		void cftleaf(int n, int isplt, double a, int nw, double w);
2245		void cftmdl1(int n, double a, double w);
2246		int isplt, j, k, m, n, n0, nw;
2247		double a, w;
2248
2249		n0 = ((cdft_arg_t *) p)->n0;
2250		n = ((cdft_arg_t *) p)->n;
2251		a = ((cdft_arg_t *) p)->a;
2252		nw = ((cdft_arg_t *) p)->nw;
2253		w = ((cdft_arg_t *) p)->w;
2254		m = n0;
2255		while (m > 512) {
2256		m >>= 2;
2257		cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2258		}
2259		cftleaf(m, 1, &a[n - m], nw, w);
2260		k = 0;
2261		for (j = n - m; j > 0; j -= m) {
2262		k++;
2263		isplt = cfttree(m, j, k, a, nw, w);
2264		cftleaf(m, isplt, &a[j - m], nw, w);
2265		}
2266		return (void *) 0;
2267		}
2268
2269
2270		void cftrec2_th(void p)
2271		{
2272		int cfttree(int n, int j, int k, double a, int nw, double w);
2273		void cftleaf(int n, int isplt, double a, int nw, double w);
2274		void cftmdl2(int n, double a, double w);
2275		int isplt, j, k, m, n, n0, nw;
2276		double a, w;
2277
2278		n0 = ((cdft_arg_t *) p)->n0;
2279		n = ((cdft_arg_t *) p)->n;
2280		a = ((cdft_arg_t *) p)->a;
2281		nw = ((cdft_arg_t *) p)->nw;
2282		w = ((cdft_arg_t *) p)->w;
2283		k = 1;
2284		m = n0;
2285		while (m > 512) {
2286		m >>= 2;
2287		k <<= 2;
2288		cftmdl2(m, &a[n - m], &w[nw - m]);
2289		}
2290		cftleaf(m, 0, &a[n - m], nw, w);
2291		k >>= 1;
2292		for (j = n - m; j > 0; j -= m) {
2293		k++;
2294		isplt = cfttree(m, j, k, a, nw, w);
2295		cftleaf(m, isplt, &a[j - m], nw, w);
2296		}
2297		return (void *) 0;
2298		}
2299		#endif /* USE_CDFT_THREADS */
2300
2301
2302		void cftrec4(int n, double a, int nw, double w)
2303	0	{
2304	0	int cfttree(int n, int j, int k, double a, int nw, double w);
2305	0	void cftleaf(int n, int isplt, double a, int nw, double w);
2306	0	void cftmdl1(int n, double a, double w);
2307	0	int isplt, j, k, m;
2308
2309	0	m = n;
2310	0	while (m > 512) {
2311	0	m >>= 2;
2312	0	cftmdl1(m, &a[n - m], &w[nw - (m >> 1)]);
2313	0	}
2314	0	cftleaf(m, 1, &a[n - m], nw, w);
2315	0	k = 0;
2316	0	for (j = n - m; j > 0; j -= m) {
2317	0	k++;
2318	0	isplt = cfttree(m, j, k, a, nw, w);
2319	0	cftleaf(m, isplt, &a[j - m], nw, w);
2320	0	}
2321	0	}
2322
2323
2324		int cfttree(int n, int j, int k, double a, int nw, double w)
2325	0	{
2326	0	void cftmdl1(int n, double a, double w);
2327	0	void cftmdl2(int n, double a, double w);
2328	0	int i, isplt, m;
2329
2330	0	if ((k & 3) != 0) {
2331	0	isplt = k & 1;
2332	0	if (isplt != 0) {
2333	0	cftmdl1(n, &a[j - n], &w[nw - (n >> 1)]);
2334	0	} else {
2335	0	cftmdl2(n, &a[j - n], &w[nw - n]);
2336	0	}
2337	0	} else {
2338	0	m = n;
2339	0	for (i = k; (i & 3) == 0; i >>= 2) {
2340	0	m <<= 2;
2341	0	}
2342	0	isplt = i & 1;
2343	0	if (isplt != 0) {
2344	0	while (m > 128) {
2345	0	cftmdl1(m, &a[j - m], &w[nw - (m >> 1)]);
2346	0	m >>= 2;
2347	0	}
2348	0	} else {
2349	0	while (m > 128) {
2350	0	cftmdl2(m, &a[j - m], &w[nw - m]);
2351	0	m >>= 2;
2352	0	}
2353	0	}
2354	0	}
2355	0	return isplt;
2356	0	}
2357
2358
2359		void cftleaf(int n, int isplt, double a, int nw, double w)
2360	0	{
2361	0	void cftmdl1(int n, double a, double w);
2362	0	void cftmdl2(int n, double a, double w);
2363	0	void cftf161(double a, double w);
2364	0	void cftf162(double a, double w);
2365	0	void cftf081(double a, double w);
2366	0	void cftf082(double a, double w);
2367
2368	0	if (n == 512) {
2369	0	cftmdl1(128, a, &w[nw - 64]);
2370	0	cftf161(a, &w[nw - 8]);
2371	0	cftf162(&a[32], &w[nw - 32]);
2372	0	cftf161(&a[64], &w[nw - 8]);
2373	0	cftf161(&a[96], &w[nw - 8]);
2374	0	cftmdl2(128, &a[128], &w[nw - 128]);
2375	0	cftf161(&a[128], &w[nw - 8]);
2376	0	cftf162(&a[160], &w[nw - 32]);
2377	0	cftf161(&a[192], &w[nw - 8]);
2378	0	cftf162(&a[224], &w[nw - 32]);
2379	0	cftmdl1(128, &a[256], &w[nw - 64]);
2380	0	cftf161(&a[256], &w[nw - 8]);
2381	0	cftf162(&a[288], &w[nw - 32]);
2382	0	cftf161(&a[320], &w[nw - 8]);
2383	0	cftf161(&a[352], &w[nw - 8]);
2384	0	if (isplt != 0) {
2385	0	cftmdl1(128, &a[384], &w[nw - 64]);
2386	0	cftf161(&a[480], &w[nw - 8]);
2387	0	} else {
2388	0	cftmdl2(128, &a[384], &w[nw - 128]);
2389	0	cftf162(&a[480], &w[nw - 32]);
2390	0	}
2391	0	cftf161(&a[384], &w[nw - 8]);
2392	0	cftf162(&a[416], &w[nw - 32]);
2393	0	cftf161(&a[448], &w[nw - 8]);
2394	0	} else {
2395	0	cftmdl1(64, a, &w[nw - 32]);
2396	0	cftf081(a, &w[nw - 8]);
2397	0	cftf082(&a[16], &w[nw - 8]);
2398	0	cftf081(&a[32], &w[nw - 8]);
2399	0	cftf081(&a[48], &w[nw - 8]);
2400	0	cftmdl2(64, &a[64], &w[nw - 64]);
2401	0	cftf081(&a[64], &w[nw - 8]);
2402	0	cftf082(&a[80], &w[nw - 8]);
2403	0	cftf081(&a[96], &w[nw - 8]);
2404	0	cftf082(&a[112], &w[nw - 8]);
2405	0	cftmdl1(64, &a[128], &w[nw - 32]);
2406	0	cftf081(&a[128], &w[nw - 8]);
2407	0	cftf082(&a[144], &w[nw - 8]);
2408	0	cftf081(&a[160], &w[nw - 8]);
2409	0	cftf081(&a[176], &w[nw - 8]);
2410	0	if (isplt != 0) {
2411	0	cftmdl1(64, &a[192], &w[nw - 32]);
2412	0	cftf081(&a[240], &w[nw - 8]);
2413	0	} else {
2414	0	cftmdl2(64, &a[192], &w[nw - 64]);
2415	0	cftf082(&a[240], &w[nw - 8]);
2416	0	}
2417	0	cftf081(&a[192], &w[nw - 8]);
2418	0	cftf082(&a[208], &w[nw - 8]);
2419	0	cftf081(&a[224], &w[nw - 8]);
2420	0	}
2421	0	}
2422
2423
2424		void cftmdl1(int n, double a, double w)
2425	0	{
2426	0	int j, j0, j1, j2, j3, k, m, mh;
2427	0	double wn4r, wk1r, wk1i, wk3r, wk3i;
2428	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
2429
2430	0	mh = n >> 3;
2431	0	m = 2 * mh;
2432	0	j1 = m;
2433	0	j2 = j1 + m;
2434	0	j3 = j2 + m;
2435	0	x0r = a[0] + a[j2];
2436	0	x0i = a[1] + a[j2 + 1];
2437	0	x1r = a[0] - a[j2];
2438	0	x1i = a[1] - a[j2 + 1];
2439	0	x2r = a[j1] + a[j3];
2440	0	x2i = a[j1 + 1] + a[j3 + 1];
2441	0	x3r = a[j1] - a[j3];
2442	0	x3i = a[j1 + 1] - a[j3 + 1];
2443	0	a[0] = x0r + x2r;
2444	0	a[1] = x0i + x2i;
2445	0	a[j1] = x0r - x2r;
2446	0	a[j1 + 1] = x0i - x2i;
2447	0	a[j2] = x1r - x3i;
2448	0	a[j2 + 1] = x1i + x3r;
2449	0	a[j3] = x1r + x3i;
2450	0	a[j3 + 1] = x1i - x3r;
2451	0	wn4r = w[1];
2452	0	k = 0;
2453	0	for (j = 2; j < mh; j += 2) {
2454	0	k += 4;
2455	0	wk1r = w[k];
2456	0	wk1i = w[k + 1];
2457	0	wk3r = w[k + 2];
2458	0	wk3i = w[k + 3];
2459	0	j1 = j + m;
2460	0	j2 = j1 + m;
2461	0	j3 = j2 + m;
2462	0	x0r = a[j] + a[j2];
2463	0	x0i = a[j + 1] + a[j2 + 1];
2464	0	x1r = a[j] - a[j2];
2465	0	x1i = a[j + 1] - a[j2 + 1];
2466	0	x2r = a[j1] + a[j3];
2467	0	x2i = a[j1 + 1] + a[j3 + 1];
2468	0	x3r = a[j1] - a[j3];
2469	0	x3i = a[j1 + 1] - a[j3 + 1];
2470	0	a[j] = x0r + x2r;
2471	0	a[j + 1] = x0i + x2i;
2472	0	a[j1] = x0r - x2r;
2473	0	a[j1 + 1] = x0i - x2i;
2474	0	x0r = x1r - x3i;
2475	0	x0i = x1i + x3r;
2476	0	a[j2] = wk1r * x0r - wk1i * x0i;
2477	0	a[j2 + 1] = wk1r * x0i + wk1i * x0r;
2478	0	x0r = x1r + x3i;
2479	0	x0i = x1i - x3r;
2480	0	a[j3] = wk3r * x0r + wk3i * x0i;
2481	0	a[j3 + 1] = wk3r * x0i - wk3i * x0r;
2482	0	j0 = m - j;
2483	0	j1 = j0 + m;
2484	0	j2 = j1 + m;
2485	0	j3 = j2 + m;
2486	0	x0r = a[j0] + a[j2];
2487	0	x0i = a[j0 + 1] + a[j2 + 1];
2488	0	x1r = a[j0] - a[j2];
2489	0	x1i = a[j0 + 1] - a[j2 + 1];
2490	0	x2r = a[j1] + a[j3];
2491	0	x2i = a[j1 + 1] + a[j3 + 1];
2492	0	x3r = a[j1] - a[j3];
2493	0	x3i = a[j1 + 1] - a[j3 + 1];
2494	0	a[j0] = x0r + x2r;
2495	0	a[j0 + 1] = x0i + x2i;
2496	0	a[j1] = x0r - x2r;
2497	0	a[j1 + 1] = x0i - x2i;
2498	0	x0r = x1r - x3i;
2499	0	x0i = x1i + x3r;
2500	0	a[j2] = wk1i * x0r - wk1r * x0i;
2501	0	a[j2 + 1] = wk1i * x0i + wk1r * x0r;
2502	0	x0r = x1r + x3i;
2503	0	x0i = x1i - x3r;
2504	0	a[j3] = wk3i * x0r + wk3r * x0i;
2505	0	a[j3 + 1] = wk3i * x0i - wk3r * x0r;
2506	0	}
2507	0	j0 = mh;
2508	0	j1 = j0 + m;
2509	0	j2 = j1 + m;
2510	0	j3 = j2 + m;
2511	0	x0r = a[j0] + a[j2];
2512	0	x0i = a[j0 + 1] + a[j2 + 1];
2513	0	x1r = a[j0] - a[j2];
2514	0	x1i = a[j0 + 1] - a[j2 + 1];
2515	0	x2r = a[j1] + a[j3];
2516	0	x2i = a[j1 + 1] + a[j3 + 1];
2517	0	x3r = a[j1] - a[j3];
2518	0	x3i = a[j1 + 1] - a[j3 + 1];
2519	0	a[j0] = x0r + x2r;
2520	0	a[j0 + 1] = x0i + x2i;
2521	0	a[j1] = x0r - x2r;
2522	0	a[j1 + 1] = x0i - x2i;
2523	0	x0r = x1r - x3i;
2524	0	x0i = x1i + x3r;
2525	0	a[j2] = wn4r * (x0r - x0i);
2526	0	a[j2 + 1] = wn4r * (x0i + x0r);
2527	0	x0r = x1r + x3i;
2528	0	x0i = x1i - x3r;
2529	0	a[j3] = -wn4r * (x0r + x0i);
2530	0	a[j3 + 1] = -wn4r * (x0i - x0r);
2531	0	}
2532
2533
2534		void cftmdl2(int n, double a, double w)
2535	0	{
2536	0	int j, j0, j1, j2, j3, k, kr, m, mh;
2537	0	double wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
2538	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, y0r, y0i, y2r, y2i;
2539
2540	0	mh = n >> 3;
2541	0	m = 2 * mh;
2542	0	wn4r = w[1];
2543	0	j1 = m;
2544	0	j2 = j1 + m;
2545	0	j3 = j2 + m;
2546	0	x0r = a[0] - a[j2 + 1];
2547	0	x0i = a[1] + a[j2];
2548	0	x1r = a[0] + a[j2 + 1];
2549	0	x1i = a[1] - a[j2];
2550	0	x2r = a[j1] - a[j3 + 1];
2551	0	x2i = a[j1 + 1] + a[j3];
2552	0	x3r = a[j1] + a[j3 + 1];
2553	0	x3i = a[j1 + 1] - a[j3];
2554	0	y0r = wn4r * (x2r - x2i);
2555	0	y0i = wn4r * (x2i + x2r);
2556	0	a[0] = x0r + y0r;
2557	0	a[1] = x0i + y0i;
2558	0	a[j1] = x0r - y0r;
2559	0	a[j1 + 1] = x0i - y0i;
2560	0	y0r = wn4r * (x3r - x3i);
2561	0	y0i = wn4r * (x3i + x3r);
2562	0	a[j2] = x1r - y0i;
2563	0	a[j2 + 1] = x1i + y0r;
2564	0	a[j3] = x1r + y0i;
2565	0	a[j3 + 1] = x1i - y0r;
2566	0	k = 0;
2567	0	kr = 2 * m;
2568	0	for (j = 2; j < mh; j += 2) {
2569	0	k += 4;
2570	0	wk1r = w[k];
2571	0	wk1i = w[k + 1];
2572	0	wk3r = w[k + 2];
2573	0	wk3i = w[k + 3];
2574	0	kr -= 4;
2575	0	wd1i = w[kr];
2576	0	wd1r = w[kr + 1];
2577	0	wd3i = w[kr + 2];
2578	0	wd3r = w[kr + 3];
2579	0	j1 = j + m;
2580	0	j2 = j1 + m;
2581	0	j3 = j2 + m;
2582	0	x0r = a[j] - a[j2 + 1];
2583	0	x0i = a[j + 1] + a[j2];
2584	0	x1r = a[j] + a[j2 + 1];
2585	0	x1i = a[j + 1] - a[j2];
2586	0	x2r = a[j1] - a[j3 + 1];
2587	0	x2i = a[j1 + 1] + a[j3];
2588	0	x3r = a[j1] + a[j3 + 1];
2589	0	x3i = a[j1 + 1] - a[j3];
2590	0	y0r = wk1r * x0r - wk1i * x0i;
2591	0	y0i = wk1r * x0i + wk1i * x0r;
2592	0	y2r = wd1r * x2r - wd1i * x2i;
2593	0	y2i = wd1r * x2i + wd1i * x2r;
2594	0	a[j] = y0r + y2r;
2595	0	a[j + 1] = y0i + y2i;
2596	0	a[j1] = y0r - y2r;
2597	0	a[j1 + 1] = y0i - y2i;
2598	0	y0r = wk3r * x1r + wk3i * x1i;
2599	0	y0i = wk3r * x1i - wk3i * x1r;
2600	0	y2r = wd3r * x3r + wd3i * x3i;
2601	0	y2i = wd3r * x3i - wd3i * x3r;
2602	0	a[j2] = y0r + y2r;
2603	0	a[j2 + 1] = y0i + y2i;
2604	0	a[j3] = y0r - y2r;
2605	0	a[j3 + 1] = y0i - y2i;
2606	0	j0 = m - j;
2607	0	j1 = j0 + m;
2608	0	j2 = j1 + m;
2609	0	j3 = j2 + m;
2610	0	x0r = a[j0] - a[j2 + 1];
2611	0	x0i = a[j0 + 1] + a[j2];
2612	0	x1r = a[j0] + a[j2 + 1];
2613	0	x1i = a[j0 + 1] - a[j2];
2614	0	x2r = a[j1] - a[j3 + 1];
2615	0	x2i = a[j1 + 1] + a[j3];
2616	0	x3r = a[j1] + a[j3 + 1];
2617	0	x3i = a[j1 + 1] - a[j3];
2618	0	y0r = wd1i * x0r - wd1r * x0i;
2619	0	y0i = wd1i * x0i + wd1r * x0r;
2620	0	y2r = wk1i * x2r - wk1r * x2i;
2621	0	y2i = wk1i * x2i + wk1r * x2r;
2622	0	a[j0] = y0r + y2r;
2623	0	a[j0 + 1] = y0i + y2i;
2624	0	a[j1] = y0r - y2r;
2625	0	a[j1 + 1] = y0i - y2i;
2626	0	y0r = wd3i * x1r + wd3r * x1i;
2627	0	y0i = wd3i * x1i - wd3r * x1r;
2628	0	y2r = wk3i * x3r + wk3r * x3i;
2629	0	y2i = wk3i * x3i - wk3r * x3r;
2630	0	a[j2] = y0r + y2r;
2631	0	a[j2 + 1] = y0i + y2i;
2632	0	a[j3] = y0r - y2r;
2633	0	a[j3 + 1] = y0i - y2i;
2634	0	}
2635	0	wk1r = w[m];
2636	0	wk1i = w[m + 1];
2637	0	j0 = mh;
2638	0	j1 = j0 + m;
2639	0	j2 = j1 + m;
2640	0	j3 = j2 + m;
2641	0	x0r = a[j0] - a[j2 + 1];
2642	0	x0i = a[j0 + 1] + a[j2];
2643	0	x1r = a[j0] + a[j2 + 1];
2644	0	x1i = a[j0 + 1] - a[j2];
2645	0	x2r = a[j1] - a[j3 + 1];
2646	0	x2i = a[j1 + 1] + a[j3];
2647	0	x3r = a[j1] + a[j3 + 1];
2648	0	x3i = a[j1 + 1] - a[j3];
2649	0	y0r = wk1r * x0r - wk1i * x0i;
2650	0	y0i = wk1r * x0i + wk1i * x0r;
2651	0	y2r = wk1i * x2r - wk1r * x2i;
2652	0	y2i = wk1i * x2i + wk1r * x2r;
2653	0	a[j0] = y0r + y2r;
2654	0	a[j0 + 1] = y0i + y2i;
2655	0	a[j1] = y0r - y2r;
2656	0	a[j1 + 1] = y0i - y2i;
2657	0	y0r = wk1i * x1r - wk1r * x1i;
2658	0	y0i = wk1i * x1i + wk1r * x1r;
2659	0	y2r = wk1r * x3r - wk1i * x3i;
2660	0	y2i = wk1r * x3i + wk1i * x3r;
2661	0	a[j2] = y0r - y2r;
2662	0	a[j2 + 1] = y0i - y2i;
2663	0	a[j3] = y0r + y2r;
2664	0	a[j3 + 1] = y0i + y2i;
2665	0	}
2666
2667
2668		void cftfx41(int n, double a, int nw, double w)
2669	0	{
2670	0	void cftf161(double a, double w);
2671	0	void cftf162(double a, double w);
2672	0	void cftf081(double a, double w);
2673	0	void cftf082(double a, double w);
2674
2675	0	if (n == 128) {
2676	0	cftf161(a, &w[nw - 8]);
2677	0	cftf162(&a[32], &w[nw - 32]);
2678	0	cftf161(&a[64], &w[nw - 8]);
2679	0	cftf161(&a[96], &w[nw - 8]);
2680	0	} else {
2681	0	cftf081(a, &w[nw - 8]);
2682	0	cftf082(&a[16], &w[nw - 8]);
2683	0	cftf081(&a[32], &w[nw - 8]);
2684	0	cftf081(&a[48], &w[nw - 8]);
2685	0	}
2686	0	}
2687
2688
2689		void cftf161(double a, double w)
2690	0	{
2691	0	double wn4r, wk1r, wk1i,
2692	0	x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
2693	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2694	0	y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2695	0	y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2696	0	y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2697
2698	0	wn4r = w[1];
2699	0	wk1r = w[2];
2700	0	wk1i = w[3];
2701	0	x0r = a[0] + a[16];
2702	0	x0i = a[1] + a[17];
2703	0	x1r = a[0] - a[16];
2704	0	x1i = a[1] - a[17];
2705	0	x2r = a[8] + a[24];
2706	0	x2i = a[9] + a[25];
2707	0	x3r = a[8] - a[24];
2708	0	x3i = a[9] - a[25];
2709	0	y0r = x0r + x2r;
2710	0	y0i = x0i + x2i;
2711	0	y4r = x0r - x2r;
2712	0	y4i = x0i - x2i;
2713	0	y8r = x1r - x3i;
2714	0	y8i = x1i + x3r;
2715	0	y12r = x1r + x3i;
2716	0	y12i = x1i - x3r;
2717	0	x0r = a[2] + a[18];
2718	0	x0i = a[3] + a[19];
2719	0	x1r = a[2] - a[18];
2720	0	x1i = a[3] - a[19];
2721	0	x2r = a[10] + a[26];
2722	0	x2i = a[11] + a[27];
2723	0	x3r = a[10] - a[26];
2724	0	x3i = a[11] - a[27];
2725	0	y1r = x0r + x2r;
2726	0	y1i = x0i + x2i;
2727	0	y5r = x0r - x2r;
2728	0	y5i = x0i - x2i;
2729	0	x0r = x1r - x3i;
2730	0	x0i = x1i + x3r;
2731	0	y9r = wk1r * x0r - wk1i * x0i;
2732	0	y9i = wk1r * x0i + wk1i * x0r;
2733	0	x0r = x1r + x3i;
2734	0	x0i = x1i - x3r;
2735	0	y13r = wk1i * x0r - wk1r * x0i;
2736	0	y13i = wk1i * x0i + wk1r * x0r;
2737	0	x0r = a[4] + a[20];
2738	0	x0i = a[5] + a[21];
2739	0	x1r = a[4] - a[20];
2740	0	x1i = a[5] - a[21];
2741	0	x2r = a[12] + a[28];
2742	0	x2i = a[13] + a[29];
2743	0	x3r = a[12] - a[28];
2744	0	x3i = a[13] - a[29];
2745	0	y2r = x0r + x2r;
2746	0	y2i = x0i + x2i;
2747	0	y6r = x0r - x2r;
2748	0	y6i = x0i - x2i;
2749	0	x0r = x1r - x3i;
2750	0	x0i = x1i + x3r;
2751	0	y10r = wn4r * (x0r - x0i);
2752	0	y10i = wn4r * (x0i + x0r);
2753	0	x0r = x1r + x3i;
2754	0	x0i = x1i - x3r;
2755	0	y14r = wn4r * (x0r + x0i);
2756	0	y14i = wn4r * (x0i - x0r);
2757	0	x0r = a[6] + a[22];
2758	0	x0i = a[7] + a[23];
2759	0	x1r = a[6] - a[22];
2760	0	x1i = a[7] - a[23];
2761	0	x2r = a[14] + a[30];
2762	0	x2i = a[15] + a[31];
2763	0	x3r = a[14] - a[30];
2764	0	x3i = a[15] - a[31];
2765	0	y3r = x0r + x2r;
2766	0	y3i = x0i + x2i;
2767	0	y7r = x0r - x2r;
2768	0	y7i = x0i - x2i;
2769	0	x0r = x1r - x3i;
2770	0	x0i = x1i + x3r;
2771	0	y11r = wk1i * x0r - wk1r * x0i;
2772	0	y11i = wk1i * x0i + wk1r * x0r;
2773	0	x0r = x1r + x3i;
2774	0	x0i = x1i - x3r;
2775	0	y15r = wk1r * x0r - wk1i * x0i;
2776	0	y15i = wk1r * x0i + wk1i * x0r;
2777	0	x0r = y12r - y14r;
2778	0	x0i = y12i - y14i;
2779	0	x1r = y12r + y14r;
2780	0	x1i = y12i + y14i;
2781	0	x2r = y13r - y15r;
2782	0	x2i = y13i - y15i;
2783	0	x3r = y13r + y15r;
2784	0	x3i = y13i + y15i;
2785	0	a[24] = x0r + x2r;
2786	0	a[25] = x0i + x2i;
2787	0	a[26] = x0r - x2r;
2788	0	a[27] = x0i - x2i;
2789	0	a[28] = x1r - x3i;
2790	0	a[29] = x1i + x3r;
2791	0	a[30] = x1r + x3i;
2792	0	a[31] = x1i - x3r;
2793	0	x0r = y8r + y10r;
2794	0	x0i = y8i + y10i;
2795	0	x1r = y8r - y10r;
2796	0	x1i = y8i - y10i;
2797	0	x2r = y9r + y11r;
2798	0	x2i = y9i + y11i;
2799	0	x3r = y9r - y11r;
2800	0	x3i = y9i - y11i;
2801	0	a[16] = x0r + x2r;
2802	0	a[17] = x0i + x2i;
2803	0	a[18] = x0r - x2r;
2804	0	a[19] = x0i - x2i;
2805	0	a[20] = x1r - x3i;
2806	0	a[21] = x1i + x3r;
2807	0	a[22] = x1r + x3i;
2808	0	a[23] = x1i - x3r;
2809	0	x0r = y5r - y7i;
2810	0	x0i = y5i + y7r;
2811	0	x2r = wn4r * (x0r - x0i);
2812	0	x2i = wn4r * (x0i + x0r);
2813	0	x0r = y5r + y7i;
2814	0	x0i = y5i - y7r;
2815	0	x3r = wn4r * (x0r - x0i);
2816	0	x3i = wn4r * (x0i + x0r);
2817	0	x0r = y4r - y6i;
2818	0	x0i = y4i + y6r;
2819	0	x1r = y4r + y6i;
2820	0	x1i = y4i - y6r;
2821	0	a[8] = x0r + x2r;
2822	0	a[9] = x0i + x2i;
2823	0	a[10] = x0r - x2r;
2824	0	a[11] = x0i - x2i;
2825	0	a[12] = x1r - x3i;
2826	0	a[13] = x1i + x3r;
2827	0	a[14] = x1r + x3i;
2828	0	a[15] = x1i - x3r;
2829	0	x0r = y0r + y2r;
2830	0	x0i = y0i + y2i;
2831	0	x1r = y0r - y2r;
2832	0	x1i = y0i - y2i;
2833	0	x2r = y1r + y3r;
2834	0	x2i = y1i + y3i;
2835	0	x3r = y1r - y3r;
2836	0	x3i = y1i - y3i;
2837	0	a[0] = x0r + x2r;
2838	0	a[1] = x0i + x2i;
2839	0	a[2] = x0r - x2r;
2840	0	a[3] = x0i - x2i;
2841	0	a[4] = x1r - x3i;
2842	0	a[5] = x1i + x3r;
2843	0	a[6] = x1r + x3i;
2844	0	a[7] = x1i - x3r;
2845	0	}
2846
2847
2848		void cftf162(double a, double w)
2849	0	{
2850	0	double wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i,
2851	0	x0r, x0i, x1r, x1i, x2r, x2i,
2852	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
2853	0	y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i,
2854	0	y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i,
2855	0	y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i;
2856
2857	0	wn4r = w[1];
2858	0	wk1r = w[4];
2859	0	wk1i = w[5];
2860	0	wk3r = w[6];
2861	0	wk3i = -w[7];
2862	0	wk2r = w[8];
2863	0	wk2i = w[9];
2864	0	x1r = a[0] - a[17];
2865	0	x1i = a[1] + a[16];
2866	0	x0r = a[8] - a[25];
2867	0	x0i = a[9] + a[24];
2868	0	x2r = wn4r * (x0r - x0i);
2869	0	x2i = wn4r * (x0i + x0r);
2870	0	y0r = x1r + x2r;
2871	0	y0i = x1i + x2i;
2872	0	y4r = x1r - x2r;
2873	0	y4i = x1i - x2i;
2874	0	x1r = a[0] + a[17];
2875	0	x1i = a[1] - a[16];
2876	0	x0r = a[8] + a[25];
2877	0	x0i = a[9] - a[24];
2878	0	x2r = wn4r * (x0r - x0i);
2879	0	x2i = wn4r * (x0i + x0r);
2880	0	y8r = x1r - x2i;
2881	0	y8i = x1i + x2r;
2882	0	y12r = x1r + x2i;
2883	0	y12i = x1i - x2r;
2884	0	x0r = a[2] - a[19];
2885	0	x0i = a[3] + a[18];
2886	0	x1r = wk1r * x0r - wk1i * x0i;
2887	0	x1i = wk1r * x0i + wk1i * x0r;
2888	0	x0r = a[10] - a[27];
2889	0	x0i = a[11] + a[26];
2890	0	x2r = wk3i * x0r - wk3r * x0i;
2891	0	x2i = wk3i * x0i + wk3r * x0r;
2892	0	y1r = x1r + x2r;
2893	0	y1i = x1i + x2i;
2894	0	y5r = x1r - x2r;
2895	0	y5i = x1i - x2i;
2896	0	x0r = a[2] + a[19];
2897	0	x0i = a[3] - a[18];
2898	0	x1r = wk3r * x0r - wk3i * x0i;
2899	0	x1i = wk3r * x0i + wk3i * x0r;
2900	0	x0r = a[10] + a[27];
2901	0	x0i = a[11] - a[26];
2902	0	x2r = wk1r * x0r + wk1i * x0i;
2903	0	x2i = wk1r * x0i - wk1i * x0r;
2904	0	y9r = x1r - x2r;
2905	0	y9i = x1i - x2i;
2906	0	y13r = x1r + x2r;
2907	0	y13i = x1i + x2i;
2908	0	x0r = a[4] - a[21];
2909	0	x0i = a[5] + a[20];
2910	0	x1r = wk2r * x0r - wk2i * x0i;
2911	0	x1i = wk2r * x0i + wk2i * x0r;
2912	0	x0r = a[12] - a[29];
2913	0	x0i = a[13] + a[28];
2914	0	x2r = wk2i * x0r - wk2r * x0i;
2915	0	x2i = wk2i * x0i + wk2r * x0r;
2916	0	y2r = x1r + x2r;
2917	0	y2i = x1i + x2i;
2918	0	y6r = x1r - x2r;
2919	0	y6i = x1i - x2i;
2920	0	x0r = a[4] + a[21];
2921	0	x0i = a[5] - a[20];
2922	0	x1r = wk2i * x0r - wk2r * x0i;
2923	0	x1i = wk2i * x0i + wk2r * x0r;
2924	0	x0r = a[12] + a[29];
2925	0	x0i = a[13] - a[28];
2926	0	x2r = wk2r * x0r - wk2i * x0i;
2927	0	x2i = wk2r * x0i + wk2i * x0r;
2928	0	y10r = x1r - x2r;
2929	0	y10i = x1i - x2i;
2930	0	y14r = x1r + x2r;
2931	0	y14i = x1i + x2i;
2932	0	x0r = a[6] - a[23];
2933	0	x0i = a[7] + a[22];
2934	0	x1r = wk3r * x0r - wk3i * x0i;
2935	0	x1i = wk3r * x0i + wk3i * x0r;
2936	0	x0r = a[14] - a[31];
2937	0	x0i = a[15] + a[30];
2938	0	x2r = wk1i * x0r - wk1r * x0i;
2939	0	x2i = wk1i * x0i + wk1r * x0r;
2940	0	y3r = x1r + x2r;
2941	0	y3i = x1i + x2i;
2942	0	y7r = x1r - x2r;
2943	0	y7i = x1i - x2i;
2944	0	x0r = a[6] + a[23];
2945	0	x0i = a[7] - a[22];
2946	0	x1r = wk1i * x0r + wk1r * x0i;
2947	0	x1i = wk1i * x0i - wk1r * x0r;
2948	0	x0r = a[14] + a[31];
2949	0	x0i = a[15] - a[30];
2950	0	x2r = wk3i * x0r - wk3r * x0i;
2951	0	x2i = wk3i * x0i + wk3r * x0r;
2952	0	y11r = x1r + x2r;
2953	0	y11i = x1i + x2i;
2954	0	y15r = x1r - x2r;
2955	0	y15i = x1i - x2i;
2956	0	x1r = y0r + y2r;
2957	0	x1i = y0i + y2i;
2958	0	x2r = y1r + y3r;
2959	0	x2i = y1i + y3i;
2960	0	a[0] = x1r + x2r;
2961	0	a[1] = x1i + x2i;
2962	0	a[2] = x1r - x2r;
2963	0	a[3] = x1i - x2i;
2964	0	x1r = y0r - y2r;
2965	0	x1i = y0i - y2i;
2966	0	x2r = y1r - y3r;
2967	0	x2i = y1i - y3i;
2968	0	a[4] = x1r - x2i;
2969	0	a[5] = x1i + x2r;
2970	0	a[6] = x1r + x2i;
2971	0	a[7] = x1i - x2r;
2972	0	x1r = y4r - y6i;
2973	0	x1i = y4i + y6r;
2974	0	x0r = y5r - y7i;
2975	0	x0i = y5i + y7r;
2976	0	x2r = wn4r * (x0r - x0i);
2977	0	x2i = wn4r * (x0i + x0r);
2978	0	a[8] = x1r + x2r;
2979	0	a[9] = x1i + x2i;
2980	0	a[10] = x1r - x2r;
2981	0	a[11] = x1i - x2i;
2982	0	x1r = y4r + y6i;
2983	0	x1i = y4i - y6r;
2984	0	x0r = y5r + y7i;
2985	0	x0i = y5i - y7r;
2986	0	x2r = wn4r * (x0r - x0i);
2987	0	x2i = wn4r * (x0i + x0r);
2988	0	a[12] = x1r - x2i;
2989	0	a[13] = x1i + x2r;
2990	0	a[14] = x1r + x2i;
2991	0	a[15] = x1i - x2r;
2992	0	x1r = y8r + y10r;
2993	0	x1i = y8i + y10i;
2994	0	x2r = y9r - y11r;
2995	0	x2i = y9i - y11i;
2996	0	a[16] = x1r + x2r;
2997	0	a[17] = x1i + x2i;
2998	0	a[18] = x1r - x2r;
2999	0	a[19] = x1i - x2i;
3000	0	x1r = y8r - y10r;
3001	0	x1i = y8i - y10i;
3002	0	x2r = y9r + y11r;
3003	0	x2i = y9i + y11i;
3004	0	a[20] = x1r - x2i;
3005	0	a[21] = x1i + x2r;
3006	0	a[22] = x1r + x2i;
3007	0	a[23] = x1i - x2r;
3008	0	x1r = y12r - y14i;
3009	0	x1i = y12i + y14r;
3010	0	x0r = y13r + y15i;
3011	0	x0i = y13i - y15r;
3012	0	x2r = wn4r * (x0r - x0i);
3013	0	x2i = wn4r * (x0i + x0r);
3014	0	a[24] = x1r + x2r;
3015	0	a[25] = x1i + x2i;
3016	0	a[26] = x1r - x2r;
3017	0	a[27] = x1i - x2i;
3018	0	x1r = y12r + y14i;
3019	0	x1i = y12i - y14r;
3020	0	x0r = y13r - y15i;
3021	0	x0i = y13i + y15r;
3022	0	x2r = wn4r * (x0r - x0i);
3023	0	x2i = wn4r * (x0i + x0r);
3024	0	a[28] = x1r - x2i;
3025	0	a[29] = x1i + x2r;
3026	0	a[30] = x1r + x2i;
3027	0	a[31] = x1i - x2r;
3028	0	}
3029
3030
3031		void cftf081(double a, double w)
3032	0	{
3033	0	double wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i,
3034	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3035	0	y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3036
3037	0	wn4r = w[1];
3038	0	x0r = a[0] + a[8];
3039	0	x0i = a[1] + a[9];
3040	0	x1r = a[0] - a[8];
3041	0	x1i = a[1] - a[9];
3042	0	x2r = a[4] + a[12];
3043	0	x2i = a[5] + a[13];
3044	0	x3r = a[4] - a[12];
3045	0	x3i = a[5] - a[13];
3046	0	y0r = x0r + x2r;
3047	0	y0i = x0i + x2i;
3048	0	y2r = x0r - x2r;
3049	0	y2i = x0i - x2i;
3050	0	y1r = x1r - x3i;
3051	0	y1i = x1i + x3r;
3052	0	y3r = x1r + x3i;
3053	0	y3i = x1i - x3r;
3054	0	x0r = a[2] + a[10];
3055	0	x0i = a[3] + a[11];
3056	0	x1r = a[2] - a[10];
3057	0	x1i = a[3] - a[11];
3058	0	x2r = a[6] + a[14];
3059	0	x2i = a[7] + a[15];
3060	0	x3r = a[6] - a[14];
3061	0	x3i = a[7] - a[15];
3062	0	y4r = x0r + x2r;
3063	0	y4i = x0i + x2i;
3064	0	y6r = x0r - x2r;
3065	0	y6i = x0i - x2i;
3066	0	x0r = x1r - x3i;
3067	0	x0i = x1i + x3r;
3068	0	x2r = x1r + x3i;
3069	0	x2i = x1i - x3r;
3070	0	y5r = wn4r * (x0r - x0i);
3071	0	y5i = wn4r * (x0r + x0i);
3072	0	y7r = wn4r * (x2r - x2i);
3073	0	y7i = wn4r * (x2r + x2i);
3074	0	a[8] = y1r + y5r;
3075	0	a[9] = y1i + y5i;
3076	0	a[10] = y1r - y5r;
3077	0	a[11] = y1i - y5i;
3078	0	a[12] = y3r - y7i;
3079	0	a[13] = y3i + y7r;
3080	0	a[14] = y3r + y7i;
3081	0	a[15] = y3i - y7r;
3082	0	a[0] = y0r + y4r;
3083	0	a[1] = y0i + y4i;
3084	0	a[2] = y0r - y4r;
3085	0	a[3] = y0i - y4i;
3086	0	a[4] = y2r - y6i;
3087	0	a[5] = y2i + y6r;
3088	0	a[6] = y2r + y6i;
3089	0	a[7] = y2i - y6r;
3090	0	}
3091
3092
3093		void cftf082(double a, double w)
3094	0	{
3095	0	double wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i,
3096	0	y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i,
3097	0	y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i;
3098
3099	0	wn4r = w[1];
3100	0	wk1r = w[2];
3101	0	wk1i = w[3];
3102	0	y0r = a[0] - a[9];
3103	0	y0i = a[1] + a[8];
3104	0	y1r = a[0] + a[9];
3105	0	y1i = a[1] - a[8];
3106	0	x0r = a[4] - a[13];
3107	0	x0i = a[5] + a[12];
3108	0	y2r = wn4r * (x0r - x0i);
3109	0	y2i = wn4r * (x0i + x0r);
3110	0	x0r = a[4] + a[13];
3111	0	x0i = a[5] - a[12];
3112	0	y3r = wn4r * (x0r - x0i);
3113	0	y3i = wn4r * (x0i + x0r);
3114	0	x0r = a[2] - a[11];
3115	0	x0i = a[3] + a[10];
3116	0	y4r = wk1r * x0r - wk1i * x0i;
3117	0	y4i = wk1r * x0i + wk1i * x0r;
3118	0	x0r = a[2] + a[11];
3119	0	x0i = a[3] - a[10];
3120	0	y5r = wk1i * x0r - wk1r * x0i;
3121	0	y5i = wk1i * x0i + wk1r * x0r;
3122	0	x0r = a[6] - a[15];
3123	0	x0i = a[7] + a[14];
3124	0	y6r = wk1i * x0r - wk1r * x0i;
3125	0	y6i = wk1i * x0i + wk1r * x0r;
3126	0	x0r = a[6] + a[15];
3127	0	x0i = a[7] - a[14];
3128	0	y7r = wk1r * x0r - wk1i * x0i;
3129	0	y7i = wk1r * x0i + wk1i * x0r;
3130	0	x0r = y0r + y2r;
3131	0	x0i = y0i + y2i;
3132	0	x1r = y4r + y6r;
3133	0	x1i = y4i + y6i;
3134	0	a[0] = x0r + x1r;
3135	0	a[1] = x0i + x1i;
3136	0	a[2] = x0r - x1r;
3137	0	a[3] = x0i - x1i;
3138	0	x0r = y0r - y2r;
3139	0	x0i = y0i - y2i;
3140	0	x1r = y4r - y6r;
3141	0	x1i = y4i - y6i;
3142	0	a[4] = x0r - x1i;
3143	0	a[5] = x0i + x1r;
3144	0	a[6] = x0r + x1i;
3145	0	a[7] = x0i - x1r;
3146	0	x0r = y1r - y3i;
3147	0	x0i = y1i + y3r;
3148	0	x1r = y5r - y7r;
3149	0	x1i = y5i - y7i;
3150	0	a[8] = x0r + x1r;
3151	0	a[9] = x0i + x1i;
3152	0	a[10] = x0r - x1r;
3153	0	a[11] = x0i - x1i;
3154	0	x0r = y1r + y3i;
3155	0	x0i = y1i - y3r;
3156	0	x1r = y5r + y7r;
3157	0	x1i = y5i + y7i;
3158	0	a[12] = x0r - x1i;
3159	0	a[13] = x0i + x1r;
3160	0	a[14] = x0r + x1i;
3161	0	a[15] = x0i - x1r;
3162	0	}
3163
3164
3165		void cftf040(double *a)
3166	0	{
3167	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3168
3169	0	x0r = a[0] + a[4];
3170	0	x0i = a[1] + a[5];
3171	0	x1r = a[0] - a[4];
3172	0	x1i = a[1] - a[5];
3173	0	x2r = a[2] + a[6];
3174	0	x2i = a[3] + a[7];
3175	0	x3r = a[2] - a[6];
3176	0	x3i = a[3] - a[7];
3177	0	a[0] = x0r + x2r;
3178	0	a[1] = x0i + x2i;
3179	0	a[2] = x1r - x3i;
3180	0	a[3] = x1i + x3r;
3181	0	a[4] = x0r - x2r;
3182	0	a[5] = x0i - x2i;
3183	0	a[6] = x1r + x3i;
3184	0	a[7] = x1i - x3r;
3185	0	}
3186
3187
3188		void cftb040(double *a)
3189	0	{
3190	0	double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
3191
3192	0	x0r = a[0] + a[4];
3193	0	x0i = a[1] + a[5];
3194	0	x1r = a[0] - a[4];
3195	0	x1i = a[1] - a[5];
3196	0	x2r = a[2] + a[6];
3197	0	x2i = a[3] + a[7];
3198	0	x3r = a[2] - a[6];
3199	0	x3i = a[3] - a[7];
3200	0	a[0] = x0r + x2r;
3201	0	a[1] = x0i + x2i;
3202	0	a[2] = x1r + x3i;
3203	0	a[3] = x1i - x3r;
3204	0	a[4] = x0r - x2r;
3205	0	a[5] = x0i - x2i;
3206	0	a[6] = x1r - x3i;
3207	0	a[7] = x1i + x3r;
3208	0	}
3209
3210
3211		void cftx020(double *a)
3212	0	{
3213	0	double x0r, x0i;
3214
3215	0	x0r = a[0] - a[2];
3216	0	x0i = a[1] - a[3];
3217	0	a[0] += a[2];
3218	0	a[1] += a[3];
3219	0	a[2] = x0r;
3220	0	a[3] = x0i;
3221	0	}
3222
3223
3224		void rftfsub(int n, double a, int nc, double c)
3225	0	{
3226	0	int j, k, kk, ks, m;
3227	0	double wkr, wki, xr, xi, yr, yi;
3228
3229	0	m = n >> 1;
3230	0	ks = 2 * nc / m;
3231	0	kk = 0;
3232	0	for (j = 2; j < m; j += 2) {
3233	0	k = n - j;
3234	0	kk += ks;
3235	0	wkr = 0.5 - c[nc - kk];
3236	0	wki = c[kk];
3237	0	xr = a[j] - a[k];
3238	0	xi = a[j + 1] + a[k + 1];
3239	0	yr = wkr * xr - wki * xi;
3240	0	yi = wkr * xi + wki * xr;
3241	0	a[j] -= yr;
3242	0	a[j + 1] -= yi;
3243	0	a[k] += yr;
3244	0	a[k + 1] -= yi;
3245	0	}
3246	0	}
3247
3248
3249		void rftbsub(int n, double a, int nc, double c)
3250	0	{
3251	0	int j, k, kk, ks, m;
3252	0	double wkr, wki, xr, xi, yr, yi;
3253
3254	0	m = n >> 1;
3255	0	ks = 2 * nc / m;
3256	0	kk = 0;
3257	0	for (j = 2; j < m; j += 2) {
3258	0	k = n - j;
3259	0	kk += ks;
3260	0	wkr = 0.5 - c[nc - kk];
3261	0	wki = c[kk];
3262	0	xr = a[j] - a[k];
3263	0	xi = a[j + 1] + a[k + 1];
3264	0	yr = wkr * xr + wki * xi;
3265	0	yi = wkr * xi - wki * xr;
3266	0	a[j] -= yr;
3267	0	a[j + 1] -= yi;
3268	0	a[k] += yr;
3269	0	a[k + 1] -= yi;
3270	0	}
3271	0	}
3272
3273
3274		void dctsub(int n, double a, int nc, double c)
3275	0	{
3276	0	int j, k, kk, ks, m;
3277	0	double wkr, wki, xr;
3278
3279	0	m = n >> 1;
3280	0	ks = nc / n;
3281	0	kk = 0;
3282	0	for (j = 1; j < m; j++) {
3283	0	k = n - j;
3284	0	kk += ks;
3285	0	wkr = c[kk] - c[nc - kk];
3286	0	wki = c[kk] + c[nc - kk];
3287	0	xr = wki * a[j] - wkr * a[k];
3288	0	a[j] = wkr * a[j] + wki * a[k];
3289	0	a[k] = xr;
3290	0	}
3291	0	a[m] *= c[0];
3292	0	}
3293
3294
3295		void dstsub(int n, double a, int nc, double c)
3296	0	{
3297	0	int j, k, kk, ks, m;
3298	0	double wkr, wki, xr;
3299
3300	0	m = n >> 1;
3301	0	ks = nc / n;
3302	0	kk = 0;
3303	0	for (j = 1; j < m; j++) {
3304	0	k = n - j;
3305	0	kk += ks;
3306	0	wkr = c[kk] - c[nc - kk];
3307	0	wki = c[kk] + c[nc - kk];
3308	0	xr = wki * a[k] - wkr * a[j];
3309	0	a[k] = wkr * a[k] + wki * a[j];
3310	0	a[j] = xr;
3311	0	}
3312	0	a[m] *= c[0];
3313	0	}
3314

Coverage Report

Created: 2024-05-04 12:45