Annotation of embedaddon/strongswan/src/libstrongswan/plugins/bliss/bliss_huffman.c, revision 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>