d742c0ac467918fc42fdc3284770576c5f53e5df
[strongswan.git] / src / libstrongswan / math / libnttfft / ntt_fft.c
1 /*
2 * Copyright (C) 2014-2016 Andreas Steffen
3 * HSR Hochschule fuer Technik Rapperswil
4 *
5 * This program is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License as published by the
7 * Free Software Foundation; either version 2 of the License, or (at your
8 * option) any later version. See <http://www.fsf.org/copyleft/gpl.txt>.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 * for more details.
14 */
15
16 #include "ntt_fft.h"
17 #include "ntt_fft_reduce.h"
18
19 typedef struct private_ntt_fft_t private_ntt_fft_t;
20
21 /**
22 * Private data structure for ntt_fft_t object
23 */
24 struct private_ntt_fft_t {
25
26 /**
27 * Public interface.
28 */
29 ntt_fft_t public;
30
31 /**
32 * FFT parameter set used as constants
33 */
34 ntt_fft_params_t *p;
35
36 };
37
38 METHOD(ntt_fft_t, get_size, uint16_t,
39 private_ntt_fft_t *this)
40 {
41 return this->p->n;
42 }
43
44 METHOD(ntt_fft_t, get_modulus, uint16_t,
45 private_ntt_fft_t *this)
46 {
47 return this->p->q;
48 }
49
50 /**
51 * Do an FFT butterfly operation
52 *
53 * x[i1] ---|+|------- x[i1]
54 * \/
55 * /\ w[iw]
56 * x[i2] ---|-|--|*|-- x[i2]
57 *
58 */
59 static void butterfly(private_ntt_fft_t *this, uint32_t *x, int i1,int i2, int iw)
60 {
61 uint32_t xp, xm;
62
63 xp = x[i1] + x[i2];
64 xm = x[i1] + (this->p->q - x[i2]);
65 if (xp >= this->p->q)
66 {
67 xp -= this->p->q;
68 }
69 x[i1] = xp;
70 x[i2] = ntt_fft_mreduce(xm * this->p->wr[iw], this->p);
71 }
72
73 /**
74 * Trivial butterfly operation of last FFT stage
75 */
76 static void butterfly_last(private_ntt_fft_t *this, uint32_t *x, int i1)
77 {
78 uint32_t xp, xm;
79 int i2 = i1 + 1;
80
81 xp = x[i1] + x[i2];
82 xm = x[i1] + (this->p->q - x[i2]);
83 if (xp >= this->p->q)
84 {
85 xp -= this->p->q;
86 }
87 if (xm >= this->p->q)
88 {
89 xm -= this->p->q;
90 }
91 x[i1] = xp;
92 x[i2] = xm;
93 }
94
95 METHOD(ntt_fft_t, transform, void,
96 private_ntt_fft_t *this, uint32_t *a, uint32_t *b, bool inverse)
97 {
98 int stage, i, j, k, m, n, s, t, iw, i_rev;
99 uint32_t tmp;
100
101 /* we are going to use the transform size n a lot */
102 n = this->p->n;
103 s = this->p->s;
104
105 if (!inverse)
106 {
107 /* apply linear phase needed for negative wrapped convolution */
108 for (i = 0; i < n; i++)
109 {
110 b[i] = ntt_fft_mreduce(a[i] * this->p->wf[s*i], this->p);
111 }
112 }
113 else if (a != b)
114 {
115 /* copy if input and output array are not the same */
116 for (i = 0; i < n; i++)
117 {
118 b[i] = a[i];
119 }
120 }
121
122 m = n;
123 k = 1;
124
125 for (stage = this->p->stages; stage > 0; stage--)
126 {
127 m >>= 1;
128 t = 0;
129
130 for (j = 0; j < k; j++)
131 {
132 if (stage == 1)
133 {
134 butterfly_last(this, b, t);
135 }
136 else
137 {
138 for (i = 0; i < m; i++)
139 {
140 iw = s * (inverse ? (n - i * k) : (i * k));
141 butterfly(this, b, t + i, t + i + m, iw);
142 }
143 }
144 t += 2*m;
145 }
146 k <<= 1;
147 }
148
149 /* Sort output in bit-reverse order */
150 for (i = 0; i < n; i++)
151 {
152 i_rev = this->p->rev[i];
153
154 if (i_rev > i)
155 {
156 tmp = b[i];
157 b[i] = b[i_rev];
158 b[i_rev] = tmp;
159 }
160 }
161
162 /**
163 * Compensate the linear phase needed for negative wrapped convolution
164 * and normalize the output array with 1/n mod q after the inverse FFT.
165 */
166 if (inverse)
167 {
168 for (i = 0; i < n; i++)
169 {
170 b[i] = ntt_fft_mreduce(b[i] * this->p->wi[i], this->p);
171 }
172 }
173 }
174
175 METHOD(ntt_fft_t, destroy, void,
176 private_ntt_fft_t *this)
177 {
178 free(this);
179 }
180
181 /**
182 * See header.
183 */
184 ntt_fft_t *ntt_fft_create(ntt_fft_params_t *params)
185 {
186 private_ntt_fft_t *this;
187
188 INIT(this,
189 .public = {
190 .get_size = _get_size,
191 .get_modulus = _get_modulus,
192 .transform = _transform,
193 .destroy = _destroy,
194 },
195 .p = params,
196 );
197
198 return &this->public;
199 }