Automatic generation of optimized Huffman codes
authorAndreas Steffen <andreas.steffen@strongswan.org>
Tue, 9 Dec 2014 06:48:51 +0000 (07:48 +0100)
committerAndreas Steffen <andreas.steffen@strongswan.org>
Tue, 9 Dec 2014 10:58:18 +0000 (11:58 +0100)
src/libstrongswan/plugins/bliss/.gitignore
src/libstrongswan/plugins/bliss/Makefile.am
src/libstrongswan/plugins/bliss/bliss_huffman.c [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_huffman_code.c [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_huffman_code.h [new file with mode: 0644]
src/libstrongswan/plugins/bliss/bliss_param_set.c
src/libstrongswan/plugins/bliss/bliss_param_set.h

index 94f77f0..fcce5a9 100644 (file)
@@ -1 +1,5 @@
 bliss_tests
+bliss_huffman
+bliss_huffman_code_1.c
+bliss_huffman_code_3.c
+bliss_huffman_code_4.c
index 3db4819..c4db154 100644 (file)
@@ -21,10 +21,31 @@ libstrongswan_bliss_la_SOURCES = \
        bliss_bitpacker.h bliss_bitpacker.c \
        bliss_fft.h bliss_fft.c \
        bliss_fft_params.h bliss_fft_params.c \
+       bliss_huffman_code.h bliss_huffman_code.c \
+       bliss_huffman_code_1.c bliss_huffman_code_3.c bliss_huffman_code_4.c \
        bliss_sampler.h bliss_sampler.c
 
 libstrongswan_bliss_la_LDFLAGS = -module -avoid-version
 
+noinst_PROGRAMS = bliss_huffman
+
+bliss_huffman_SOURCES = bliss_huffman.c
+bliss_huffman_LDADD = \
+       $(top_builddir)/src/libstrongswan/libstrongswan.la -lm \
+       bliss_param_set.o bliss_fft_params.o
+
+bliss_huffman_code_1.c :       bliss_huffman bliss_huffman_code.h
+       $(AM_V_GEN) \
+       ./bliss_huffman 1 8 > $@
+
+bliss_huffman_code_3.c :       bliss_huffman bliss_huffman_code.h
+       $(AM_V_GEN) \
+       ./bliss_huffman 3 16 > $@
+
+bliss_huffman_code_4.c :       bliss_huffman bliss_huffman_code.h
+       $(AM_V_GEN) \
+       ./bliss_huffman 4 32 > $@
+
 TESTS = bliss_tests
 
 check_PROGRAMS = $(TESTS)
diff --git a/src/libstrongswan/plugins/bliss/bliss_huffman.c b/src/libstrongswan/plugins/bliss/bliss_huffman.c
new file mode 100644 (file)
index 0000000..890dfba
--- /dev/null
@@ -0,0 +1,408 @@
+/*
+ * 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_param_set.h"
+
+#include <library.h>
+
+#include <stdio.h>
+#include <math.h>
+
+typedef struct tuple_t tuple_t;
+
+struct tuple_t {
+       int8_t z1;
+       int8_t z2;
+       uint16_t index;
+       uint16_t bits;
+       uint32_t code;
+};
+
+typedef struct node_t node_t;
+
+struct node_t {
+       node_t *l;
+       node_t *r;
+       tuple_t *tuple;
+       double p;
+       uint16_t depth;
+       uint16_t index;
+};
+
+static void print_node(node_t *node)
+{
+       if (node->tuple)
+       {
+               fprintf(stderr, "(%1d,%2d)", node->tuple->z1, node->tuple->z2);
+       }
+       else
+       {
+               fprintf(stderr, "      ");
+       }
+       fprintf(stderr, "  %18.16f\n", node->p);
+}
+
+static double code_node(node_t *node, int *index, uint8_t bits, uint32_t code)
+{
+       double code_length = 0;
+
+       node->index = (*index)++;
+
+       if (node->tuple)
+       {
+               node->tuple->code = code;
+               node->tuple->bits = bits;
+               code_length += node->p * bits;
+       }
+       if (node->l)
+       {
+               code_length += code_node(node->l, index, bits + 1, (code << 1));
+       }
+       if (node->r)
+       {
+               code_length += code_node(node->r, index, bits + 1, (code << 1) + 1);
+       }
+
+       return code_length;
+
+}
+
+#define NO_NODE                -1
+#define NO_TUPLE       -1
+
+static void write_node(node_t *node)
+{
+       int16_t node_0, node_1, tuple;
+
+       node_0 = node->l ? node->l->index : NO_NODE;
+       node_1 = node->r ? node->r->index : NO_NODE;
+       tuple = node->tuple ? node->tuple->index : NO_TUPLE;
+       printf("\t{ %3d, %3d, %3d },  /* %3d: ", node_0, node_1, tuple, node->index);
+
+       if (node->tuple)
+       {
+               printf("(%d,%2d) %2u bit%s ", node->tuple->z1, node->tuple->z2,
+                          node->tuple->bits, (node->tuple->bits == 1) ? " " : "s");
+       }
+       printf("*/\n");
+
+       if (node->l)
+       {
+               write_node(node->l);
+       }
+       if (node->r)
+       {
+               write_node(node->r);
+       }
+}
+
+static void write_code_tables(int bliss_type, int n_z1, int n_z2, node_t *nodes,
+                                                         tuple_t **tuples)
+{
+       int index, i, k;
+       uint32_t bit;
+       double code_length;
+
+       printf("/*\n");
+       printf(" * Copyright (C) 2014 Andreas Steffen\n");
+       printf(" * HSR Hochschule fuer Technik Rapperswil\n");
+       printf(" *\n");
+       printf(" * Optimum Huffman code for %N signatures\n",
+                       bliss_param_set_id_names, bliss_type);
+       printf(" *\n");
+       printf(" * This file has been automatically generated by the"
+                  " bliss_huffman utility\n");
+       printf(" * Do not edit manually!\n");
+       printf(" */\n\n");
+       printf("#include \"bliss_huffman_code.h\"\n\n");
+
+       printf("static bliss_huffman_code_node_t nodes[] = {\n");
+       index = 0;
+       code_length = code_node(nodes, &index, 0, 0);
+       write_node(nodes);
+       printf("};\n\n");
+
+       printf("static bliss_huffman_code_tuple_t tuples[] = {\n");
+       index = 0;
+       for (i = 0; i < n_z1; i++)
+       {
+               if (i > 0)
+               {
+                       printf("\n");
+               }
+               for (k = 1 - n_z2; k < n_z2; k++)
+               {
+                       printf("\t{ %5u, %2u },  /* %3d: (%1d,%2d) ",
+                                               tuples[index]->code, tuples[index]->bits, index, i, k);
+                       bit = 1 << (tuples[index]->bits - 1);
+                       while (bit)
+                       {
+                               printf("%s", (tuples[index]->code & bit) ? "1" : "0");
+                               bit >>= 1;
+                       }
+                       printf(" */\n");
+                       index++;
+               }
+       }
+       printf("};\n\n");
+       printf("/* code_length = %6.4f bits/tuple (%d bits) */\n\n",
+                          code_length, (int)(512 * code_length + 1));
+
+       printf("bliss_huffman_code_t bliss_huffman_code_%d = {\n", bliss_type);
+       printf("\t.n_z1 = %d,\n", n_z1);
+       printf("\t.n_z2 = %d,\n", n_z2);
+       printf("\t.tuples = tuples,\n");
+       printf("\t.nodes  = nodes\n");
+       printf("};\n");
+}
+
+static void destroy_node(node_t *node)
+{
+       if (node->l)
+       {
+               destroy_node(node->l);
+       }
+       if (node->r)
+       {
+               destroy_node(node->r);
+       }
+       free(node->tuple);
+       free(node);
+}
+
+/**
+ * Generate a Huffman code for the optimum encoding of BLISS signatures
+ */
+int main(int argc, char *argv[])
+{
+       bliss_param_set_t *set;
+       int dx, bliss_type, depth = 1, groups, groups_left, pairs = 1;
+       int i_max = 9, k_max = 8, index_max = (2*k_max - 1) * i_max;
+       int i, i_top, k, k_top; 
+       uint16_t index;
+       double p, p_z1[i_max], p_z2[k_max], x_z1[i_max], x_z2[k_max];
+       double t, x, x0, p_sum, entropy = 0, erf_i, erf_k, erf_0 = 0;
+       tuple_t *tuple, *tuples[index_max];
+       node_t *node, *node_l, *node_r, *nodes = NULL;
+       linked_list_t *node_list;
+       enumerator_t *enumerator;
+
+       if (argc < 2)
+       {
+               fprintf(stderr, "usage: bliss_huffman <bliss type> [<pairs>]\n");
+               exit(1);
+       }
+       if (argc > 2)
+       {
+               pairs = atoi(argv[2]);
+       }
+       fprintf(stderr, "%d code pairs with constant length\n\n", pairs);
+       groups_left = groups = pairs >> 1;
+
+       library_init(NULL, "bliss_huffman");
+       lib->plugins->load(lib->plugins, "bliss");
+       atexit(library_deinit);
+
+       bliss_type = atoi(argv[1]);
+       set = bliss_param_set_get_by_id(bliss_type);
+       if (!set)
+       {
+               fprintf(stderr, "bliss type %d unsupported\n");
+               exit(1);
+       }
+       
+       t = 1/(sqrt(2) * set->sigma);
+
+       /* Probability distribution for z1 */
+       i_top = (set->B_inf + 255) / 256;
+       p_sum = 0;
+       x = 0;
+
+       for (i = 0; i < i_top; i++)
+       {
+               x = min(x + 256, set->B_inf);
+               erf_i = erf(t*x);
+               p_z1[i] = erf_i - erf_0;
+               p_sum += p_z1[i];
+               erf_0 = erf_i;
+               x_z1[i] = x;
+       }
+
+       /* Normalize and print the probability distribution for z1 */
+       fprintf(stderr, "i   p_z1[i]\n");
+       x0 = 0;
+
+       for (i = 0; i < i_top; i++)
+       {
+               p_z1[i] /= p_sum;
+               fprintf(stderr, "%1d   %18.16f     %4.0f .. %4.0f\n", i, p_z1[i],
+                                               x0, x_z1[i]);
+               x0 = x_z1[i];
+       }
+       fprintf(stderr, "\n");
+
+       /* Probability distribution for z2 */
+       dx = 1 << set->d;
+       k_top = 1 + set->B_inf / dx;
+       x = (double)(dx/2) - 0.5;
+       p_sum = 0;
+
+       for (k = 0; k < k_top; k++)
+       {
+               
+               erf_k = erf(t*x) / 2;
+               p_z2[k] = (k == 0) ? 2*erf_k : erf_k - erf_0;
+               p_sum +=  (k == 0) ? p_z2[k] : 2*p_z2[k];
+               erf_0 = erf_k;
+               x_z2[k] = x;
+               x += dx;
+       }
+
+       /* Normalize the probability distribution for z2 */
+       for (k = 0; k < k_top; k++)
+       {
+               p_z2[k] /= p_sum;
+       }
+       
+       /* Print the probability distribution for z2 */
+       fprintf(stderr, " k  p_z2[k]  dx = %d\n", dx);
+
+       for (k = 1 - k_top; k < k_top; k++)
+       {
+               
+               fprintf(stderr, "%2d  %18.16f  ",k, p_z2[abs(k)]);
+               if (k < 0)
+               {
+                       fprintf(stderr, "%7.1f ..%7.1f\n", -x_z2[-k], -x_z2[-k-1]);
+               } 
+               else if (k == 0)
+               {
+                       fprintf(stderr, "%7.1f ..%7.1f\n", -x_z2[k], x_z2[k]);
+               }
+               else
+               {
+                       fprintf(stderr, "%7.1f ..%7.1f\n", x_z2[k-1], x_z2[k]);
+               } 
+       }
+       fprintf(stderr, "\n");
+
+       /* Compute probabilities of tuples (z1, z2) */
+       node_list =linked_list_create();
+       fprintf(stderr, "(i, k)  p\n");
+       p_sum =0;
+       index = 0;
+
+       for (i = 0; i < i_top; i++)
+       {
+               for (k = 1 - k_top; k < k_top; k++)
+               {
+                       p = p_z1[i] * p_z2[abs(k)];
+                       fprintf(stderr, "(%1d,%2d)  %18.16f\n", i, k, p);
+                       p_sum += p;
+                       entropy += -log(p) * p;
+
+                       tuple = malloc_thing(tuple_t);
+                       node = malloc_thing(node_t);
+                       tuple->z1 = i;
+                       tuple->z2 = k;
+                       tuple->index = index;
+                       tuples[index++] = tuple;
+                       node->p = p;
+                       node->tuple = tuple;
+                       node->depth = 0;
+                       node->r = NULL;
+                       node->l = NULL;
+                       node_list->insert_last(node_list, node);
+               }
+       }
+       entropy /= log(2);
+       fprintf(stderr, "        %18.16f, entropy = %6.4f bits/tuple (%d bits)\n\n",
+                                                        p_sum, entropy, (int)(512 * entropy));
+
+       /* Build Huffman tree */
+       while (node_list->get_count(node_list) > 1)
+       {
+               node_r = node_l = NULL;
+               enumerator = node_list->create_enumerator(node_list);
+
+               while (enumerator->enumerate(enumerator, &node))
+               {
+                       if (pairs > 0)
+                       {
+                               if (!node->tuple)
+                               {
+                                       continue;
+                               }
+                       }
+                       else if (groups_left > 0)
+                       {
+                               if (node->tuple || node->depth != depth)
+                               {
+                                       continue;
+                               }
+                       }
+                       if (node_r == NULL || node->p < node_r->p)
+                       {
+                               node_l = node_r;
+                               node_r = node;
+                       }
+                       else if (node_l == NULL || node->p < node_l->p)
+                       {
+                               node_l = node;
+                       }
+               }
+               enumerator->destroy(enumerator);
+
+               node = malloc_thing(node_t);
+               node->l = node_l;
+               node->r = node_r;
+               node->p = node_l->p + node_r->p;
+               node->depth = 1 + max(node_l->depth, node_r->depth);
+               node->tuple = NULL;
+
+               print_node(node_r);
+               print_node(node_l);
+               fprintf(stderr, "        %18.16f", node->p);
+
+               node_list->remove(node_list, node_l, NULL);
+               node_list->remove(node_list, node_r, NULL);
+               node_list->insert_last(node_list, node);
+
+
+               if (pairs > 0)
+               {
+                       pairs--;
+               }
+               else if (groups > 0)
+               {
+                       if (--groups_left == 0)
+                       {
+                               groups >>= 1;
+                               groups_left = groups;
+                               depth++;
+                       }
+               }
+               fprintf(stderr, " %3d\n\n", node_list->get_count(node_list));
+       }
+
+       node_list->remove_first(node_list, (void**)&nodes);
+       node_list->destroy(node_list);
+
+       write_code_tables(bliss_type, i_top, k_top, nodes, tuples);
+
+       destroy_node(nodes);
+       exit(0);
+}
+
diff --git a/src/libstrongswan/plugins/bliss/bliss_huffman_code.c b/src/libstrongswan/plugins/bliss/bliss_huffman_code.c
new file mode 100644 (file)
index 0000000..68a799b
--- /dev/null
@@ -0,0 +1,39 @@
+/*
+ * 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_huffman_code.h"
+
+extern bliss_huffman_code_t bliss_huffman_code_1;
+extern bliss_huffman_code_t bliss_huffman_code_3;
+extern bliss_huffman_code_t bliss_huffman_code_4;
+
+/**
+ * See header.
+ */
+bliss_huffman_code_t* bliss_huffman_code_get_by_id(bliss_param_set_id_t id)
+{
+       switch (id)
+       {
+               case BLISS_I:
+                       return &bliss_huffman_code_1;
+               case BLISS_III:
+                       return &bliss_huffman_code_3;
+               case BLISS_IV:
+                       return &bliss_huffman_code_4;
+               default:
+                       return NULL;
+       }
+}
+
diff --git a/src/libstrongswan/plugins/bliss/bliss_huffman_code.h b/src/libstrongswan/plugins/bliss/bliss_huffman_code.h
new file mode 100644 (file)
index 0000000..ede769e
--- /dev/null
@@ -0,0 +1,77 @@
+/*
+ * 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_huffman_code bliss_huffman_code
+ * @{ @ingroup bliss_p
+ */
+
+#ifndef BLISS_HUFFMAN_CODE_H_
+#define BLISS_HUFFMAN_CODE_H_
+
+#include "bliss_param_set.h"
+
+#include <library.h>
+
+typedef struct bliss_huffman_code_t bliss_huffman_code_t;
+typedef struct bliss_huffman_code_tuple_t bliss_huffman_code_tuple_t;
+typedef struct bliss_huffman_code_node_t bliss_huffman_code_node_t;
+
+struct bliss_huffman_code_tuple_t {
+       uint32_t code;
+       uint16_t bits;
+};
+
+struct bliss_huffman_code_node_t {
+       int16_t node_0;
+       int16_t node_1;
+       int16_t tuple;
+};
+
+/**
+ * Defines the Huffman code for the optimum encoding of a BLISS signature
+ */
+struct bliss_huffman_code_t {
+
+       /**
+        * Range of z1:  0..n_z1-1
+        */
+       uint16_t n_z1;
+
+       /**
+        * Range of z2:  -n_z2..n_z2
+        */
+       uint16_t n_z2;
+
+       /**
+        * Table of tuple codewords
+        */
+       bliss_huffman_code_tuple_t *tuples;
+
+       /**
+        * Table of binary decision nodes
+        */
+       bliss_huffman_code_node_t *nodes;
+};
+
+/**
+ * Get Optimum Huffman code for BLISS signature given by BLISS parameter set ID
+ *
+ * @param id   BLISS parameter set ID
+ * @return             Optimum Huffman code for BLISS signature
+*/
+bliss_huffman_code_t* bliss_huffman_code_get_by_id(bliss_param_set_id_t id);
+
+#endif /** BLISS_HUFFMAN_CODE_H_ @}*/
index 90a5bcc..cbdb56d 100644 (file)
@@ -132,6 +132,7 @@ static bliss_param_set_t bliss_param_sets[] = {
                .non_zero2 = 0,
                .kappa = 23,
                .nks_max = 46479,
+               .sigma = 215,
                .k_sigma = 254,
                .k_sigma_bits = 8,
                .c = c_bliss_i,
@@ -160,6 +161,7 @@ static bliss_param_set_t bliss_param_sets[] = {
                .non_zero2 = 16,
                .kappa = 30,
                .nks_max = 128626,
+               .sigma = 250,
                .k_sigma = 295,
                .k_sigma_bits = 9,
                .c = c_bliss_iii,
@@ -188,6 +190,7 @@ static bliss_param_set_t bliss_param_sets[] = {
                .non_zero2 = 31,
                .kappa = 39,
                .nks_max = 244669,
+               .sigma = 271,
                .k_sigma = 320,
                .k_sigma_bits = 9,
                .c = c_bliss_iv,
index 168b6f3..8c73afa 100644 (file)
@@ -25,6 +25,7 @@ typedef enum bliss_param_set_id_t bliss_param_set_id_t;
 typedef struct bliss_param_set_t bliss_param_set_t;
 
 #include "bliss_fft_params.h"
+#include "bliss_huffman_code.h"
 
 #include <library.h>
 
@@ -111,6 +112,11 @@ struct bliss_param_set_t {
        uint32_t nks_max;
 
        /**
+        * Standard deviation sigma
+        */
+       uint16_t sigma;
+
+       /**
         *  k_sigma = ceiling[ sqrt(2*ln 2) * sigma ]
         */
        uint16_t k_sigma;
@@ -164,6 +170,7 @@ struct bliss_param_set_t {
         * B_verify bound
         */
        uint32_t B_l2;
+
 };
 
 /**