split-radix FFT
[libav.git] / libavcodec / i386 / fft_3dn2.c
1 /*
2 * FFT/MDCT transform with Extended 3DNow! optimizations
3 * Copyright (c) 2006 Zuxy MENG Jie, Loren Merritt
4 * Based on fft_sse.c copyright (c) 2002 Fabrice Bellard.
5 *
6 * This file is part of FFmpeg.
7 *
8 * FFmpeg is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * FFmpeg is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with FFmpeg; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 */
22
23 #include "libavutil/x86_cpu.h"
24 #include "libavcodec/dsputil.h"
25
26 #ifdef EMULATE_3DNOWEXT
27 #define ff_fft_calc_3dn2 ff_fft_calc_3dn
28 #define ff_fft_dispatch_3dn2 ff_fft_dispatch_3dn
29 #define ff_fft_dispatch_interleave_3dn2 ff_fft_dispatch_interleave_3dn
30 #define ff_imdct_calc_3dn2 ff_imdct_calc_3dn
31 #define ff_imdct_half_3dn2 ff_imdct_half_3dn
32 #endif
33
34 void ff_fft_dispatch_3dn2(FFTComplex *z, int nbits);
35 void ff_fft_dispatch_interleave_3dn2(FFTComplex *z, int nbits);
36
37 void ff_fft_calc_3dn2(FFTContext *s, FFTComplex *z)
38 {
39 int n = 1<<s->nbits;
40 int i;
41 ff_fft_dispatch_interleave_3dn2(z, s->nbits);
42 asm volatile("femms");
43 if(n <= 8)
44 for(i=0; i<n; i+=2)
45 FFSWAP(FFTSample, z[i].im, z[i+1].re);
46 }
47
48 static void imdct_3dn2(MDCTContext *s, const FFTSample *input, FFTSample *tmp)
49 {
50 long n4, n2, n;
51 x86_reg k;
52 const uint16_t *revtab = s->fft.revtab;
53 const FFTSample *tcos = s->tcos;
54 const FFTSample *tsin = s->tsin;
55 const FFTSample *in1, *in2;
56 FFTComplex *z = (FFTComplex *)tmp;
57
58 n = 1 << s->nbits;
59 n2 = n >> 1;
60 n4 = n >> 2;
61
62 /* pre rotation */
63 in1 = input;
64 in2 = input + n2 - 1;
65 for(k = 0; k < n4; k++) {
66 // FIXME a single block is faster, but gcc 2.95 and 3.4.x on 32bit can't compile it
67 asm volatile(
68 "movd %0, %%mm0 \n\t"
69 "movd %2, %%mm1 \n\t"
70 "punpckldq %1, %%mm0 \n\t"
71 "punpckldq %3, %%mm1 \n\t"
72 "movq %%mm0, %%mm2 \n\t"
73 "pfmul %%mm1, %%mm0 \n\t"
74 "pswapd %%mm1, %%mm1 \n\t"
75 "pfmul %%mm1, %%mm2 \n\t"
76 "pfpnacc %%mm2, %%mm0 \n\t"
77 ::"m"(in2[-2*k]), "m"(in1[2*k]),
78 "m"(tcos[k]), "m"(tsin[k])
79 );
80 asm volatile(
81 "movq %%mm0, %0 \n\t"
82 :"=m"(z[revtab[k]])
83 );
84 }
85
86 ff_fft_calc_3dn2(&s->fft, z);
87
88 /* post rotation + reordering */
89 for(k = 0; k < n4; k++) {
90 asm volatile(
91 "movq %0, %%mm0 \n\t"
92 "movd %1, %%mm1 \n\t"
93 "punpckldq %2, %%mm1 \n\t"
94 "movq %%mm0, %%mm2 \n\t"
95 "pfmul %%mm1, %%mm0 \n\t"
96 "pswapd %%mm1, %%mm1 \n\t"
97 "pfmul %%mm1, %%mm2 \n\t"
98 "pfpnacc %%mm2, %%mm0 \n\t"
99 "movq %%mm0, %0 \n\t"
100 :"+m"(z[k])
101 :"m"(tcos[k]), "m"(tsin[k])
102 );
103 }
104 }
105
106 void ff_imdct_calc_3dn2(MDCTContext *s, FFTSample *output,
107 const FFTSample *input, FFTSample *tmp)
108 {
109 x86_reg k;
110 long n8, n2, n;
111 FFTComplex *z = (FFTComplex *)tmp;
112
113 n = 1 << s->nbits;
114 n2 = n >> 1;
115 n8 = n >> 3;
116
117 imdct_3dn2(s, input, tmp);
118
119 k = n-8;
120 asm volatile("movd %0, %%mm7" ::"r"(1<<31));
121 asm volatile(
122 "1: \n\t"
123 "movq (%4,%0), %%mm0 \n\t" // z[n8+k]
124 "neg %0 \n\t"
125 "pswapd -8(%4,%0), %%mm1 \n\t" // z[n8-1-k]
126 "movq %%mm0, %%mm2 \n\t"
127 "pxor %%mm7, %%mm2 \n\t"
128 "punpckldq %%mm1, %%mm2 \n\t"
129 "pswapd %%mm2, %%mm3 \n\t"
130 "punpckhdq %%mm1, %%mm0 \n\t"
131 "pswapd %%mm0, %%mm4 \n\t"
132 "pxor %%mm7, %%mm0 \n\t"
133 "pxor %%mm7, %%mm4 \n\t"
134 "movq %%mm3, -8(%3,%0) \n\t" // output[n-2-2*k] = { z[n8-1-k].im, -z[n8+k].re }
135 "movq %%mm4, -8(%2,%0) \n\t" // output[n2-2-2*k]= { -z[n8-1-k].re, z[n8+k].im }
136 "neg %0 \n\t"
137 "movq %%mm0, (%1,%0) \n\t" // output[2*k] = { -z[n8+k].im, z[n8-1-k].re }
138 "movq %%mm2, (%2,%0) \n\t" // output[n2+2*k] = { -z[n8+k].re, z[n8-1-k].im }
139 "sub $8, %0 \n\t"
140 "jge 1b \n\t"
141 :"+r"(k)
142 :"r"(output), "r"(output+n2), "r"(output+n), "r"(z+n8)
143 :"memory"
144 );
145 asm volatile("femms");
146 }
147
148 void ff_imdct_half_3dn2(MDCTContext *s, FFTSample *output,
149 const FFTSample *input, FFTSample *tmp)
150 {
151 x86_reg j, k;
152 long n8, n4, n;
153 FFTComplex *z = (FFTComplex *)tmp;
154
155 n = 1 << s->nbits;
156 n4 = n >> 2;
157 n8 = n >> 3;
158
159 imdct_3dn2(s, input, tmp);
160
161 j = -n;
162 k = n-8;
163 asm volatile("movd %0, %%mm7" ::"r"(1<<31));
164 asm volatile(
165 "1: \n\t"
166 "movq (%3,%1), %%mm0 \n\t" // z[n8+k]
167 "pswapd (%3,%0), %%mm1 \n\t" // z[n8-1-k]
168 "movq %%mm0, %%mm2 \n\t"
169 "punpckldq %%mm1, %%mm0 \n\t"
170 "punpckhdq %%mm2, %%mm1 \n\t"
171 "pxor %%mm7, %%mm0 \n\t"
172 "pxor %%mm7, %%mm1 \n\t"
173 "movq %%mm0, (%2,%1) \n\t" // output[n4+2*k] = { -z[n8+k].re, z[n8-1-k].im }
174 "movq %%mm1, (%2,%0) \n\t" // output[n4-2-2*k] = { -z[n8-1-k].re, z[n8+k].im }
175 "sub $8, %1 \n\t"
176 "add $8, %0 \n\t"
177 "jl 1b \n\t"
178 :"+r"(j), "+r"(k)
179 :"r"(output+n4), "r"(z+n8)
180 :"memory"
181 );
182 asm volatile("femms");
183 }
184