src/libstrongswan/plugins/bliss/bliss_huffman.c

   1 /*
   2  * Copyright (C) 2014 Tobias Brunner
   3  * Copyright (C) 2014 Andreas Steffen
   4  *
   5  * Copyright (C) secunet Security Networks AG
   6  *
   7  * This program is free software; you can redistribute it and/or modify it
   8  * under the terms of the GNU General Public License as published by the
   9  * Free Software Foundation; either version 2 of the License, or (at your
  10  * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
  11  *
  12  * This program is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15  * for more details.
  16  */
  17
  18 #include "bliss_param_set.h"
  19
  20 #include <library.h>
  21
  22 #include <stdio.h>
  23 #include <math.h>
  24
  25 typedef struct tuple_t tuple_t;
  26
  27 struct tuple_t {
  28         int8_t z1;
  29         int8_t z2;
  30         uint16_t index;
  31         uint16_t bits;
  32         uint32_t code;
  33 };
  34
  35 typedef struct node_t node_t;
  36
  37 struct node_t {
  38         node_t *next;
  39         node_t *l;
  40         node_t *r;
  41         tuple_t *tuple;
  42         double p;
  43         uint16_t depth;
  44         uint16_t index;
  45 };
  46
  47 static void print_node(node_t *node)
  48 {
  49         if (node->tuple)
  50         {
  51                 fprintf(stderr, "(%1d,%2d)", node->tuple->z1, node->tuple->z2);
  52         }
  53         else
  54         {
  55                 fprintf(stderr, "      ");
  56         }
  57         fprintf(stderr, "  %18.16f\n", node->p);
  58 }
  59
  60 static double code_node(node_t *node, int *index, uint8_t bits, uint32_t code)
  61 {
  62         double code_length = 0;
  63
  64         node->index = (*index)++;
  65
  66         if (node->tuple)
  67         {
  68                 node->tuple->code = code;
  69                 node->tuple->bits = bits;
  70                 code_length += node->p * bits;
  71         }
  72         if (node->l)
  73         {
  74                 code_length += code_node(node->l, index, bits + 1, (code << 1));
  75         }
  76         if (node->r)
  77         {
  78                 code_length += code_node(node->r, index, bits + 1, (code << 1) + 1);
  79         }
  80
  81         return code_length;
  82
  83 }
  84
  85 static void write_node(node_t *node)
  86 {
  87         int16_t node_0, node_1, tuple;
  88
  89         node_0 = node->l ? node->l->index : BLISS_HUFFMAN_CODE_NO_NODE;
  90         node_1 = node->r ? node->r->index : BLISS_HUFFMAN_CODE_NO_NODE;
  91         tuple = node->tuple ? node->tuple->index : BLISS_HUFFMAN_CODE_NO_TUPLE;
  92
  93         printf("\t{ %3d, %3d, %3d },  /* %3d: ", node_0, node_1, tuple, node->index);
  94
  95         if (node->tuple)
  96         {
  97                 printf("(%d,%2d) %2u bit%s ", node->tuple->z1, node->tuple->z2,
  98                            node->tuple->bits, (node->tuple->bits == 1) ? " " : "s");
  99         }
 100         printf("*/\n");
 101
 102         if (node->l)
 103         {
 104                 write_node(node->l);
 105         }
 106         if (node->r)
 107         {
 108                 write_node(node->r);
 109         }
 110 }
 111
 112 static void write_header(void)
 113 {
 114         printf("/*\n");
 115         printf(" * Copyright (C) 2014 Andreas Steffen\n");
 116         printf(" *\n");
 117         printf(" * Optimum Huffman code for BLISS-X signatures\n");
 118         printf(" *\n");
 119         printf(" * This file has been automatically generated by the"
 120                    " bliss_huffman utility\n");
 121         printf(" * Do not edit manually!\n");
 122         printf(" */\n\n");
 123 };
 124
 125 static void write_code_tables(int bliss_type, int n_z1, int n_z2, node_t *nodes,
 126                                                           tuple_t **tuples)
 127 {
 128         int index, i, k;
 129         uint32_t bit;
 130         double code_length;
 131
 132         printf("#include \"bliss_huffman_code.h\"\n\n");
 133
 134         printf("static bliss_huffman_code_node_t nodes[] = {\n");
 135         index = 0;
 136         code_length = code_node(nodes, &index, 0, 0);
 137         write_node(nodes);
 138         printf("};\n\n");
 139
 140         printf("static bliss_huffman_code_tuple_t tuples[] = {\n");
 141         index = 0;
 142         for (i = 0; i < n_z1; i++)
 143         {
 144                 if (i > 0)
 145                 {
 146                         printf("\n");
 147                 }
 148                 for (k = 1 - n_z2; k < n_z2; k++)
 149                 {
 150                         printf("\t{ %5u, %2u },  /* %3d: (%1d,%2d) ",
 151                                                 tuples[index]->code, tuples[index]->bits, index, i, k);
 152                         bit = 1 << (tuples[index]->bits - 1);
 153                         while (bit)
 154                         {
 155                                 printf("%s", (tuples[index]->code & bit) ? "1" : "0");
 156                                 bit >>= 1;
 157                         }
 158                         printf(" */\n");
 159                         index++;
 160                 }
 161         }
 162         printf("};\n\n");
 163         printf("/* code_length = %6.4f bits/tuple (%d bits) */\n\n",
 164                            code_length, (int)(512 * code_length + 1));
 165
 166         printf("bliss_huffman_code_t bliss_huffman_code_%d = {\n", bliss_type);
 167         printf("\t.n_z1 = %d,\n", n_z1);
 168         printf("\t.n_z2 = %d,\n", n_z2);
 169         printf("\t.tuples = tuples,\n");
 170         printf("\t.nodes  = nodes\n");
 171         printf("};\n");
 172 }
 173
 174 static void destroy_node(node_t *node)
 175 {
 176         if (node->l)
 177         {
 178                 destroy_node(node->l);
 179         }
 180         if (node->r)
 181         {
 182                 destroy_node(node->r);
 183         }
 184         free(node->tuple);
 185         free(node);
 186 }
 187
 188 static void remove_node(node_t *list, node_t **last, node_t *node)
 189 {
 190         node_t *current, *prev;
 191
 192         for (current = list->next, prev = list; current;
 193                  prev = current, current = current->next)
 194         {
 195                 if (current == node)
 196                 {
 197                         prev->next = current->next;
 198                         if (*last == current)
 199                         {
 200                                 *last = prev->next ?: prev;
 201                         }
 202                         break;
 203                 }
 204         }
 205 }
 206
 207 /**
 208  * Generate a Huffman code for the optimum encoding of BLISS signatures
 209  */
 210 int main(int argc, char *argv[])
 211 {
 212         const bliss_param_set_t *set;
 213         int dx, bliss_type, depth = 1, groups, groups_left, pairs = 1;
 214         int i_max = 9, k_max = 8, index_max = (2*k_max - 1) * i_max;
 215         int i, i_top, k, k_top;
 216         uint16_t index;
 217         double p, p_z1[i_max], p_z2[k_max], x_z1[i_max], x_z2[k_max];
 218         double t, x, x0, p_sum, entropy = 0, erf_i, erf_k, erf_0 = 0;
 219         tuple_t *tuple, *tuples[index_max];
 220         node_t *node, *node_l, *node_r, *nodes = NULL;
 221         node_t *node_list, *node_last;
 222
 223         if (argc < 2)
 224         {
 225                 fprintf(stderr, "usage: bliss_huffman <bliss type> [<pairs>]\n");
 226                 exit(1);
 227         }
 228         if (argc > 2)
 229         {
 230                 pairs = atoi(argv[2]);
 231         }
 232         fprintf(stderr, "%d code pairs with constant length\n\n", pairs);
 233         groups_left = groups = pairs >> 1;
 234
 235         bliss_type = atoi(argv[1]);
 236         set = bliss_param_set_get_by_id(bliss_type);
 237         if (!set)
 238         {
 239                 fprintf(stderr, "bliss type %d unsupported\n", bliss_type);
 240                 exit(1);
 241         }
 242         write_header();
 243         printf("/*\n");
 244         printf(" * Design: sigma = %u\n", set->sigma);
 245         printf(" *\n");
 246
 247         t = 1/(sqrt(2) * set->sigma);
 248
 249         /* Probability distribution for z1 */
 250         i_top = (set->B_inf + 255) / 256;
 251         p_sum = 0;
 252         x = 0;
 253
 254         for (i = 0; i < i_top; i++)
 255         {
 256                 x = min(x + 256, set->B_inf);
 257                 erf_i = erf(t*x);
 258                 p_z1[i] = erf_i - erf_0;
 259                 p_sum += p_z1[i];
 260                 erf_0 = erf_i;
 261                 x_z1[i] = x;
 262         }
 263
 264         /* Normalize and print the probability distribution for z1 */
 265         printf(" *   i  p_z1[i]\n");
 266         x0 = 0;
 267
 268         for (i = 0; i < i_top; i++)
 269         {
 270                 p_z1[i] /= p_sum;
 271                 printf(" *  %2d  %18.16f      %4.0f .. %4.0f\n", i, p_z1[i], x0, x_z1[i]);
 272                 x0 = x_z1[i];
 273         }
 274         printf(" *\n");
 275
 276         /* Probability distribution for z2 */
 277         dx = 1 << set->d;
 278         k_top = 1 + set->B_inf / dx;
 279         x = (dx >> 1) - 0.5;
 280         p_sum = 0;
 281
 282         for (k = 0; k < k_top; k++)
 283         {
 284
 285                 erf_k = erf(t*x) / 2;
 286                 p_z2[k] = (k == 0) ? 2*erf_k : erf_k - erf_0;
 287                 p_sum +=  (k == 0) ? p_z2[k] : 2*p_z2[k];
 288                 erf_0 = erf_k;
 289                 x_z2[k] = x;
 290                 x += dx;
 291         }
 292
 293         /* Normalize the probability distribution for z2 */
 294         for (k = 0; k < k_top; k++)
 295         {
 296                 p_z2[k] /= p_sum;
 297         }
 298
 299         /* Print the probability distribution for z2 */
 300         printf(" *   k  p_z2[k]  dx = %d\n", dx);
 301
 302         for (k = 1 - k_top; k < k_top; k++)
 303         {
 304
 305                 printf(" *  %2d  %18.16f  ",k, p_z2[abs(k)]);
 306                 if (k < 0)
 307                 {
 308                         printf(" %7.1f ..%7.1f\n", -x_z2[-k], -x_z2[-k-1]);
 309                 }
 310                 else if (k == 0)
 311                 {
 312                         printf(" %7.1f ..%7.1f\n", -x_z2[k], x_z2[k]);
 313                 }
 314                 else
 315                 {
 316                         printf(" %7.1f ..%7.1f\n", x_z2[k-1], x_z2[k]);
 317                 }
 318         }
 319         printf(" *\n");
 320
 321         /* Compute probabilities of tuples (z1, z2) */
 322         INIT(node_list);
 323         node_last = node_list;
 324         printf(" *  (i, k)  p\n");
 325         p_sum =0;
 326         index = 0;
 327
 328         for (i = 0; i < i_top; i++)
 329         {
 330                 for (k = 1 - k_top; k < k_top; k++)
 331                 {
 332                         p = p_z1[i] * p_z2[abs(k)];
 333                         printf(" *  (%1d,%2d)  %18.16f\n", i, k, p);
 334                         p_sum += p;
 335                         entropy += -log(p) * p;
 336
 337                         INIT(tuple,
 338                                 .z1 = i,
 339                                 .z2 = k,
 340                                 .index = index,
 341                         );
 342                         tuples[index++] = tuple;
 343
 344                         INIT(node,
 345                                 .p = p,
 346                                 .tuple = tuple,
 347                         );
 348                         node_last->next = node;
 349                         node_last = node;
 350                 }
 351                 printf(" *\n");
 352         }
 353         entropy /= log(2);
 354         printf(" *  p_sum   %18.16f\n", p_sum);
 355         printf(" *\n");
 356         printf(" * entropy = %6.4f bits/tuple (%d bits)\n",
 357                            entropy, (int)(512 * entropy));
 358         printf(" */\n\n");
 359
 360         /* Build Huffman tree */
 361         while (node_list->next != node_last)
 362         {
 363                 node_r = node_l = NULL;
 364
 365                 for (node = node_list->next; node; node = node->next)
 366                 {
 367                         if (pairs > 0)
 368                         {
 369                                 if (!node->tuple)
 370                                 {
 371                                         continue;
 372                                 }
 373                         }
 374                         else if (groups_left > 0)
 375                         {
 376                                 if (node->tuple || node->depth != depth)
 377                                 {
 378                                         continue;
 379                                 }
 380                         }
 381                         if (node_r == NULL || node->p < node_r->p)
 382                         {
 383                                 node_l = node_r;
 384                                 node_r = node;
 385                         }
 386                         else if (node_l == NULL || node->p < node_l->p)
 387                         {
 388                                 node_l = node;
 389                         }
 390                 }
 391
 392                 INIT(node,
 393                         .l = node_l,
 394                         .r = node_r,
 395                         .p = node_l->p + node_r->p,
 396                         .depth = 1 + max(node_l->depth, node_r->depth),
 397                         .tuple = NULL,
 398                 );
 399                 print_node(node_r);
 400                 print_node(node_l);
 401                 fprintf(stderr, "        %18.16f", node->p);
 402
 403                 remove_node(node_list, &node_last, node_l);
 404                 remove_node(node_list, &node_last, node_r);
 405                 node_last->next = node;
 406                 node_last = node;
 407
 408                 if (pairs > 0)
 409                 {
 410                         pairs--;
 411                 }
 412                 else if (groups > 0)
 413                 {
 414                         if (--groups_left == 0)
 415                         {
 416                                 groups >>= 1;
 417                                 groups_left = groups;
 418                                 depth++;
 419                         }
 420                 }
 421                 fprintf(stderr, "\n\n");
 422         }
 423
 424
 425         nodes = node_list->next;
 426
 427         write_code_tables(bliss_type, i_top, k_top, nodes, tuples);
 428
 429         destroy_node(nodes);
 430         destroy_node(node_list);
 431         exit(0);
 432 }
 433