![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ntt_fft.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / strongswan / src / libstrongswan / math / libnttfft|
Return to ntt_fft.c CVS log | Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / strongswan / src / libstrongswan / math / libnttfft |
1.1 misho 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: }