Implemented Number Theoretic Transform using the FFT algorithm
authorAndreas Steffen <andreas.steffen@strongswan.org>
Sat, 25 Oct 2014 17:09:41 +0000 (19:09 +0200)
committerAndreas Steffen <andreas.steffen@strongswan.org>
Sat, 29 Nov 2014 13:51:14 +0000 (14:51 +0100)
By pre-multiplying the input arrays with a linear phase the
fast multiplication via FFT and inverse FFT computes a negative
wrapped convolution corresponding to a modulus of x^n+1.

src/libstrongswan/plugins/bliss/Makefile.am
src/libstrongswan/plugins/bliss/bliss_fft.c [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_fft.h [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_fft_params.c [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_fft_params.h [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_private_key.c
src/libstrongswan/plugins/bliss/bliss_private_key.h
src/libstrongswan/plugins/bliss/bliss_public_key.h

index 3075688..94eb6da 100644 (file)
@@ -14,6 +14,8 @@ endif
 libstrongswan_bliss_la_SOURCES = \
        bliss_plugin.h bliss_plugin.c \
        bliss_private_key.h bliss_private_key.c \
-       bliss_public_key.h bliss_public_key.c
+       bliss_public_key.h bliss_public_key.c \
+       bliss_fft.h bliss_fft.c \
+       bliss_fft_params.h bliss_fft_params.c
 
 libstrongswan_bliss_la_LDFLAGS = -module -avoid-version
diff --git a/src/libstrongswan/plugins/bliss/bliss_fft.c b/src/libstrongswan/plugins/bliss/bliss_fft.c
new file mode 100644 (file)
index 0000000..033c214
--- /dev/null
@@ -0,0 +1,199 @@
+/*
+ * Copyright (C) 2014 Andreas Steffen
+ * HSR Hochschule fuer Technik Rapperswil
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2 of the License, or (at your
+ * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+#include "bliss_fft.h"
+
+typedef struct private_bliss_fft_t private_bliss_fft_t;
+
+/**
+ * Private data structure for bliss_fft_t object
+ */
+struct private_bliss_fft_t {
+       /**
+        * Public interface.
+        */
+       bliss_fft_t public;
+
+       /**
+        * FFT parameter set used as constants
+        */
+       bliss_fft_params_t *p;
+
+};
+
+METHOD(bliss_fft_t, get_size, uint16_t,
+       private_bliss_fft_t *this)
+{
+       return this->p->n;
+}
+
+METHOD(bliss_fft_t, get_modulus, uint16_t,
+       private_bliss_fft_t *this)
+{
+       return this->p->q;
+}
+
+/**
+ * Do an FFT butterfly operation
+ *
+ * x[i1] ---|+|------- x[i1]
+ *        \/
+ *        /\    w[iw]  
+ * x[i2] ---|-|--|*|-- x[i2]
+ *
+ */
+static void butterfly(private_bliss_fft_t *this, uint32_t *x, int i1,int i2,
+                                                                                                                         int iw)
+{
+       uint32_t xp, xm;
+
+       xp = x[i1] + x[i2];
+       xm = x[i1] + (this->p->q - x[i2]);
+       if (xp >= this->p->q)
+       {
+               xp -= this->p->q;
+       }
+       x[i1] =  xp;
+       x[i2] = (xm * this->p->w[iw]) % this->p->q;
+}
+
+/**
+ * Trivial butterfly operation of last FFT stage
+ */
+static void butterfly_last(private_bliss_fft_t *this, uint32_t *x, int i1)
+{
+       uint32_t xp, xm;
+       int i2 = i1 + 1;
+
+       xp = x[i1] + x[i2];
+       xm = x[i1] + (this->p->q - x[i2]);
+       if (xp >= this->p->q)
+       {
+               xp -= this->p->q;
+       }
+       if (xm >= this->p->q)
+       {
+               xm -= this->p->q;
+       }
+       x[i1] = xp;
+       x[i2] = xm;
+}
+
+METHOD(bliss_fft_t, transform, void,
+       private_bliss_fft_t *this, uint32_t *a, uint32_t *b, bool inverse)
+{
+       int stage, i, j, k, m, n, t, iw, i_rev;
+       uint16_t q;
+       uint32_t tmp;
+
+       /* we are going to use the transform size n and the modulus q a lot */
+       n = this->p->n;
+       q = this->p->q;
+
+       if (!inverse)
+       {
+               /* apply linear phase needed for negative wrapped convolution */
+               for (i = 0; i < n; i++)
+               {
+                       b[i] = (a[i] * this->p->w[i]) % q;
+               }
+       }
+       else if (a != b)
+       {
+               /* copy if input and output array are not the same */
+               for (i = 0; i < n; i++)
+               {
+                       b[i] = a[i];
+               }
+       }
+
+       m = n;
+       k = 1;
+
+       for (stage = this->p->stages; stage > 0; stage--)
+       {
+               m >>= 1;
+               t = 0;
+
+               for (j = 0; j < k; j++)
+               {
+                       if (stage == 1)
+                       {
+                               butterfly_last(this, b, t);
+                       }
+                       else
+                       {
+                               for (i = 0; i < m; i++)
+                               {
+                                       iw = 2 * (inverse ? (n - i * k) : (i * k));
+                                       butterfly(this, b, t + i, t + i + m, iw);
+                               }                               
+                       }
+                       t += 2*m;
+               }
+               k <<= 1;
+       }
+
+       /* Sort output in bit-reverse order */
+       for (i = 0; i < n; i++)
+       {
+               i_rev = this->p->rev[i];
+
+               if (i_rev > i)
+               {
+                       tmp = b[i];
+                       b[i] = b[i_rev];
+                       b[i_rev] = tmp;
+               }
+       }
+
+       /**
+        * Compensate the linear phase needed for negative wrapped convolution
+        * and normalize the output array with 1/n mod q after the inverse FFT. 
+        */
+       if (inverse)
+       {
+               for (i = 0; i < n; i++)
+               {
+                       b[i] = (((b[i] * this->p->w[2*n - i]) % q) * this->p->n_inv) % q;
+               }
+       }
+}
+
+METHOD(bliss_fft_t, destroy, void,
+       private_bliss_fft_t *this)
+{
+       free(this);
+}
+
+/**
+ * See header.
+ */
+bliss_fft_t *bliss_fft_create(bliss_fft_params_t *params)
+{
+       private_bliss_fft_t *this;
+
+       INIT(this,
+               .public = {
+                       .get_size = _get_size,
+                       .get_modulus = _get_modulus,
+                       .transform = _transform,
+                       .destroy = _destroy,
+               },
+               .p = params,
+       );
+
+       return &this->public;
+}
diff --git a/src/libstrongswan/plugins/bliss/bliss_fft.h b/src/libstrongswan/plugins/bliss/bliss_fft.h
new file mode 100644 (file)
index 0000000..a79edd2
--- /dev/null
@@ -0,0 +1,71 @@
+/*
+ * Copyright (C) 2014 Andreas Steffen
+ * HSR Hochschule fuer Technik Rapperswil
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2 of the License, or (at your
+ * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+/**
+ * @defgroup bliss_fft bliss_fft
+ * @{ @ingroup bliss_p
+ */
+
+#ifndef BLISS_FFT_H_
+#define BLISS_FFT_H_
+
+#include "bliss_fft_params.h"
+
+#include <library.h>
+
+typedef struct bliss_fft_t bliss_fft_t;
+
+/**
+ * Implements a Number Theoretic Transform (NTT) via the FFT algorithm
+ */
+struct bliss_fft_t {
+
+       /**
+        * Get the size of the Number Theoretic Transform
+        *
+        * @result                      Transform size
+        */
+       uint16_t (*get_size)(bliss_fft_t *this);
+
+       /**
+        * Get the prime modulus of the Number Theoretic Transform
+        *
+        * @result                      Prime modulus
+        */
+       uint16_t (*get_modulus)(bliss_fft_t *this);
+
+       /**
+        * Compute the [inverse] NTT of a polynomial
+        *
+        * @param a                     Coefficient of input polynomial
+        * @param b                     Coefficient of output polynomial
+        * @param inverse       TRUE if the inverse NTT has to be computed
+        */
+       void (*transform)(bliss_fft_t *this, uint32_t *a, uint32_t *b, bool inverse);
+
+       /**
+        * Destroy bliss_fft_t object
+        */
+       void (*destroy)(bliss_fft_t *this);
+};
+
+/**
+ * Create a bliss_fft_t object for a given FFT parameter set
+ *
+ * @param params               FFT parameters
+ */
+bliss_fft_t *bliss_fft_create(bliss_fft_params_t *params);
+
+#endif /** BLISS_FFT_H_ @}*/
diff --git a/src/libstrongswan/plugins/bliss/bliss_fft_params.c b/src/libstrongswan/plugins/bliss/bliss_fft_params.c
new file mode 100644 (file)
index 0000000..c892c06
--- /dev/null
@@ -0,0 +1,215 @@
+/*
+ * Copyright (C) 2014 Andreas Steffen
+ * HSR Hochschule fuer Technik Rapperswil
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2 of the License, or (at your
+ * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+#include "bliss_fft_params.h"
+
+/**
+ * FFT parameters for q = 12289 and 2n = 1024
+ */
+static uint16_t w_12289_1024[] = {
+           1,    49,  2401,  7048,  1260,   295,  2166,  7822,  2319,  3030,
+        1002, 12231,  9447,  8210,  9042,   654,  7468,  9551,  1017,   677,
+        8595,  3329,  3364,  5079,  3091,  3991, 11224,  9260, 11336,  2459,
+        9890,  5339,  3542,  1512,   354,  5057,  2013,   325,  3636,  6118,
+        4846,  3963,  9852,  3477, 10616,  4046,  1630,  6136,  5728, 10314,
+        1537,  1579,  3637,  6167,  7247, 11011, 11112,  3772,   493, 11868,
+        3949,  9166,  6730, 10256, 10984,  9789,   390,  6821,  2426,  8273,
+       12129,  4449,  9088,  2908,  7313,  1956,  9821,  1958,  9919,  6760,
+       11726,  9280,    27,  1323,  3382,  5961,  9442,  7965,  9326,  2281,
+        1168,  8076,  2476, 10723,  9289,   468, 10643,  5369,  5012, 12097,
+
+        2881,  5990, 10863,  3860,  4805,  1954,  9723,  9445,  8112,  4240,
+       11136,  4948,  8961,  8974,  9611,  3957,  9558,  1360,  5195,  8775,
+       12149,  5429,  7952,  8689,  7935,  7856,  3985, 10930,  7143,  5915,
+        7188,  8120,  4632,  5766, 12176,  6752, 11334,  2361,  5088,  3532,
+        1022,   922,  8311,  1702,  9664,  6554,  1632,  6234, 10530, 12121,
+        4057,  2169,  7969,  9522, 11885,  4782,   827,  3656,  7098,  3710,
+        9744, 10474,  9377,  4780,   729, 11143,  5291,  1190,  9154,  6142,
+        6022,   142,  6958,  9139,  5407,  6874,  5023,   347,  4714,  9784,
+         145,  7105,  4053,  1973, 10654,  5908,  6845,  3602,  4452,  9235,
+       10111,  3879,  5736, 10706,  8456,  8807,  1428,  8527, 12286, 12142,
+
+        5086,  3434,  8509, 11404,  5791,  1112,  5332,  3199,  9283,   174,
+        8526, 12237,  9741, 10327,  2174,  8214,  9238, 10258, 11082,  2302,
+        2197,  9341,  3016,   316,  3195,  9087,  2859,  4912,  7197,  8561,
+        1663,  7753, 11227,  9407,  6250, 11314,  1381,  6224, 10040,   400,
+     7311,  1858,  5019,   151,  7399,  6170,  7394,  5925,  7678,  7552,
+        1378,  6077,  2837,  3834,  3531,   973, 10810,  1263,   442,  9369,
+        4388,  6099,  3915,  7500, 11119,  4115,  5011, 12048,   480, 11231,
+        9603,  3565,  2639,  6421,  7404,  6415,  7110,  4298,  1689,  9027,
+       12208,  8320,  2143,  6695,  8541,   683,  8889,  5446,  8785,   350,
+        4861,  4698,  9000, 10885,  4938,  8471,  9542,   576,  3646,  6608,
+
+        4278,   709, 10163,  6427,  7698,  8532,   242, 11858,  3459,  9734,
+        9984,  9945,  8034,   418,  8193,  8209,  8993, 10542,   420,  8291,
+         722, 10800,   773,  1010,   334,  4077,  3149,  6833,  3014,   218,
+       10682,  7280,   339,  4322,  2865,  5206,  9314,  1693,  9223,  9523,
+       11934,  7183,  7875,  4916,  7393,  5876,  5277,   504,   118,  5782,
+         671,  8301,  1212, 10232,  9808,  1321,  3284,  1159,  7635,  5445,
+        8736, 10238, 10102,  3438,  8705,  8719,  9405,  6152,  6512, 11863,
+        3704,  9450,  8357,  3956,  9509, 11248, 10436,  7515, 11854,  3263,
+         130,  6370,  4905,  6854,  4043,  1483, 11222,  9162,  6534,   652,
+        7370,  4749, 11499, 10446,  8005, 11286,     9,   441,  9320,  1987,
+
+       11340,  2655,  7205,  8953,  8582,  2692,  9018, 11767, 11289,   156,
+        7644,  5886,  5767, 12225,  9153,  6093,  3621,  5383,  5698,  8844,
+        3241, 11341,  2704,  9606,  3712,  9842,  2987, 11184,  7300,  1319,
+        3186,  8646,  5828,  2925,  8146,  5906,  6747, 11089,  2645,  6715,
+        9521, 11836,  2381,  6068,  2396,  6803,  1544,  1922,  8155,  6347,
+        3778,   787,  1696,  9370,  4437,  8500, 10963,  8760, 11414,  6281,
+         544,  2078,  3510, 12233,  9545,   723, 10849,  3174,  8058,  1594,
+        4372,  5315,  2366,  5333,  3248, 11684,  7222,  9786,   243, 11907,
+        5860,  4493, 11244, 10240, 10200,  8240, 10512, 11239,  9995, 10484,
+        9867,  4212,  9764, 11454,  8241, 10561,  1351,  4754, 11744, 10162,
+
+        6378,  5297,  1484, 11271, 11563,  1293,  1912,  7665,  6915,  7032,
+         476, 11035, 12288, 12240,  9888,  5241, 11029, 11994, 10123,  4467,
+        9970,  9259, 11287,    58,  2842,  4079,  3247, 11635,  4821,  2738,
+       11272, 11612,  3694,  8960,  8925,  7210,  9198,  8298,  1065,  3029,
+         953,  9830,  2399,  6950,  8747, 10777, 11935,  7232, 10276, 11964,
+        8653,  6171,  7443,  8326,  2437,  8812,  1673,  8243, 10659,  6153,
+        6561,  1975, 10752, 10710,  8652,  6122,  5042,  1278,  1177,  8517,
+       11796,   421,  8340,  3123,  5559,  2033,  1305,  2500, 11899,  5468,
+        9863,  4016,   160,  7840,  3201,  9381,  4976, 10333,  2468, 10331,
+        2370,  5529,  563,   3009, 12262, 10966,  8907,  6328,  2847,  4324,
+
+        2963, 10008, 11121,  4213,  9813,  1566,  3000, 11821,  1646,  6920,
+        7277,   192,  9408,  6299,  1426,  8429,  7484, 10335,  2566,  2844,
+        4177,  8049,  1153,  7341,  3328,  3315,  2678,  8332,  2731, 10929,
+        7094,  3514,   140,  6860,  4337,  3600,  4354,  4433,  8304,  1359,
+        5146,  6374,  5101,  4169,  7657,  6523,   113,  5537,   955,  9928,
+        7201,  8757, 11267, 11367,  3978, 10587,  2625,  5735, 10657,  6055,
+        1759,   168,  8232, 10120,  4320,  2767,   404,  7507, 11462,  8633,
+        5191,  8579,  2545,  1815,  2912,  7509, 11560,  1146,  6998, 11099,
+        3135,  6147,  6267, 12147,  5331,  3150,  6882,  5415,  7266, 11942,
+        7575,  2505, 12144,  5184,  8236, 10316,  1635,  6381,  5444,  8687,
+
+        7837,  3054,  2178,  8410,  6553,  1583,  3833,  3482, 10861,  3762,
+           3,   147,  7203,  8855,  3780,   885,  6498, 11177,  6957,  9090,
+        3006, 12115,  3763,    52,  2548,  1962, 10115,  4075,  3051,  2031,
+        1207,  9987, 10092,  2948,  9273, 11973,  9094,  3202,  9430,  7377,
+        5092,  3728, 10626,  4536,  1062,  2882,  6039,   975, 10908,  6065,
+        2249, 11889,  4978, 10431,  7270, 12138,  4890,  6119,  4895,  6364,
+        4611,  4737, 10911,  6212,  9452,  8455,  8758, 11316,  1479, 11026,
+       11847,  2920,  7901,  6190,  8374,  4789,  1170,  8174,  7278,   241,
+       11809,  1058,  2686,  8724,  9650,  5868,  4885,  5874,  5179,  7991,
+       10600,  3262,    81,  3969, 10146,  5594,  3748, 11606,  3400,  6843,
+
+        3504, 11939,  7428,  7591,  3289,  1404,  7351,  3818,  2747, 11713,
+        8643,  5681,  8011, 11580,  2126,  5862,  4591,  3757, 12047,   431,
+        8830,  2555,  2305,  2344,  4255, 11871,  4096,  4080,  3296,  1747,
+       11869,  3998, 11567,  1489, 11516, 11279, 11955,  8212,  9140,  5456,
+        9275, 12071,  1607,  5009, 11950,  7967,  9424,  7083,  2975, 10596,
+        3066,  2766,   355,  5106,  4414,  7373,  4896,  6413,  7012, 11785,
+       12171,  6507, 11618,  3988, 11077,  2057,  2481, 10968,  9005, 11130,
+        4654,  6844,  3553,  2051,  2187,  8851,  3584,  3570,  2884,  6137,
+        5777,   426,  8585,  2839,  3932,  8333,  2780,  1041,  1853,  4774,
+         435,  9026, 12159,  5919,  7384,  5435,  8246, 10806,  1067,  3127,
+
+        5755, 11637,  4919,  7540,   790,  1843,  4284,  1003, 12280, 11848,
+        2969, 10302,   949,  9634,  5084,  3336,  3707,  9597,  3271,   522,
+        1000, 12133,  4645,  6403,  6522,    64,  3136,  6196,  8668,  6906,
+        6591,  3445,  9048,   948,  9585,  2683,  8577,  2447,  9302,  1105,
+        4989, 10970,  9103,  3643,  6461,  9364,  4143,  6383,  5542,  1200,
+        9644,  5574,  2768,   453,  9908,  6221,  9893,  5486, 10745, 10367,
+        4134,  5942,  8511, 11502, 10593,  2919,  7852,  3789,  1326,  3529,
+         875,  6008, 11745, 10211,  8779,    56,  2744, 11566,  1440,  9115,
+        4231, 10695,  7917,  6974,  9923,  6956,  9041,   605,  5067,  2503,
+       12046,   382,  6429,  7796,  1045,  2049,  2089,  4049,  1777,  1050,
+
+        2294,  1805,  2422,  8077,  2525,   835,  4048,  1728, 10938,  7535,
+         545,  2127,  5911,  6992, 10805,  1018,   726, 10996, 10377,  4624,
+         5374, 5257, 11813,  1254,     1
+};
+
+/**
+ * Bit-reversed indices for n = 512
+ */
+static uint16_t rev_512[] = {
+         0, 256, 128, 384,  64, 320, 192, 448,  32, 288, 
+       160, 416,  96, 352, 224, 480,  16, 272, 144, 400,
+        80, 336, 208, 464,  48, 304, 176, 432, 112, 368,
+       240, 496,   8, 264, 136, 392,  72, 328, 200, 456,
+        40, 296, 168, 424, 104, 360, 232, 488,  24, 280,
+       152, 408,  88, 344, 216, 472,  56, 312, 184, 440,
+       120, 376, 248, 504,   4, 260, 132, 388,  68, 324,
+       196, 452,  36, 292, 164, 420, 100, 356, 228, 484,
+        20, 276, 148, 404,  84, 340, 212, 468,  52, 308,
+       180, 436, 116, 372, 244, 500,  12, 268, 140, 396,
+
+        76, 332, 204, 460,  44, 300, 172, 428, 108, 364,
+       236, 492,  28, 284, 156, 412,  92, 348, 220, 476,
+        60, 316, 188, 444, 124, 380, 252, 508,   2, 258,
+       130, 386,  66, 322, 194, 450,  34, 290, 162, 418,
+        98, 354, 226, 482,  18, 274, 146, 402,  82, 338,
+       210, 466,  50, 306, 178, 434, 114, 370, 242, 498,
+        10, 266, 138, 394,  74, 330, 202, 458,  42, 298,
+       170, 426, 106, 362, 234, 490,  26, 282, 154, 410,
+        90, 346, 218, 474,  58, 314, 186, 442, 122, 378,
+       250, 506,   6, 262, 134, 390,  70, 326, 198, 454,
+
+        38, 294, 166, 422, 102, 358, 230, 486,  22, 278,
+       150, 406,  86, 342, 214, 470,  54, 310, 182, 438,
+       118, 374, 246, 502,  14, 270, 142, 398,  78, 334,
+       206, 462,  46, 302, 174, 430, 110, 366, 238, 494,
+        30, 286, 158, 414,  94, 350, 222, 478,  62, 318,
+       190, 446, 126, 382, 254, 510,   1, 257, 129, 385,
+        65, 321, 193, 449,  33, 289, 161, 417,  97, 353,
+       225, 481,  17, 273, 145, 401,  81, 337, 209, 465,
+        49, 305, 177, 433, 113, 369, 241, 497,   9, 265,
+       137, 393,  73, 329, 201, 457,  41, 297, 169, 425,
+
+       105, 361, 233, 489,  25, 281, 153, 409,  89, 345,
+       217, 473,  57, 313, 185, 441, 121, 377, 249, 505,
+         5, 261, 133, 389,  69, 325, 197, 453,  37, 293,
+       165, 421, 101, 357, 229, 485,  21, 277, 149, 405,
+        85, 341, 213, 469,  53, 309, 181, 437, 117, 373,
+       245, 501,  13, 269, 141, 397,  77, 333, 205, 461,
+        45, 301, 173, 429, 109, 365, 237, 493,  29, 285,
+       157, 413,  93, 349, 221, 477,  61, 317, 189, 445,
+       125, 381, 253, 509,   3, 259, 131, 387,  67, 323,
+       195, 451,  35, 291, 163, 419,  99, 355, 227, 483,
+
+        19, 275, 147, 403,  83, 339, 211, 467,  51, 307,
+       179, 435, 115, 371, 243, 499,  11, 267, 139, 395,
+        75, 331, 203, 459,  43, 299, 171, 427, 107, 363,
+       235, 491,  27, 283, 155, 411,  91, 347, 219, 475,
+        59, 315, 187, 443, 123, 379, 251, 507,   7, 263,
+       135, 391,  71, 327, 199, 455,  39, 295, 167, 423,
+       103, 359, 231, 487,  23, 279, 151, 407,  87, 343,
+       215, 471,  55, 311, 183, 439, 119, 375, 247, 503,
+        15, 271, 143, 399,  79, 335, 207, 463,  47, 303,
+       175, 431, 111, 367, 239, 495,  31, 287, 159, 415,
+
+        95, 351, 223, 479,  63, 319, 191, 447, 127, 383,
+       255, 511
+};
+
+bliss_fft_params_t bliss_fft_12289_512 = {
+       12289, 512, 12265, 9, w_12289_1024, rev_512
+};
+
+/**
+ * FFT parameters for q = 17 and n = 16
+ */
+static uint16_t w_17_16[] = {
+       1, 3, 9, 10, 13, 5, 15, 11, 16, 14, 8, 7, 4, 12, 2, 6, 1 };
+
+/**
+ * Bit-reversed indices for n = 8
+ */
+static uint16_t rev_8[] = { 0, 4, 2, 6, 1, 5, 3, 7 };
+
+bliss_fft_params_t bliss_fft_17_8 = { 17, 8, 15, 3, w_17_16, rev_8 };
diff --git a/src/libstrongswan/plugins/bliss/bliss_fft_params.h b/src/libstrongswan/plugins/bliss/bliss_fft_params.h
new file mode 100644 (file)
index 0000000..31b151b
--- /dev/null
@@ -0,0 +1,75 @@
+/*
+ * Copyright (C) 2014 Andreas Steffen
+ * HSR Hochschule fuer Technik Rapperswil
+ *
+ * This program is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2 of the License, or (at your
+ * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+/**
+ * @defgroup bliss_fft_params bliss_fft_params
+ * @{ @ingroup bliss_p
+ */
+
+#ifndef BLISS_FFT_PARAMS_H_
+#define BLISS_FFT_PARAMS_H_
+
+#include <library.h>
+
+typedef struct bliss_fft_params_t bliss_fft_params_t;
+
+/**
+ * Defines the parameters for an NTT computed via the FFT algorithm
+ */
+struct bliss_fft_params_t {
+
+       /**
+        * Prime modulus
+        */
+       uint16_t q;
+
+       /**
+        * Size of the FFT with the condition k * n = q-1
+        */
+       uint16_t n;
+
+       /**
+        * Inverse of n mod q used for normalization of the FFT
+        */
+       uint16_t n_inv;
+
+       /**
+        * Number of FFT stages  stages = log2(n)
+        */
+       uint16_t stages;
+
+       /**
+        * FFT twiddle factors (n-th roots of unity)
+        */
+       uint16_t *w;
+
+       /**
+        * FFT bit reversal
+        */
+       uint16_t *rev;
+
+};
+
+/**
+ * FFT parameters for q = 12289 and n = 512
+ */
+extern bliss_fft_params_t bliss_fft_12289_512;
+
+/**
+ * FFT parameters for q = 17 and n = 8
+ */
+extern bliss_fft_params_t bliss_fft_17_8;
+
+#endif /** BLISS_FFT_PARAMS_H_ @}*/
index 7d3cee9..dcf1b7d 100644 (file)
  */
 
 #include "bliss_private_key.h"
+#include "bliss_fft.h"
+
+#define _GNU_SOURCE
+#include <stdlib.h>
 
 typedef struct private_bliss_private_key_t private_bliss_private_key_t;
 
@@ -149,12 +153,257 @@ static private_bliss_private_key_t *bliss_private_key_create_empty(void)
 }
 
 /**
+ * Compute the scalar product of a vector x with a negative wrapped vector y
+ */
+static int16_t wrapped_product(int16_t *x, int16_t *y, int n, int shift)
+{
+       int16_t product = 0;
+       int i;
+
+       for (i = 0; i < n - shift; i++)
+       {
+               product += x[i] * y[i + shift];
+       }
+       for (i = n - shift; i < n; i++)
+       {
+               product -= x[i] * y[i + shift - n];
+       }
+       return product;
+}
+
+/**
+ * Apply a negative wrapped rotation to a vector x
+ */
+static void wrap(int16_t *x, int n, int shift, int16_t *x_wrapped)
+{
+       int i;
+
+       for (i = 0; i < n - shift; i++)
+       {
+               x_wrapped[i + shift] = x[i];
+       }
+       for (i = n - shift; i < n; i++)
+       {
+               x_wrapped[i + shift - n] = -x[i];
+       }
+}
+
+static int compare(const int16_t *a, const int16_t *b)
+{
+       int16_t temp = *a - *b;
+
+       if (temp > 0)
+       {
+               return 1;
+       }
+       else if (temp < 0)
+       {
+               return -1;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+/**
+ * Compute the Nk(S) norm of S = (s1, s2)
+ */
+static uint32_t nks_norm(int16_t *s1, int16_t *s2, int n)
+{
+       int16_t t[n], t_wrapped[n], max_kappa[n];
+       uint32_t nks = 0;
+       int i, j, kappa = 23;
+
+       for (i = 0; i < n; i++)
+       {
+               t[i] = wrapped_product(s1, s1, n, i) + wrapped_product(s2, s2, n, i);
+               DBG1(DBG_LIB, "t[%d] = %5d", i, t[i]);
+       }
+
+       for (i = 0; i < n; i++)
+       {
+               wrap(t, n, i, t_wrapped);
+               qsort(t_wrapped, n, sizeof(int16_t), (__compar_fn_t)compare);
+               max_kappa[i] = 0;
+
+               for (j = 1; j <= kappa; j++)
+               {
+                       max_kappa[i] += t_wrapped[n - j];
+               }
+               DBG1(DBG_LIB, "max_kappa[%d] = %5d", i, max_kappa[i]);
+       }
+       qsort(max_kappa, n, sizeof(int16_t), (__compar_fn_t)compare);
+
+       for (i = 1; i <= kappa; i++)
+       {
+               nks += max_kappa[n - i];
+       }
+       return nks;
+}
+
+/**
+ * Compute the inverse x1 of x modulo q as x^(-1) = x^(q-2) mod q
+ */
+static uint32_t invert(uint32_t x, uint16_t q)
+{
+       uint32_t x1, x2;
+       uint16_t q2;
+       int i, i_max;
+
+       q2 = q - 2;
+       x1 = (q2 & 1) ? x : 1;
+       x2 = x;
+       i_max = 15;
+
+       while ((q2 & (1 << i_max)) == 0)
+       {
+               i_max--;
+       }
+       for (i = 1; i <= i_max; i++)
+       {
+               x2 = (x2 * x2) % q;
+
+               if (q2 & (1 << i))
+               {
+                       x1 = (x1 * x2) % q;
+               }
+       }
+
+       return x1;
+}
+
+/**
  * See header.
  */
 bliss_private_key_t *bliss_private_key_gen(key_type_t type, va_list args)
 {
        private_bliss_private_key_t *this;
        u_int key_size = 1;
+       int i;
+       uint32_t *a, *A, *F, *G, nks;
+       uint16_t q, n, l2_norm;
+       bliss_fft_t *fft;
+
+       int16_t f[] = {
+                0,  0,  0,  0,  1,  1,  0, -1,  0,  1, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+                0,  0, -1,  0,  0,  0, -1,  1,  0,  0, 
+                1,  1, -1,  0,  1,  1,  0,  0,  0,  0, 
+                0,  0,  0,  0,  1,  0, -1,  0, -1,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+               -1,  0,  1,  0,  1,  0,  0,  0,  0,  0, 
+                0,  0,  1,  0,  0,  1,  0, -1,  0,  0, 
+                0, -1,  0,  0,  1, -1,  0,  0,  0,  0, 
+                0,  0,  0,  0,  1,  0,  0,  0,  0,  0, 
+
+                0,  0,  0,  0,  0,  1,  0, -1,  0,  1, 
+                1, -1,  0,  1,  0,  1,  0,  0,  1,  0, 
+               -1,  0,  0,  0,  1, -1,  1,  0,  1,  0, 
+                0,  0,  0, -1, -1,  0,  0,  0, -1,  0, 
+               -1,  0,  0,  0,  1,  0,  1,  0,  0,  1, 
+                0, -1,  0,  0,  0, -1,  0, -1,  0,  0, 
+               -1,  0,  0, -1,  1,  1, -1,  0,  0, -1, 
+                0,  0,  1, -1,  0, -1,  0,  0,  1,  0, 
+                0,  1,  1,  0,  0,  0, -1,  0,  0, -1, 
+                0,  0,  0, -1,  1,  0,  0, -1,  0,  1, 
+
+                0,  0,  1, -1, -1,  0,  0,  0,  0,  1, 
+                0,  0,  0,  0,  0,  0,  0,  0, -1,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  1, 
+                1,  0, -1,  0,  0,  0,  0,  0,  1,  0, 
+                0,  1,  0,  1,  0,  0,  0,  0,  0,  1, 
+                0,  0,  0,  0,  0,  1,  0,  0, -1,  0, 
+                1, -1,  0,  0,  0,  0,  1,  0, -1,  0, 
+               -1,  0, -1,  1,  0,  0,  1,  1,  0,  0, 
+                0,  0,  0,  0,  0,  1,  0,  0, -1,  0, 
+                0,  1, -1,  0,  1,  0,  0,  0,  1, -1, 
+
+               -1,  0,  0,  0, -1,  0,  1,  0,  0,  0, 
+                0,  0,  1, -1,  1,  0,  0,  0,  0,  0, 
+                0, -1,  0,  1,  0,  0,  0,  0,  0, -1, 
+                0,  0,  0,  0,  0,  0,  0,  1,  0,  1, 
+                0,  0,  0,  0, -1,  1,  0,  0,  0,  0, 
+                0,  0,  0,  0, -1,  0,  0,  0,  0,  1, 
+                0, -1,  0,  0,  0,  1,  0,  1, -1,  1, 
+                0,  0,  0,  0,  1,  1,  0,  1,  0,  0, 
+                0, -1,  1,  0,  1, -1,  0,  0, -1,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+
+                0,  0,  1,  0,  0,  0,  0,  0, -1,  0, 
+                0, -1,  0,  0,  0,  1,  0,  0,  0,  1, 
+                0,  0,  0,  0,  0,  0,  0, -1,  0,  0, 
+                0,  0, -1,  1,  0,  1,  0,  0, -1,  0, 
+                0, -1,  0,  0,  0,  0,  0,  0,  0,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+               -1,  0,  0,  1,  0, -1,  1,  0,  1,  0, 
+               -1,  1,  0,  0,  0,  0,  0,  0,  0, -1, 
+                0,  0,  0,  0,  0,  1,  1,  0, -1,  0, 
+               -1,  0,  0,  0, -1, -1,  0,  0,  0,  0, 
+
+                0,  0, -1, -1,  0,  1, -1,  0,  0,  0, 
+                0, -1
+       };
+
+int16_t g[] = {
+               -1,  0,  0,  0,  0,  0,  0,  0,  1,  0, 
+                0,  0,  0,  0,  0,  0,  0,  1,  1,  0, 
+                1,  0,  0,  0,  1,  0, -1,  0,  0,  0, 
+                0,  0,  0,  0, -1,  1,  0,  0,  0,  0, 
+                0,  0, -1,  1,  0,  0,  0,  0,  0, -1, 
+               -1,  0, -1,  0,  1, -1,  0,  0,  0,  0, 
+                0,  0,  0,  0,  1, -1,  0,  0,  1,  0, 
+                0,  1,  0, -1,  1,  0,  0,  0,  0,  0, 
+                0,  0,  0,  0, -1,  1,  0, -1,  0,  1, 
+                0,  1,  0,  0,  1,  0, -1,  1,  0,  0, 
+
+                0,  0,  0,  0,  0,  0,  1,  1,  0,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  1,  0, 
+                0,  0,  0,  0,  0,  0,  1,  0,  0,  1, 
+                1,  1,  1,  0,  0,  0,  0,  0, -1,  0, 
+                1,  1,  0,  0,  1,  0, -1,  1,  0,  0, 
+               -1,  0,  0, -1,  0,  0,  0,  0,  0,  0, 
+                0,  0,  0, -1,  0, -1,  1,  0,  0, -1, 
+                0,  0, -1,  0,  0,  0,  0,  0,  0, -1, 
+                0,  0,  1,  0,  0,  0,  0, -1,  0,  1, 
+                0,  0, -1,  0,  0,  0,  0,  0,  0,  0, 
+
+                0,  0,  0,  0,  0,  0,  0,  1,  0,  1, 
+                1,  0, -1,  0,  0,  0,  0,  1,  1,  0, 
+                0, -1,  1,  1,  0,  0,  0, -1,  1,  0, 
+                0,  1,  0,  0,  0,  0, -1,  1,  0,  0, 
+               -1,  0,  1,  0,  0,  0, -1, -1,  0,  0, 
+                0, -1,  1,  1,  0,  1,  0,  0,  0,  0, 
+                0,  0,  1,  0,  0,  1,  1,  1,  1,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+                0,  0,  0,  0,  1,  0,  0,  0,  0,  1, 
+                0, -1, -1,  0, -1,  0, -1, -1,  1, -1, 
+
+               -1,  0,  0,  0, -1,  0,  0,  0,  0,  0, 
+               -1,  1,  0, -1,  0,  0,  0,  0,  0,  0, 
+                1,  0,  0,  0,  0,  0, -1,  0,  0, -1, 
+                1,  0,  0, -1,  0,  0,  0,  0,  0, -1, 
+                0,  0, -1,  0,  0,  0, -1,  0,  0,  0, 
+                0,  0,  0,  1,  1,  0,  0,  0,  0,  0, 
+               -1,  0,  0,  0,  0,  0, -1,  0,  0,  0, 
+                1,  1,  1,  0,  0,  0,  0, -1, -1,  0, 
+               -1,  0,  1,  0,  0,  0,  0,  0, -1,  0, 
+               -1,  0,  0,  0,  0,  1,  0,  0, -1, -1, 
+
+                1,  0,  0,  0,  0,  0,  0,  0, -1,  0, 
+                0,  0,  1, -1,  0,  0,  1,  0,  0,  1, 
+               -1,  0,  1,  0,  1,  1,  0,  0,  0,  0, 
+                0,  0,  0,  0,  0,  0,  0,  0,  0,  0, 
+                0, -1,  0,  0,  0,  0,  1,  1,  0,  0, 
+                1,  0, -1,  0,  0,  1,  1,  0,  0,  0, 
+                0,  0, -1,  0,  0,  0,  0,  0,  1,  0, 
+                0,  0,  1, -1, -1,  0,  0,  0,  0,  1, 
+               -1, -1,  1,  0,  0,  0,  0,  0,  0,  1, 
+                0,  0,  0,  1,  0,  0,  0,  1,  1,  0, 
+
+               -1,  0,  1,  0,  0,  0,  0,  0,  0,  0, 
+                0, -1 
+};
 
        while (TRUE)
        {
@@ -180,6 +429,61 @@ bliss_private_key_t *bliss_private_key_gen(key_type_t type, va_list args)
        this = bliss_private_key_create_empty();
        this->key_size = key_size;
 
+       /* We derive the public key from the private key using the FFT */
+       fft = bliss_fft_create(&bliss_fft_12289_512);
+       n = fft->get_size(fft);
+       q = fft->get_modulus(fft);
+
+       /* Compute 2g + 1 */
+       for (i = 0; i < n; i++)
+       {
+               g[i] *= 2;
+       }
+       g[0] += 1;
+
+       l2_norm = wrapped_product(f, f, n, 0) + wrapped_product(g, g, n, 0);
+       nks = nks_norm(f, g, n);
+       DBG1(DBG_LIB, "L2 norm of s1||s2: %d, Nk(S) = %u", l2_norm, nks);
+
+       F = malloc(n * sizeof(uint32_t));
+       G = malloc(n * sizeof(uint32_t));
+       A = malloc(n * sizeof(uint32_t));
+       a = malloc(n * sizeof(uint32_t));
+
+       /* Convert signed arrays to unsigned arrays before FFT */
+       for (i = 0; i < n; i++)
+       {
+               F[i] = (f[i] < 0) ? f[i] + q : f[i];
+               G[i] = (g[i] < 0) ? g[i] + q : g[i];
+       }
+       fft->transform(fft, F, F, FALSE);
+       fft->transform(fft, G, G, FALSE);
+
+       for (i = 0; i < n; i++)
+       {
+               if (F[i] == 0)
+               {
+                       DBG1(DBG_LIB, "F[%d] is zero", i);
+               }
+               A[i] = invert(F[i], q);
+               A[i] = (G[i] * A[i]) % q;
+       }
+       fft->transform(fft, A, a, TRUE);
+
+       DBG1(DBG_LIB, "   i   f   g     a     F     G     A");
+       for (i = 0; i < n; i++)
+       {
+               DBG1(DBG_LIB, "%4d %3d %3d %5u %5u %5u %5u",
+                        i, f[i], g[i], a[i], F[i], G[i], A[i]);
+       }
+
+       /* Cleanup */
+       fft->destroy(fft);
+       free(a);
+       free(A);
+       free(F);
+       free(G);
+
        return &this->public;
 }
 
index c9b93a6..cb4ff80 100644 (file)
@@ -15,7 +15,7 @@
 
 /**
  * @defgroup bliss_private_key bliss_private_key
- * @{ @ingroup gmp_p
+ * @{ @ingroup bliss_p
  */
 
 #ifndef BLISS_PRIVATE_KEY_H_
index 71e0e4e..e211f4b 100644 (file)
@@ -15,7 +15,7 @@
 
 /**
  * @defgroup bliss_public_key bliss_public_key
- * @{ @ingroup gmp_p
+ * @{ @ingroup bliss_p
  */
 
 #ifndef BLISS_PUBLIC_KEY_H_