Annotation of embedaddon/strongswan/src/libstrongswan/plugins/bliss/bliss_huffman.c, revision 1.1.1.1

1.1       misho       1: /*
                      2:  * Copyright (C) 2014 Tobias Brunner
                      3:  * Copyright (C) 2014 Andreas Steffen
                      4:  * HSR Hochschule fuer Technik Rapperswil
                      5:  *
                      6:  * This program is free software; you can redistribute it and/or modify it
                      7:  * under the terms of the GNU General Public License as published by the
                      8:  * Free Software Foundation; either version 2 of the License, or (at your
                      9:  * option) any later version.  See <http://www.fsf.org/copyleft/gpl.txt>.
                     10:  *
                     11:  * This program is distributed in the hope that it will be useful, but
                     12:  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
                     13:  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
                     14:  * for more details.
                     15:  */
                     16: 
                     17: #include "bliss_param_set.h"
                     18: 
                     19: #include <library.h>
                     20: 
                     21: #include <stdio.h>
                     22: #include <math.h>
                     23: 
                     24: typedef struct tuple_t tuple_t;
                     25: 
                     26: struct tuple_t {
                     27:        int8_t z1;
                     28:        int8_t z2;
                     29:        uint16_t index;
                     30:        uint16_t bits;
                     31:        uint32_t code;
                     32: };
                     33: 
                     34: typedef struct node_t node_t;
                     35: 
                     36: struct node_t {
                     37:        node_t *next;
                     38:        node_t *l;
                     39:        node_t *r;
                     40:        tuple_t *tuple;
                     41:        double p;
                     42:        uint16_t depth;
                     43:        uint16_t index;
                     44: };
                     45: 
                     46: static void print_node(node_t *node)
                     47: {
                     48:        if (node->tuple)
                     49:        {
                     50:                fprintf(stderr, "(%1d,%2d)", node->tuple->z1, node->tuple->z2);
                     51:        }
                     52:        else
                     53:        {
                     54:                fprintf(stderr, "      ");
                     55:        }
                     56:        fprintf(stderr, "  %18.16f\n", node->p);
                     57: }
                     58: 
                     59: static double code_node(node_t *node, int *index, uint8_t bits, uint32_t code)
                     60: {
                     61:        double code_length = 0;
                     62: 
                     63:        node->index = (*index)++;
                     64: 
                     65:        if (node->tuple)
                     66:        {
                     67:                node->tuple->code = code;
                     68:                node->tuple->bits = bits;
                     69:                code_length += node->p * bits;
                     70:        }
                     71:        if (node->l)
                     72:        {
                     73:                code_length += code_node(node->l, index, bits + 1, (code << 1));
                     74:        }
                     75:        if (node->r)
                     76:        {
                     77:                code_length += code_node(node->r, index, bits + 1, (code << 1) + 1);
                     78:        }
                     79: 
                     80:        return code_length;
                     81: 
                     82: }
                     83: 
                     84: static void write_node(node_t *node)
                     85: {
                     86:        int16_t node_0, node_1, tuple;
                     87: 
                     88:        node_0 = node->l ? node->l->index : BLISS_HUFFMAN_CODE_NO_NODE;
                     89:        node_1 = node->r ? node->r->index : BLISS_HUFFMAN_CODE_NO_NODE;
                     90:        tuple = node->tuple ? node->tuple->index : BLISS_HUFFMAN_CODE_NO_TUPLE;
                     91: 
                     92:        printf("\t{ %3d, %3d, %3d },  /* %3d: ", node_0, node_1, tuple, node->index);
                     93: 
                     94:        if (node->tuple)
                     95:        {
                     96:                printf("(%d,%2d) %2u bit%s ", node->tuple->z1, node->tuple->z2,
                     97:                           node->tuple->bits, (node->tuple->bits == 1) ? " " : "s");
                     98:        }
                     99:        printf("*/\n");
                    100: 
                    101:        if (node->l)
                    102:        {
                    103:                write_node(node->l);
                    104:        }
                    105:        if (node->r)
                    106:        {
                    107:                write_node(node->r);
                    108:        }
                    109: }
                    110: 
                    111: static void write_header(void)
                    112: {
                    113:        printf("/*\n");
                    114:        printf(" * Copyright (C) 2014 Andreas Steffen\n");
                    115:        printf(" * HSR Hochschule fuer Technik Rapperswil\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: 

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>