1 #ifndef BALANCING_H_
2 #define BALANCING_H_
3
4 #include <stdint.h>
5 #include "macros.h" // for MIN
6 #include "mod_ul.h"
7
8 // #include "mf.h"
9
10 /* This is used only for sanity checking after dispatch. */
11 /* The primes are 2^64-59 and 2^63-25 */
12 #define DUMMY_VECTOR_COORD_VALUE(j) \
13 (UINT64_C(0xffffffffffffffc5) / (1 + (j)) \
14 + UINT64_C(0x7fffffffffffffe7) * (1 + (j)))
15 #define DUMMY_VECTOR_COORD_VALUE2(j) \
16 (DUMMY_VECTOR_COORD_VALUE((j) ^ 0xbeef) ^ ((DUMMY_VECTOR_COORD_VALUE((j) ^ 0xcafe) << 32)))
17
18 /* here's some magma code for generating these arbitrary vectors, but
19 * it's not checked.
20 nrows:=65536;
21 vv:=func<x|VectorSpace(GF(2),64)!Intseq(x mod 2^64,2,64)>;
22 p:=2^64-59;
23 q:=2^63-25;
24 i2p:=func<a|Polynomial(GF(2),Intseq(a,2))>;
25 p2i:=func<p|Seqint(ChangeUniverse(Eltseq(p),Integers()),2)>;
26 ixor:=func<a,b|p2i(i2p(a)+i2p(b))>;
27 vs:=func<v|vv(Seqint(ChangeUniverse(Eltseq(v),Integers()),2)*2^32)>;
28 arbitrary1:=func<j|vv(p div j + q * j)>;
29 arbitrary2:=[arbitrary1(1+ixor(j, 0xbeef)) + vs(arbitrary1(1+ixor(j, 0xcafe))):j in [0..nrows-1]];
30 */
31
32 /* balancing structure */
33
34 #define FLAG_COLPERM 1
35 #define FLAG_ROWPERM 2
36 #define FLAG_REPLICATE 8 /* work with square matrices (padding if
37 need be), and replicate one
38 permutation to the other side, so that
39 we get conjugated permutations.
40 */
41
42 #define BALANCING_MAGIC UINT32_C(0xba1a0000)
43
44 struct balancing_header_s {
45 uint32_t zero; /* previous versions had neither zero nor magic.
46 By enforcing a zero field here, we make sure
47 that older and newer code will choke on
48 incompatible balancing files */
49 uint32_t magic;
50 uint32_t nh;
51 uint32_t nv;
52 uint32_t nrows;
53 uint32_t ncols;
54 uint32_t nzrows;
55 uint32_t nzcols;
56 uint64_t ncoeffs;
57 uint32_t checksum;
58 uint32_t flags;
59 // pshuf indicates two integers a,b such that the COLUMN i of the input
60 // matrix is in fact mapped to column a*i+b mod n in the matrix we work
61 // with. pshuf_inv indicates the inverse permutation. a and b do
62 // **NOT** depend on nh and nv. Therefore all balancing are
63 // compatible even though a permutation is being used.
64 // a == rhuf[0] is the smallest integer greater than floor(sqrt(n)) which
65 // is coprime to n. b is 42. pshuf_inv[0,1] are computed accordingly.
66 //
67 // This means that on all occasions, we are in fact working with the
68 // matrix M*S, for some permutation S. As long as we are concerned
69 // with left nullspace, this does not make any difference, but it
70 // does for right nullspace !
71 /* Note: pshuf and pshuf_inv have to do with what is called
72 * "shuffled-product" elsewhere (and has actually become the default)
73 */
74 uint32_t pshuf[2];
75 uint32_t pshuf_inv[2];
76 };
77 typedef struct balancing_header_s balancing_header[1];
78 typedef struct balancing_header_s * balancing_header_ptr;
79
80 struct balancing_s {
81 balancing_header h;
82 uint32_t trows; // number of rows including padding.
83 uint32_t tcols; // number of cols including padding.
84 uint32_t * rowperm; // row index for new mat. --> row index for old mat.
85 uint32_t * colperm; // might be equal to colperm.
86 };
87 typedef struct balancing_s balancing[1];
88 typedef struct balancing_s * balancing_ptr;
89
90 #ifdef __cplusplus
91 extern "C" {
92 #endif
93
94 /* Once the flags and perm[] fields have been provided, the caller must
95 * call _finalize() in order to 1) update the trows and tcols fields 2)
96 * compute the checksum of the balancing.
97 */
98 extern void balancing_set_row_col_count(balancing_ptr bal);
99 extern void balancing_finalize(balancing_ptr bal);
100 extern void balancing_write_inner(balancing_ptr bal, const char *);
101 extern void balancing_write(balancing_ptr bal, const char * , const char *);
102 extern void balancing_read(balancing_ptr bal, const char *);
103 extern void balancing_read_header(balancing_ptr bal, const char * filename);
104 extern void balancing_clear(balancing_ptr bal);
105 extern void balancing_init(balancing_ptr bal);
106
107 /* helper for the functions below */
balancing_index_shuffle_common_(unsigned long r,unsigned long n,uint32_t * shuf)108 static inline unsigned long balancing_index_shuffle_common_(unsigned long r, unsigned long n, uint32_t * shuf)
109 {
110 modulusul_t M;
111 modul_initmod_ul(M, n);
112 residueul_t x,a,b;
113 modul_init_noset0(a, M);
114 modul_init_noset0(b, M);
115 modul_init_noset0(x, M);
116 modul_set_ul(a, shuf[0], M);
117 modul_set_ul(b, shuf[1], M);
118 modul_set_ul(x, r, M);
119 modul_mul(x, x, a, M);
120 modul_add(x, x, b, M);
121 unsigned long r1 = modul_get_ul(x, M);
122 modul_clear(a, M);
123 modul_clear(b, M);
124 modul_clear(x, M);
125 modul_clearmod(M);
126 return r1;
127 }
128
129 /* These two relate to the global permutation represented by the rshuf /
130 * rshuf_inv arrays */
balancing_pre_shuffle(balancing_ptr bal,unsigned long r)131 static inline unsigned long balancing_pre_shuffle(balancing_ptr bal, unsigned long r)
132 {
133 unsigned int K = MIN(bal->h->ncols, bal->h->nrows);
134 if (r >= K) return r;
135 return balancing_index_shuffle_common_(r, K, bal->h->pshuf);
136 }
balancing_pre_unshuffle(balancing_ptr bal,unsigned long r)137 static inline unsigned long balancing_pre_unshuffle(balancing_ptr bal, unsigned long r)
138 {
139 unsigned int K = MIN(bal->h->ncols, bal->h->nrows);
140 if (r >= K) return r;
141 return balancing_index_shuffle_common_(r, K, bal->h->pshuf_inv);
142 }
143
144
145 #ifdef __cplusplus
146 }
147 #endif
148
149 #endif /* BALANCING_H_ */
150