1#include "cplx_fft_internal.h"
 2
 3void cplx_set(CPLX r, const CPLX a) {
 4  r[0] = a[0];
 5  r[1] = a[1];
 6}
 7void cplx_neg(CPLX r, const CPLX a) {
 8  r[0] = -a[0];
 9  r[1] = -a[1];
10}
11void cplx_add(CPLX r, const CPLX a, const CPLX b) {
12  r[0] = a[0] + b[0];
13  r[1] = a[1] + b[1];
14}
15void cplx_sub(CPLX r, const CPLX a, const CPLX b) {
16  r[0] = a[0] - b[0];
17  r[1] = a[1] - b[1];
18}
19void cplx_mul(CPLX r, const CPLX a, const CPLX b) {
20  double re = a[0] * b[0] - a[1] * b[1];
21  r[1] = a[0] * b[1] + a[1] * b[0];
22  r[0] = re;
23}
24
25/**
26 * @brief splits 2h evaluations of one polynomials into 2 times h evaluations of even/odd polynomial
27 * Input:  Q_0(y),...,Q_{h-1}(y),Q_0(-y),...,Q_{h-1}(-y)
28 * Output: P_0(z),...,P_{h-1}(z),P_h(z),...,P_{2h-1}(z)
29 * where Q_i(X)=P_i(X^2)+X.P_{h+i}(X^2) and y^2 = z
30 * @param h number of "coefficients" h >= 1
31 * @param data 2h complex coefficients interleaved and 256b aligned
32 * @param powom y represented as (yre,yim)
33 */
34EXPORT void cplx_split_fft_ref(int32_t h, CPLX* data, const CPLX powom) {
35  CPLX* d0 = data;
36  CPLX* d1 = data + h;
37  for (uint64_t i = 0; i < h; ++i) {
38    CPLX diff;
39    cplx_sub(diff, d0[i], d1[i]);
40    cplx_add(d0[i], d0[i], d1[i]);
41    cplx_mul(d1[i], diff, powom);
42  }
43}
44
45/**
46 * @brief Do two layers of itwiddle (i.e. split).
47 * Input/output:  d0,d1,d2,d3 of length h
48 * Algo:
49 *   itwiddle(d0,d1,om[0]),itwiddle(d2,d3,i.om[0])
50 *   itwiddle(d0,d2,om[1]),itwiddle(d1,d3,om[1])
51 * @param h number of "coefficients" h >= 1
52 * @param data 4h complex coefficients interleaved and 256b aligned
53 * @param powom om[0] (re,im) and om[1] where om[1]=om[0]^2
54 */
55EXPORT void cplx_bisplit_fft_ref(int32_t h, CPLX* data, const CPLX powom[2]) {
56  CPLX* d0 = data;
57  CPLX* d2 = data + 2 * h;
58  const CPLX* om0 = powom;
59  CPLX iom0;
60  iom0[0] = powom[0][1];
61  iom0[1] = -powom[0][0];
62  const CPLX* om1 = powom + 1;
63  cplx_split_fft_ref(h, d0, *om0);
64  cplx_split_fft_ref(h, d2, iom0);
65  cplx_split_fft_ref(2 * h, d0, *om1);
66}
67
68/**
69 * Input: Q(y),Q(-y)
70 * Output: P_0(z),P_1(z)
71 * where Q(X)=P_0(X^2)+X.P_1(X^2) and y^2 = z
72 * @param data 2 complexes coefficients interleaved and 256b aligned
73 * @param powom (z,-z) interleaved: (zre,zim,-zre,-zim)
74 */
75void split_fft_last_ref(CPLX* data, const CPLX powom) {
76  CPLX diff;
77  cplx_sub(diff, data[0], data[1]);
78  cplx_add(data[0], data[0], data[1]);
79  cplx_mul(data[1], diff, powom);
80}