code/cauchy_fermat_erasure_coding.hh at master · aicodix/code · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
/*
Cauchy Fermat Prime Field Erasure Coding

Copyright 2024 Ahmet Inan <inan@aicodix.de>
*/

#pragma once

namespace CODE {

template <typename PF, typename IO, int MAX_LEN>
struct CauchyFermatErasureCoding
{
	static_assert(MAX_LEN < int(PF::P-1), "Block length must be smaller than largest field value");
	PF temp[MAX_LEN];
	typedef unsigned used_word;
	static constexpr int used_width = 8 * sizeof(used_word);
	static constexpr int used_length = (MAX_LEN + used_width - 1) / used_width;
	used_word used_values[used_length];
	PF row_num, row_den;
	CauchyFermatErasureCoding(){}
	// $a_{ij} = \frac{1}{x_i + y_j}$
	PF cauchy_matrix(int i, int j)
	{
		PF row(i), col(j);
		return rcp(row + col);
	}
	// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
	PF inverse_cauchy_matrix(const IO *rows, int i, int j, int n)
	{
#if 0
		PF row_j(rows[j]), col_i(i);
		PF prod_xy(1), prod_x(1), prod_y(1);
		for (int k = 0; k < n; k++) {
			PF row_k(rows[k]), col_k(k);
			prod_xy *= (row_j + col_k) * (row_k + col_i);
			if (k != j)
				prod_x *= (row_j - row_k);
			if (k != i)
				prod_y *= (col_i - col_k);
		}
		return prod_xy / ((row_j + col_i) * prod_x * prod_y);
#else
		PF row_j(rows[j]), col_i(i);
		if (j == 0) {
			PF num(1), den(1);
			for (int k = 0; k < n; k++) {
				PF row_k(rows[k]), col_k(k);
				num *= row_k + col_i;
				if (k != i)
					den *= col_i - col_k;
			}
			row_num = num;
			row_den = den;
		}
		PF num(row_num), den(row_den);
		for (int k = 0; k < n; k++) {
			PF row_k(rows[k]), col_k(k);
			num *= row_j + col_k;
			if (k != j)
				den *= row_j - row_k;
		}
		return num / ((row_j + col_i) * den);
#endif
	}
	void mac(const IO *a, PF b, int len, bool first)
	{
		if (first) {
			for (int i = 0; i < len; i++)
				temp[i] = b * PF(a[i]);
		} else {
			for (int i = 0; i < len; i++)
				temp[i] += b * PF(a[i]);
		}
	}
	void mac_sub(IO *c, const IO *a, PF b, IO s, int len, bool first, bool last)
	{
		int v = PF::P-1;
		if (first && last) {
			for (int i = 0; i < len; i++)
				c[i] = (b * PF(a[i] == s ? v : a[i]))();
		} else if (first) {
			for (int i = 0; i < len; i++)
				temp[i] = b * PF(a[i] == s ? v : a[i]);
		} else if (last) {
			for (int i = 0; i < len; i++)
				c[i] = ((temp[i] + b * PF(a[i] == s ? v : a[i])))();
		} else {
			for (int i = 0; i < len; i++)
				temp[i] += b * PF(a[i] == s ? v : a[i]);
		}
	}
	int find_unused(int block_len)
	{
		int limit = (block_len + used_width - 1) / used_width;
		for (int i = 0; i < limit; ++i)
			used_values[i] = 0;
		for (int i = 0; i < block_len; ++i)
			if (int(temp[i]())/used_width < limit)
				used_values[int(temp[i]())/used_width] |= 1 << int(temp[i]())%used_width;
		int s = 0;
		while (s/used_width < limit && used_values[s/used_width] & 1 << s%used_width)
			++s;
		return s;
	}
	int encode(const IO *data, IO *block, int block_id, int block_len, int block_cnt)
	{
		assert(block_id >= block_cnt && block_id < int(PF::P) / 2);
		assert(block_len <= MAX_LEN);
		for (int k = 0; k < block_cnt; k++) {
			PF a_ik = cauchy_matrix(block_id, k);
			mac(data + block_len * k, a_ik, block_len, !k);
		}
		int sub = find_unused(block_len);
		for (int i = 0; i < block_len; ++i)
			block[i] = temp[i]() == PF::P-1 ? sub : temp[i]();
		return sub;
	}
	void decode(IO *data, const IO *blocks, const IO *block_subs, const IO *block_ids, int block_idx, int block_len, int block_cnt)
	{
		assert(block_len <= MAX_LEN);
		for (int k = 0; k < block_cnt; k++) {
			PF b_ik = inverse_cauchy_matrix(block_ids, block_idx, k, block_cnt);
			mac_sub(data, blocks + block_len * k, b_ik, block_subs[k], block_len, !k, k == block_cnt - 1);
		}
	}
};

}