shamir.c 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338
  1. /*
  2. * Implementation of the hazardous parts of the SSS library
  3. *
  4. * Copyright (c) 2017 Daan Sprenkels <hello@dsprenkels.com>
  5. * Copyright (c) 2019 SatoshiLabs
  6. *
  7. * Permission is hereby granted, free of charge, to any person obtaining
  8. * a copy of this software and associated documentation files (the "Software"),
  9. * to deal in the Software without restriction, including without limitation
  10. * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  11. * and/or sell copies of the Software, and to permit persons to whom the
  12. * Software is furnished to do so, subject to the following conditions:
  13. *
  14. * The above copyright notice and this permission notice shall be included
  15. * in all copies or substantial portions of the Software.
  16. *
  17. * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  18. * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19. * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  20. * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES
  21. * OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  22. * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23. * OTHER DEALINGS IN THE SOFTWARE.
  24. *
  25. * This code contains the actual Shamir secret sharing functionality. The
  26. * implementation of this code is based on the idea that the user likes to
  27. * generate/combine 32 shares (in GF(2^8)) at the same time, because a 256 bit
  28. * key will be exactly 32 bytes. Therefore we bitslice all the input and
  29. * unbitslice the output right before returning.
  30. *
  31. * This bitslice approach optimizes natively on all architectures that are 32
  32. * bit or more. Care is taken to use not too many registers, to ensure that no
  33. * values have to be leaked to the stack.
  34. *
  35. * All functions in this module are implemented constant time and constant
  36. * lookup operations, as all proper crypto code must be.
  37. */
  38. #include "shamir.h"
  39. #include <string.h>
  40. #include "memzero.h"
  41. static void bitslice(uint32_t r[8], const uint8_t* x, size_t len) {
  42. size_t bit_idx = 0, arr_idx = 0;
  43. uint32_t cur = 0;
  44. memset(r, 0, sizeof(uint32_t[8]));
  45. for(arr_idx = 0; arr_idx < len; arr_idx++) {
  46. cur = (uint32_t)x[arr_idx];
  47. for(bit_idx = 0; bit_idx < 8; bit_idx++) {
  48. r[bit_idx] |= ((cur >> bit_idx) & 1) << arr_idx;
  49. }
  50. }
  51. }
  52. static void unbitslice(uint8_t* r, const uint32_t x[8], size_t len) {
  53. size_t bit_idx = 0, arr_idx = 0;
  54. uint32_t cur = 0;
  55. memset(r, 0, sizeof(uint8_t) * len);
  56. for(bit_idx = 0; bit_idx < 8; bit_idx++) {
  57. cur = (uint32_t)x[bit_idx];
  58. for(arr_idx = 0; arr_idx < len; arr_idx++) {
  59. r[arr_idx] |= ((cur >> arr_idx) & 1) << bit_idx;
  60. }
  61. }
  62. }
  63. static void bitslice_setall(uint32_t r[8], const uint8_t x) {
  64. size_t idx = 0;
  65. for(idx = 0; idx < 8; idx++) {
  66. r[idx] = -((x >> idx) & 1);
  67. }
  68. }
  69. /*
  70. * Add (XOR) `r` with `x` and store the result in `r`.
  71. */
  72. static void gf256_add(uint32_t r[8], const uint32_t x[8]) {
  73. size_t idx = 0;
  74. for(idx = 0; idx < 8; idx++) r[idx] ^= x[idx];
  75. }
  76. /*
  77. * Safely multiply two bitsliced polynomials in GF(2^8) reduced by
  78. * x^8 + x^4 + x^3 + x + 1. `r` and `a` may overlap, but overlapping of `r`
  79. * and `b` will produce an incorrect result! If you need to square a polynomial
  80. * use `gf256_square` instead.
  81. */
  82. static void gf256_mul(uint32_t r[8], const uint32_t a[8], const uint32_t b[8]) {
  83. /* This function implements Russian Peasant multiplication on two
  84. * bitsliced polynomials.
  85. *
  86. * I personally think that these kinds of long lists of operations
  87. * are often a bit ugly. A double for loop would be nicer and would
  88. * take up a lot less lines of code.
  89. * However, some compilers seem to fail in optimizing these kinds of
  90. * loops. So we will just have to do this by hand.
  91. */
  92. uint32_t a2[8] = {0};
  93. memcpy(a2, a, sizeof(uint32_t[8]));
  94. r[0] = a2[0] & b[0]; /* add (assignment, because r is 0) */
  95. r[1] = a2[1] & b[0];
  96. r[2] = a2[2] & b[0];
  97. r[3] = a2[3] & b[0];
  98. r[4] = a2[4] & b[0];
  99. r[5] = a2[5] & b[0];
  100. r[6] = a2[6] & b[0];
  101. r[7] = a2[7] & b[0];
  102. a2[0] ^= a2[7]; /* reduce */
  103. a2[2] ^= a2[7];
  104. a2[3] ^= a2[7];
  105. r[0] ^= a2[7] & b[1]; /* add */
  106. r[1] ^= a2[0] & b[1];
  107. r[2] ^= a2[1] & b[1];
  108. r[3] ^= a2[2] & b[1];
  109. r[4] ^= a2[3] & b[1];
  110. r[5] ^= a2[4] & b[1];
  111. r[6] ^= a2[5] & b[1];
  112. r[7] ^= a2[6] & b[1];
  113. a2[7] ^= a2[6]; /* reduce */
  114. a2[1] ^= a2[6];
  115. a2[2] ^= a2[6];
  116. r[0] ^= a2[6] & b[2]; /* add */
  117. r[1] ^= a2[7] & b[2];
  118. r[2] ^= a2[0] & b[2];
  119. r[3] ^= a2[1] & b[2];
  120. r[4] ^= a2[2] & b[2];
  121. r[5] ^= a2[3] & b[2];
  122. r[6] ^= a2[4] & b[2];
  123. r[7] ^= a2[5] & b[2];
  124. a2[6] ^= a2[5]; /* reduce */
  125. a2[0] ^= a2[5];
  126. a2[1] ^= a2[5];
  127. r[0] ^= a2[5] & b[3]; /* add */
  128. r[1] ^= a2[6] & b[3];
  129. r[2] ^= a2[7] & b[3];
  130. r[3] ^= a2[0] & b[3];
  131. r[4] ^= a2[1] & b[3];
  132. r[5] ^= a2[2] & b[3];
  133. r[6] ^= a2[3] & b[3];
  134. r[7] ^= a2[4] & b[3];
  135. a2[5] ^= a2[4]; /* reduce */
  136. a2[7] ^= a2[4];
  137. a2[0] ^= a2[4];
  138. r[0] ^= a2[4] & b[4]; /* add */
  139. r[1] ^= a2[5] & b[4];
  140. r[2] ^= a2[6] & b[4];
  141. r[3] ^= a2[7] & b[4];
  142. r[4] ^= a2[0] & b[4];
  143. r[5] ^= a2[1] & b[4];
  144. r[6] ^= a2[2] & b[4];
  145. r[7] ^= a2[3] & b[4];
  146. a2[4] ^= a2[3]; /* reduce */
  147. a2[6] ^= a2[3];
  148. a2[7] ^= a2[3];
  149. r[0] ^= a2[3] & b[5]; /* add */
  150. r[1] ^= a2[4] & b[5];
  151. r[2] ^= a2[5] & b[5];
  152. r[3] ^= a2[6] & b[5];
  153. r[4] ^= a2[7] & b[5];
  154. r[5] ^= a2[0] & b[5];
  155. r[6] ^= a2[1] & b[5];
  156. r[7] ^= a2[2] & b[5];
  157. a2[3] ^= a2[2]; /* reduce */
  158. a2[5] ^= a2[2];
  159. a2[6] ^= a2[2];
  160. r[0] ^= a2[2] & b[6]; /* add */
  161. r[1] ^= a2[3] & b[6];
  162. r[2] ^= a2[4] & b[6];
  163. r[3] ^= a2[5] & b[6];
  164. r[4] ^= a2[6] & b[6];
  165. r[5] ^= a2[7] & b[6];
  166. r[6] ^= a2[0] & b[6];
  167. r[7] ^= a2[1] & b[6];
  168. a2[2] ^= a2[1]; /* reduce */
  169. a2[4] ^= a2[1];
  170. a2[5] ^= a2[1];
  171. r[0] ^= a2[1] & b[7]; /* add */
  172. r[1] ^= a2[2] & b[7];
  173. r[2] ^= a2[3] & b[7];
  174. r[3] ^= a2[4] & b[7];
  175. r[4] ^= a2[5] & b[7];
  176. r[5] ^= a2[6] & b[7];
  177. r[6] ^= a2[7] & b[7];
  178. r[7] ^= a2[0] & b[7];
  179. memzero(a2, sizeof(a2));
  180. }
  181. /*
  182. * Square `x` in GF(2^8) and write the result to `r`. `r` and `x` may overlap.
  183. */
  184. static void gf256_square(uint32_t r[8], const uint32_t x[8]) {
  185. uint32_t r8 = 0, r10 = 0, r12 = 0, r14 = 0;
  186. /* Use the Freshman's Dream rule to square the polynomial
  187. * Assignments are done from 7 downto 0, because this allows the user
  188. * to execute this function in-place (e.g. `gf256_square(r, r);`).
  189. */
  190. r14 = x[7];
  191. r12 = x[6];
  192. r10 = x[5];
  193. r8 = x[4];
  194. r[6] = x[3];
  195. r[4] = x[2];
  196. r[2] = x[1];
  197. r[0] = x[0];
  198. /* Reduce with x^8 + x^4 + x^3 + x + 1 until order is less than 8 */
  199. r[7] = r14; /* r[7] was 0 */
  200. r[6] ^= r14;
  201. r10 ^= r14;
  202. /* Skip, because r13 is always 0 */
  203. r[4] ^= r12;
  204. r[5] = r12; /* r[5] was 0 */
  205. r[7] ^= r12;
  206. r8 ^= r12;
  207. /* Skip, because r11 is always 0 */
  208. r[2] ^= r10;
  209. r[3] = r10; /* r[3] was 0 */
  210. r[5] ^= r10;
  211. r[6] ^= r10;
  212. r[1] = r14; /* r[1] was 0 */
  213. r[2] ^= r14; /* Substitute r9 by r14 because they will always be equal*/
  214. r[4] ^= r14;
  215. r[5] ^= r14;
  216. r[0] ^= r8;
  217. r[1] ^= r8;
  218. r[3] ^= r8;
  219. r[4] ^= r8;
  220. }
  221. /*
  222. * Invert `x` in GF(2^8) and write the result to `r`
  223. */
  224. static void gf256_inv(uint32_t r[8], uint32_t x[8]) {
  225. uint32_t y[8] = {0}, z[8] = {0};
  226. gf256_square(y, x); // y = x^2
  227. gf256_square(y, y); // y = x^4
  228. gf256_square(r, y); // r = x^8
  229. gf256_mul(z, r, x); // z = x^9
  230. gf256_square(r, r); // r = x^16
  231. gf256_mul(r, r, z); // r = x^25
  232. gf256_square(r, r); // r = x^50
  233. gf256_square(z, r); // z = x^100
  234. gf256_square(z, z); // z = x^200
  235. gf256_mul(r, r, z); // r = x^250
  236. gf256_mul(r, r, y); // r = x^254
  237. memzero(y, sizeof(y));
  238. memzero(z, sizeof(z));
  239. }
  240. bool shamir_interpolate(
  241. uint8_t* result,
  242. uint8_t result_index,
  243. const uint8_t* share_indices,
  244. const uint8_t** share_values,
  245. uint8_t share_count,
  246. size_t len) {
  247. size_t i = 0, j = 0;
  248. uint32_t x[8] = {0};
  249. uint32_t xs[share_count][8];
  250. memset(xs, 0, sizeof(xs));
  251. uint32_t ys[share_count][8];
  252. memset(ys, 0, sizeof(ys));
  253. uint32_t num[8] = {~0}; /* num is the numerator (=1) */
  254. uint32_t denom[8] = {0};
  255. uint32_t tmp[8] = {0};
  256. uint32_t secret[8] = {0};
  257. bool ret = true;
  258. if(len > SHAMIR_MAX_LEN) return false;
  259. /* Collect the x and y values */
  260. for(i = 0; i < share_count; i++) {
  261. bitslice_setall(xs[i], share_indices[i]);
  262. bitslice(ys[i], share_values[i], len);
  263. }
  264. bitslice_setall(x, result_index);
  265. for(i = 0; i < share_count; i++) {
  266. memcpy(tmp, x, sizeof(uint32_t[8]));
  267. gf256_add(tmp, xs[i]);
  268. gf256_mul(num, num, tmp);
  269. }
  270. /* Use Lagrange basis polynomials to calculate the secret coefficient */
  271. for(i = 0; i < share_count; i++) {
  272. /* The code below assumes that none of the share_indices are equal to
  273. * result_index. We need to treat that as a special case. */
  274. if(share_indices[i] != result_index) {
  275. memcpy(denom, x, sizeof(denom));
  276. gf256_add(denom, xs[i]);
  277. } else {
  278. bitslice_setall(denom, 1);
  279. gf256_add(secret, ys[i]);
  280. }
  281. for(j = 0; j < share_count; j++) {
  282. if(i == j) continue;
  283. memcpy(tmp, xs[i], sizeof(uint32_t[8]));
  284. gf256_add(tmp, xs[j]);
  285. gf256_mul(denom, denom, tmp);
  286. }
  287. if((denom[0] | denom[1] | denom[2] | denom[3] | denom[4] | denom[5] | denom[6] |
  288. denom[7]) == 0) {
  289. /* The share_indices are not unique. */
  290. ret = false;
  291. break;
  292. }
  293. gf256_inv(tmp, denom); /* inverted denominator */
  294. gf256_mul(tmp, tmp, num); /* basis polynomial */
  295. gf256_mul(tmp, tmp, ys[i]); /* scaled coefficient */
  296. gf256_add(secret, tmp);
  297. }
  298. if(ret == true) {
  299. unbitslice(result, secret, len);
  300. }
  301. memzero(x, sizeof(x));
  302. memzero(xs, sizeof(xs));
  303. memzero(ys, sizeof(ys));
  304. memzero(num, sizeof(num));
  305. memzero(denom, sizeof(denom));
  306. memzero(tmp, sizeof(tmp));
  307. memzero(secret, sizeof(secret));
  308. return ret;
  309. }