Annotation of embedaddon/strongswan/src/libstrongswan/math/libnttfft/ntt_fft.c, revision 1.1
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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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