integrity-test: check code and ro segments of libnttfft
[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 /**
20 * Described in header.
21 */
22 void libnttfft_init(void)
23 {
24 /* empty */
25 }
26
27 typedef struct private_ntt_fft_t private_ntt_fft_t;
28
29 /**
30 * Private data structure for ntt_fft_t object
31 */
32 struct private_ntt_fft_t {
33
34 /**
35 * Public interface.
36 */
37 ntt_fft_t public;
38
39 /**
40 * FFT parameter set used as constants
41 */
42 const ntt_fft_params_t *p;
43
44 };
45
46 METHOD(ntt_fft_t, get_size, uint16_t,
47 private_ntt_fft_t *this)
48 {
49 return this->p->n;
50 }
51
52 METHOD(ntt_fft_t, get_modulus, uint16_t,
53 private_ntt_fft_t *this)
54 {
55 return this->p->q;
56 }
57
58 /**
59 * Do an FFT butterfly operation
60 *
61 * x[i1] ---|+|------- x[i1]
62 * \/
63 * /\ w[iw]
64 * x[i2] ---|-|--|*|-- x[i2]
65 *
66 */
67 static void butterfly(private_ntt_fft_t *this, uint32_t *x, int i1,int i2, int iw)
68 {
69 uint32_t xp, xm;
70
71 xp = x[i1] + x[i2];
72 xm = x[i1] + (this->p->q - x[i2]);
73 if (xp >= this->p->q)
74 {
75 xp -= this->p->q;
76 }
77 x[i1] = xp;
78 x[i2] = ntt_fft_mreduce(xm * this->p->wr[iw], this->p);
79 }
80
81 /**
82 * Trivial butterfly operation of last FFT stage
83 */
84 static void butterfly_last(private_ntt_fft_t *this, uint32_t *x, int i1)
85 {
86 uint32_t xp, xm;
87 int i2 = i1 + 1;
88
89 xp = x[i1] + x[i2];
90 xm = x[i1] + (this->p->q - x[i2]);
91 if (xp >= this->p->q)
92 {
93 xp -= this->p->q;
94 }
95 if (xm >= this->p->q)
96 {
97 xm -= this->p->q;
98 }
99 x[i1] = xp;
100 x[i2] = xm;
101 }
102
103 METHOD(ntt_fft_t, transform, void,
104 private_ntt_fft_t *this, uint32_t *a, uint32_t *b, bool inverse)
105 {
106 int stage, i, j, k, m, n, s, t, iw, i_rev;
107 uint32_t tmp;
108
109 /* we are going to use the transform size n a lot */
110 n = this->p->n;
111 s = this->p->s;
112
113 if (!inverse)
114 {
115 /* apply linear phase needed for negative wrapped convolution */
116 for (i = 0; i < n; i++)
117 {
118 b[i] = ntt_fft_mreduce(a[i] * this->p->wf[s*i], this->p);
119 }
120 }
121 else if (a != b)
122 {
123 /* copy if input and output array are not the same */
124 for (i = 0; i < n; i++)
125 {
126 b[i] = a[i];
127 }
128 }
129
130 m = n;
131 k = 1;
132
133 for (stage = this->p->stages; stage > 0; stage--)
134 {
135 m >>= 1;
136 t = 0;
137
138 for (j = 0; j < k; j++)
139 {
140 if (stage == 1)
141 {
142 butterfly_last(this, b, t);
143 }
144 else
145 {
146 for (i = 0; i < m; i++)
147 {
148 iw = s * (inverse ? (n - i * k) : (i * k));
149 butterfly(this, b, t + i, t + i + m, iw);
150 }
151 }
152 t += 2*m;
153 }
154 k <<= 1;
155 }
156
157 /* Sort output in bit-reverse order */
158 for (i = 0; i < n; i++)
159 {
160 i_rev = this->p->rev[i];
161
162 if (i_rev > i)
163 {
164 tmp = b[i];
165 b[i] = b[i_rev];
166 b[i_rev] = tmp;
167 }
168 }
169
170 /**
171 * Compensate the linear phase needed for negative wrapped convolution
172 * and normalize the output array with 1/n mod q after the inverse FFT.
173 */
174 if (inverse)
175 {
176 for (i = 0; i < n; i++)
177 {
178 b[i] = ntt_fft_mreduce(b[i] * this->p->wi[i], this->p);
179 }
180 }
181 }
182
183 METHOD(ntt_fft_t, destroy, void,
184 private_ntt_fft_t *this)
185 {
186 free(this);
187 }
188
189 /**
190 * See header.
191 */
192 ntt_fft_t *ntt_fft_create(const ntt_fft_params_t *params)
193 {
194 private_ntt_fft_t *this;
195
196 INIT(this,
197 .public = {
198 .get_size = _get_size,
199 .get_modulus = _get_modulus,
200 .transform = _transform,
201 .destroy = _destroy,
202 },
203 .p = params,
204 );
205
206 return &this->public;
207 }