1 /**
2 * \file djb.h
3 *
4 * \brief Dan Bernstein's "Optimizing linear maps mod 2"
5 *
6 * This code is a port of sort1.cpp available at http://binary.cr.yp.to/linearmod2.html
7 *
8 * Given a matrix A djb_compile(A) will compute a djb_t data structure which realises A with
9 * (heuristically) (m * n)/(log m - loglog m) XORs.
10 *
11 * It makes use of a binary heap written by Martin Kunev which is available at
12 * https://gist.github.com/martinkunev/1365481
13 *
14 * \author Martin Albrecht <martinralbrecht@googlemail.com>
15 */
16
17 #ifndef M4RI_DJB_H
18 #define M4RI_DJB_H
19
20 #include <m4ri/mzd.h>
21
22 /**
23 * \brief Specify source type of addition
24 */
25
26 typedef enum {
27 source_target, //< add from target matrix
28 source_source //< add from source matrix
29 } srctyp_t;
30
31 /**
32 * \brief DJB's optimized linear maps mod 2
33 */
34
35 typedef struct {
36 rci_t nrows; /*!< Number of rows of map */
37 rci_t ncols; /*!< Number of columns of map */
38 rci_t *target; /*!< target row at index i */
39 rci_t *source; /*!< source row at index i */
40 srctyp_t *srctyp; /*!< source type at index i */
41 rci_t length; /*!< length of target, source and srctype */
42 wi_t allocated; /*!< how much did we allocate already */
43 } djb_t;
44
45 /**
46 * Standard allocation chunk
47 */
48
49 #define M4RI_DJB_BASE_SIZE 64
50
51 /**
52 * Allocate a new DJB linear map
53 *
54 * \param nrows Number of rows
55 * \param ncols Number of columns
56 */
57
djb_init(rci_t nrows,rci_t ncols)58 static inline djb_t *djb_init(rci_t nrows, rci_t ncols) {
59 /* we want to use realloc, so we call unaligned malloc */
60 djb_t *m = (djb_t*)malloc(sizeof(djb_t));
61 if (m == NULL)
62 m4ri_die("malloc failed.\n");
63
64 m->nrows = nrows;
65 m->ncols = ncols;
66 m->target = (rci_t*)malloc(sizeof(rci_t) * M4RI_DJB_BASE_SIZE);
67 m->source = (rci_t*)malloc(sizeof(rci_t) * M4RI_DJB_BASE_SIZE);
68 m->srctyp = (srctyp_t*)malloc(sizeof(srctyp_t) * M4RI_DJB_BASE_SIZE);
69 m->length = 0;
70 m->allocated = M4RI_DJB_BASE_SIZE;
71
72 if (m->target == NULL || m->source == NULL || m->srctyp == NULL)
73 m4ri_die("malloc failed.\n");
74 return m;
75 }
76
77 /**
78 * Free a DJB linear maps
79 *
80 * \param m Map
81 */
82
djb_free(djb_t * m)83 static inline void djb_free(djb_t *m) {
84 free(m->target);
85 free(m->source);
86 free(m->srctyp);
87 free(m);
88 }
89
90 /**
91 * Add a new operation out[target] ^= srctype[source] to queue.
92 *
93 * \param z DJB linear map.
94 * \param target Output index
95 * \param source Input index
96 * \param srctyp Type of input (source_source or source_target)
97 */
98
djb_push_back(djb_t * z,rci_t target,rci_t source,srctyp_t srctyp)99 static inline void djb_push_back(djb_t *z, rci_t target, rci_t source, srctyp_t srctyp) {
100 assert((target < z->nrows) &&
101 ((source < z->ncols) | (srctyp != source_source)) &&
102 ((source < z->nrows) | (srctyp != source_target)));
103 if (z->length >= z->allocated) {
104 z->allocated += M4RI_DJB_BASE_SIZE;
105 z->target = (rci_t*)realloc(z->target, z->allocated*sizeof(rci_t));
106 z->source = (rci_t*)realloc(z->source, z->allocated*sizeof(rci_t));
107 z->srctyp = (srctyp_t*)realloc(z->srctyp, z->allocated*sizeof(srctyp_t));
108 }
109 z->target[z->length] = target;
110 z->source[z->length] = source;
111 z->srctyp[z->length] = srctyp;
112 z->length++;
113 }
114
115 /**
116 * Compile a new DJB linear map from A.
117 *
118 * \param A
119 */
120
121 djb_t *djb_compile(mzd_t *A);
122
123 /**
124 * \brief W = m*V
125 *
126 * Apply the linear map m to V and write the result in W.
127 *
128 * \param z DJB linear map.
129 * \param W Output matrix
130 * \param V Input matrix
131 */
132
133 void djb_apply_mzd(djb_t *z, mzd_t *W, const mzd_t *V);
134
135
136 /**
137 * Print infomrmation on linear map mA
138 */
139
djb_info(const djb_t * z)140 static inline void djb_info(const djb_t *z) {
141 double save = (double)z->length / (double)(z->nrows * z->ncols);
142 printf("%d x %d linear map in %d xors (cost: %.5f)\n", z->nrows, z->ncols, z->length, save);
143 }
144
145
146 #endif //M4RI_DJB_H
147