You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

579 lines
22 KiB

  1. /*
  2. Public domain by Andrew M. <liquidsun@gmail.com>
  3. See: https://github.com/floodyberry/curve25519-donna
  4. 32 bit integer curve25519 implementation
  5. */
  6. typedef uint32_t bignum25519[10];
  7. typedef uint32_t bignum25519align16[12];
  8. static const uint32_t reduce_mask_25 = (1 << 25) - 1;
  9. static const uint32_t reduce_mask_26 = (1 << 26) - 1;
  10. /* out = in */
  11. DONNA_INLINE static void
  12. curve25519_copy(bignum25519 out, const bignum25519 in) {
  13. out[0] = in[0];
  14. out[1] = in[1];
  15. out[2] = in[2];
  16. out[3] = in[3];
  17. out[4] = in[4];
  18. out[5] = in[5];
  19. out[6] = in[6];
  20. out[7] = in[7];
  21. out[8] = in[8];
  22. out[9] = in[9];
  23. }
  24. /* out = a + b */
  25. DONNA_INLINE static void
  26. curve25519_add(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  27. out[0] = a[0] + b[0];
  28. out[1] = a[1] + b[1];
  29. out[2] = a[2] + b[2];
  30. out[3] = a[3] + b[3];
  31. out[4] = a[4] + b[4];
  32. out[5] = a[5] + b[5];
  33. out[6] = a[6] + b[6];
  34. out[7] = a[7] + b[7];
  35. out[8] = a[8] + b[8];
  36. out[9] = a[9] + b[9];
  37. }
  38. DONNA_INLINE static void
  39. curve25519_add_after_basic(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  40. uint32_t c;
  41. out[0] = a[0] + b[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  42. out[1] = a[1] + b[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  43. out[2] = a[2] + b[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  44. out[3] = a[3] + b[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  45. out[4] = a[4] + b[4] + c; c = (out[4] >> 26); out[4] &= reduce_mask_26;
  46. out[5] = a[5] + b[5] + c; c = (out[5] >> 25); out[5] &= reduce_mask_25;
  47. out[6] = a[6] + b[6] + c; c = (out[6] >> 26); out[6] &= reduce_mask_26;
  48. out[7] = a[7] + b[7] + c; c = (out[7] >> 25); out[7] &= reduce_mask_25;
  49. out[8] = a[8] + b[8] + c; c = (out[8] >> 26); out[8] &= reduce_mask_26;
  50. out[9] = a[9] + b[9] + c; c = (out[9] >> 25); out[9] &= reduce_mask_25;
  51. out[0] += 19 * c;
  52. }
  53. DONNA_INLINE static void
  54. curve25519_add_reduce(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  55. uint32_t c;
  56. out[0] = a[0] + b[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  57. out[1] = a[1] + b[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  58. out[2] = a[2] + b[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  59. out[3] = a[3] + b[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  60. out[4] = a[4] + b[4] + c; c = (out[4] >> 26); out[4] &= reduce_mask_26;
  61. out[5] = a[5] + b[5] + c; c = (out[5] >> 25); out[5] &= reduce_mask_25;
  62. out[6] = a[6] + b[6] + c; c = (out[6] >> 26); out[6] &= reduce_mask_26;
  63. out[7] = a[7] + b[7] + c; c = (out[7] >> 25); out[7] &= reduce_mask_25;
  64. out[8] = a[8] + b[8] + c; c = (out[8] >> 26); out[8] &= reduce_mask_26;
  65. out[9] = a[9] + b[9] + c; c = (out[9] >> 25); out[9] &= reduce_mask_25;
  66. out[0] += 19 * c;
  67. }
  68. /* multiples of p */
  69. static const uint32_t twoP0 = 0x07ffffda;
  70. static const uint32_t twoP13579 = 0x03fffffe;
  71. static const uint32_t twoP2468 = 0x07fffffe;
  72. static const uint32_t fourP0 = 0x0fffffb4;
  73. static const uint32_t fourP13579 = 0x07fffffc;
  74. static const uint32_t fourP2468 = 0x0ffffffc;
  75. /* out = a - b */
  76. DONNA_INLINE static void
  77. curve25519_sub(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  78. uint32_t c;
  79. out[0] = twoP0 + a[0] - b[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  80. out[1] = twoP13579 + a[1] - b[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  81. out[2] = twoP2468 + a[2] - b[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  82. out[3] = twoP13579 + a[3] - b[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  83. out[4] = twoP2468 + a[4] - b[4] + c;
  84. out[5] = twoP13579 + a[5] - b[5] ;
  85. out[6] = twoP2468 + a[6] - b[6] ;
  86. out[7] = twoP13579 + a[7] - b[7] ;
  87. out[8] = twoP2468 + a[8] - b[8] ;
  88. out[9] = twoP13579 + a[9] - b[9] ;
  89. }
  90. /* out = a - b, where a is the result of a basic op (add,sub) */
  91. DONNA_INLINE static void
  92. curve25519_sub_after_basic(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  93. uint32_t c;
  94. out[0] = fourP0 + a[0] - b[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  95. out[1] = fourP13579 + a[1] - b[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  96. out[2] = fourP2468 + a[2] - b[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  97. out[3] = fourP13579 + a[3] - b[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  98. out[4] = fourP2468 + a[4] - b[4] + c; c = (out[4] >> 26); out[4] &= reduce_mask_26;
  99. out[5] = fourP13579 + a[5] - b[5] + c; c = (out[5] >> 25); out[5] &= reduce_mask_25;
  100. out[6] = fourP2468 + a[6] - b[6] + c; c = (out[6] >> 26); out[6] &= reduce_mask_26;
  101. out[7] = fourP13579 + a[7] - b[7] + c; c = (out[7] >> 25); out[7] &= reduce_mask_25;
  102. out[8] = fourP2468 + a[8] - b[8] + c; c = (out[8] >> 26); out[8] &= reduce_mask_26;
  103. out[9] = fourP13579 + a[9] - b[9] + c; c = (out[9] >> 25); out[9] &= reduce_mask_25;
  104. out[0] += 19 * c;
  105. }
  106. DONNA_INLINE static void
  107. curve25519_sub_reduce(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  108. uint32_t c;
  109. out[0] = fourP0 + a[0] - b[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  110. out[1] = fourP13579 + a[1] - b[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  111. out[2] = fourP2468 + a[2] - b[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  112. out[3] = fourP13579 + a[3] - b[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  113. out[4] = fourP2468 + a[4] - b[4] + c; c = (out[4] >> 26); out[4] &= reduce_mask_26;
  114. out[5] = fourP13579 + a[5] - b[5] + c; c = (out[5] >> 25); out[5] &= reduce_mask_25;
  115. out[6] = fourP2468 + a[6] - b[6] + c; c = (out[6] >> 26); out[6] &= reduce_mask_26;
  116. out[7] = fourP13579 + a[7] - b[7] + c; c = (out[7] >> 25); out[7] &= reduce_mask_25;
  117. out[8] = fourP2468 + a[8] - b[8] + c; c = (out[8] >> 26); out[8] &= reduce_mask_26;
  118. out[9] = fourP13579 + a[9] - b[9] + c; c = (out[9] >> 25); out[9] &= reduce_mask_25;
  119. out[0] += 19 * c;
  120. }
  121. /* out = -a */
  122. DONNA_INLINE static void
  123. curve25519_neg(bignum25519 out, const bignum25519 a) {
  124. uint32_t c;
  125. out[0] = twoP0 - a[0] ; c = (out[0] >> 26); out[0] &= reduce_mask_26;
  126. out[1] = twoP13579 - a[1] + c; c = (out[1] >> 25); out[1] &= reduce_mask_25;
  127. out[2] = twoP2468 - a[2] + c; c = (out[2] >> 26); out[2] &= reduce_mask_26;
  128. out[3] = twoP13579 - a[3] + c; c = (out[3] >> 25); out[3] &= reduce_mask_25;
  129. out[4] = twoP2468 - a[4] + c; c = (out[4] >> 26); out[4] &= reduce_mask_26;
  130. out[5] = twoP13579 - a[5] + c; c = (out[5] >> 25); out[5] &= reduce_mask_25;
  131. out[6] = twoP2468 - a[6] + c; c = (out[6] >> 26); out[6] &= reduce_mask_26;
  132. out[7] = twoP13579 - a[7] + c; c = (out[7] >> 25); out[7] &= reduce_mask_25;
  133. out[8] = twoP2468 - a[8] + c; c = (out[8] >> 26); out[8] &= reduce_mask_26;
  134. out[9] = twoP13579 - a[9] + c; c = (out[9] >> 25); out[9] &= reduce_mask_25;
  135. out[0] += 19 * c;
  136. }
  137. /* out = a * b */
  138. #define curve25519_mul_noinline curve25519_mul
  139. static void
  140. curve25519_mul(bignum25519 out, const bignum25519 a, const bignum25519 b) {
  141. uint32_t r0,r1,r2,r3,r4,r5,r6,r7,r8,r9;
  142. uint32_t s0,s1,s2,s3,s4,s5,s6,s7,s8,s9;
  143. uint64_t m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,c;
  144. uint32_t p;
  145. r0 = b[0];
  146. r1 = b[1];
  147. r2 = b[2];
  148. r3 = b[3];
  149. r4 = b[4];
  150. r5 = b[5];
  151. r6 = b[6];
  152. r7 = b[7];
  153. r8 = b[8];
  154. r9 = b[9];
  155. s0 = a[0];
  156. s1 = a[1];
  157. s2 = a[2];
  158. s3 = a[3];
  159. s4 = a[4];
  160. s5 = a[5];
  161. s6 = a[6];
  162. s7 = a[7];
  163. s8 = a[8];
  164. s9 = a[9];
  165. m1 = mul32x32_64(r0, s1) + mul32x32_64(r1, s0);
  166. m3 = mul32x32_64(r0, s3) + mul32x32_64(r1, s2) + mul32x32_64(r2, s1) + mul32x32_64(r3, s0);
  167. m5 = mul32x32_64(r0, s5) + mul32x32_64(r1, s4) + mul32x32_64(r2, s3) + mul32x32_64(r3, s2) + mul32x32_64(r4, s1) + mul32x32_64(r5, s0);
  168. m7 = mul32x32_64(r0, s7) + mul32x32_64(r1, s6) + mul32x32_64(r2, s5) + mul32x32_64(r3, s4) + mul32x32_64(r4, s3) + mul32x32_64(r5, s2) + mul32x32_64(r6, s1) + mul32x32_64(r7, s0);
  169. m9 = mul32x32_64(r0, s9) + mul32x32_64(r1, s8) + mul32x32_64(r2, s7) + mul32x32_64(r3, s6) + mul32x32_64(r4, s5) + mul32x32_64(r5, s4) + mul32x32_64(r6, s3) + mul32x32_64(r7, s2) + mul32x32_64(r8, s1) + mul32x32_64(r9, s0);
  170. r1 *= 2;
  171. r3 *= 2;
  172. r5 *= 2;
  173. r7 *= 2;
  174. m0 = mul32x32_64(r0, s0);
  175. m2 = mul32x32_64(r0, s2) + mul32x32_64(r1, s1) + mul32x32_64(r2, s0);
  176. m4 = mul32x32_64(r0, s4) + mul32x32_64(r1, s3) + mul32x32_64(r2, s2) + mul32x32_64(r3, s1) + mul32x32_64(r4, s0);
  177. m6 = mul32x32_64(r0, s6) + mul32x32_64(r1, s5) + mul32x32_64(r2, s4) + mul32x32_64(r3, s3) + mul32x32_64(r4, s2) + mul32x32_64(r5, s1) + mul32x32_64(r6, s0);
  178. m8 = mul32x32_64(r0, s8) + mul32x32_64(r1, s7) + mul32x32_64(r2, s6) + mul32x32_64(r3, s5) + mul32x32_64(r4, s4) + mul32x32_64(r5, s3) + mul32x32_64(r6, s2) + mul32x32_64(r7, s1) + mul32x32_64(r8, s0);
  179. r1 *= 19;
  180. r2 *= 19;
  181. r3 = (r3 / 2) * 19;
  182. r4 *= 19;
  183. r5 = (r5 / 2) * 19;
  184. r6 *= 19;
  185. r7 = (r7 / 2) * 19;
  186. r8 *= 19;
  187. r9 *= 19;
  188. m1 += (mul32x32_64(r9, s2) + mul32x32_64(r8, s3) + mul32x32_64(r7, s4) + mul32x32_64(r6, s5) + mul32x32_64(r5, s6) + mul32x32_64(r4, s7) + mul32x32_64(r3, s8) + mul32x32_64(r2, s9));
  189. m3 += (mul32x32_64(r9, s4) + mul32x32_64(r8, s5) + mul32x32_64(r7, s6) + mul32x32_64(r6, s7) + mul32x32_64(r5, s8) + mul32x32_64(r4, s9));
  190. m5 += (mul32x32_64(r9, s6) + mul32x32_64(r8, s7) + mul32x32_64(r7, s8) + mul32x32_64(r6, s9));
  191. m7 += (mul32x32_64(r9, s8) + mul32x32_64(r8, s9));
  192. r3 *= 2;
  193. r5 *= 2;
  194. r7 *= 2;
  195. r9 *= 2;
  196. m0 += (mul32x32_64(r9, s1) + mul32x32_64(r8, s2) + mul32x32_64(r7, s3) + mul32x32_64(r6, s4) + mul32x32_64(r5, s5) + mul32x32_64(r4, s6) + mul32x32_64(r3, s7) + mul32x32_64(r2, s8) + mul32x32_64(r1, s9));
  197. m2 += (mul32x32_64(r9, s3) + mul32x32_64(r8, s4) + mul32x32_64(r7, s5) + mul32x32_64(r6, s6) + mul32x32_64(r5, s7) + mul32x32_64(r4, s8) + mul32x32_64(r3, s9));
  198. m4 += (mul32x32_64(r9, s5) + mul32x32_64(r8, s6) + mul32x32_64(r7, s7) + mul32x32_64(r6, s8) + mul32x32_64(r5, s9));
  199. m6 += (mul32x32_64(r9, s7) + mul32x32_64(r8, s8) + mul32x32_64(r7, s9));
  200. m8 += (mul32x32_64(r9, s9));
  201. r0 = (uint32_t)m0 & reduce_mask_26; c = (m0 >> 26);
  202. m1 += c; r1 = (uint32_t)m1 & reduce_mask_25; c = (m1 >> 25);
  203. m2 += c; r2 = (uint32_t)m2 & reduce_mask_26; c = (m2 >> 26);
  204. m3 += c; r3 = (uint32_t)m3 & reduce_mask_25; c = (m3 >> 25);
  205. m4 += c; r4 = (uint32_t)m4 & reduce_mask_26; c = (m4 >> 26);
  206. m5 += c; r5 = (uint32_t)m5 & reduce_mask_25; c = (m5 >> 25);
  207. m6 += c; r6 = (uint32_t)m6 & reduce_mask_26; c = (m6 >> 26);
  208. m7 += c; r7 = (uint32_t)m7 & reduce_mask_25; c = (m7 >> 25);
  209. m8 += c; r8 = (uint32_t)m8 & reduce_mask_26; c = (m8 >> 26);
  210. m9 += c; r9 = (uint32_t)m9 & reduce_mask_25; p = (uint32_t)(m9 >> 25);
  211. m0 = r0 + mul32x32_64(p,19); r0 = (uint32_t)m0 & reduce_mask_26; p = (uint32_t)(m0 >> 26);
  212. r1 += p;
  213. out[0] = r0;
  214. out[1] = r1;
  215. out[2] = r2;
  216. out[3] = r3;
  217. out[4] = r4;
  218. out[5] = r5;
  219. out[6] = r6;
  220. out[7] = r7;
  221. out[8] = r8;
  222. out[9] = r9;
  223. }
  224. /* out = in*in */
  225. static void
  226. curve25519_square(bignum25519 out, const bignum25519 in) {
  227. uint32_t r0,r1,r2,r3,r4,r5,r6,r7,r8,r9;
  228. uint32_t d6,d7,d8,d9;
  229. uint64_t m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,c;
  230. uint32_t p;
  231. r0 = in[0];
  232. r1 = in[1];
  233. r2 = in[2];
  234. r3 = in[3];
  235. r4 = in[4];
  236. r5 = in[5];
  237. r6 = in[6];
  238. r7 = in[7];
  239. r8 = in[8];
  240. r9 = in[9];
  241. m0 = mul32x32_64(r0, r0);
  242. r0 *= 2;
  243. m1 = mul32x32_64(r0, r1);
  244. m2 = mul32x32_64(r0, r2) + mul32x32_64(r1, r1 * 2);
  245. r1 *= 2;
  246. m3 = mul32x32_64(r0, r3) + mul32x32_64(r1, r2 );
  247. m4 = mul32x32_64(r0, r4) + mul32x32_64(r1, r3 * 2) + mul32x32_64(r2, r2);
  248. r2 *= 2;
  249. m5 = mul32x32_64(r0, r5) + mul32x32_64(r1, r4 ) + mul32x32_64(r2, r3);
  250. m6 = mul32x32_64(r0, r6) + mul32x32_64(r1, r5 * 2) + mul32x32_64(r2, r4) + mul32x32_64(r3, r3 * 2);
  251. r3 *= 2;
  252. m7 = mul32x32_64(r0, r7) + mul32x32_64(r1, r6 ) + mul32x32_64(r2, r5) + mul32x32_64(r3, r4 );
  253. m8 = mul32x32_64(r0, r8) + mul32x32_64(r1, r7 * 2) + mul32x32_64(r2, r6) + mul32x32_64(r3, r5 * 2) + mul32x32_64(r4, r4 );
  254. m9 = mul32x32_64(r0, r9) + mul32x32_64(r1, r8 ) + mul32x32_64(r2, r7) + mul32x32_64(r3, r6 ) + mul32x32_64(r4, r5 * 2);
  255. d6 = r6 * 19;
  256. d7 = r7 * 2 * 19;
  257. d8 = r8 * 19;
  258. d9 = r9 * 2 * 19;
  259. m0 += (mul32x32_64(d9, r1 ) + mul32x32_64(d8, r2 ) + mul32x32_64(d7, r3 ) + mul32x32_64(d6, r4 * 2) + mul32x32_64(r5, r5 * 2 * 19));
  260. m1 += (mul32x32_64(d9, r2 / 2) + mul32x32_64(d8, r3 ) + mul32x32_64(d7, r4 ) + mul32x32_64(d6, r5 * 2));
  261. m2 += (mul32x32_64(d9, r3 ) + mul32x32_64(d8, r4 * 2) + mul32x32_64(d7, r5 * 2) + mul32x32_64(d6, r6 ));
  262. m3 += (mul32x32_64(d9, r4 ) + mul32x32_64(d8, r5 * 2) + mul32x32_64(d7, r6 ));
  263. m4 += (mul32x32_64(d9, r5 * 2) + mul32x32_64(d8, r6 * 2) + mul32x32_64(d7, r7 ));
  264. m5 += (mul32x32_64(d9, r6 ) + mul32x32_64(d8, r7 * 2));
  265. m6 += (mul32x32_64(d9, r7 * 2) + mul32x32_64(d8, r8 ));
  266. m7 += (mul32x32_64(d9, r8 ));
  267. m8 += (mul32x32_64(d9, r9 ));
  268. r0 = (uint32_t)m0 & reduce_mask_26; c = (m0 >> 26);
  269. m1 += c; r1 = (uint32_t)m1 & reduce_mask_25; c = (m1 >> 25);
  270. m2 += c; r2 = (uint32_t)m2 & reduce_mask_26; c = (m2 >> 26);
  271. m3 += c; r3 = (uint32_t)m3 & reduce_mask_25; c = (m3 >> 25);
  272. m4 += c; r4 = (uint32_t)m4 & reduce_mask_26; c = (m4 >> 26);
  273. m5 += c; r5 = (uint32_t)m5 & reduce_mask_25; c = (m5 >> 25);
  274. m6 += c; r6 = (uint32_t)m6 & reduce_mask_26; c = (m6 >> 26);
  275. m7 += c; r7 = (uint32_t)m7 & reduce_mask_25; c = (m7 >> 25);
  276. m8 += c; r8 = (uint32_t)m8 & reduce_mask_26; c = (m8 >> 26);
  277. m9 += c; r9 = (uint32_t)m9 & reduce_mask_25; p = (uint32_t)(m9 >> 25);
  278. m0 = r0 + mul32x32_64(p,19); r0 = (uint32_t)m0 & reduce_mask_26; p = (uint32_t)(m0 >> 26);
  279. r1 += p;
  280. out[0] = r0;
  281. out[1] = r1;
  282. out[2] = r2;
  283. out[3] = r3;
  284. out[4] = r4;
  285. out[5] = r5;
  286. out[6] = r6;
  287. out[7] = r7;
  288. out[8] = r8;
  289. out[9] = r9;
  290. }
  291. /* out = in ^ (2 * count) */
  292. static void
  293. curve25519_square_times(bignum25519 out, const bignum25519 in, int count) {
  294. uint32_t r0,r1,r2,r3,r4,r5,r6,r7,r8,r9;
  295. uint32_t d6,d7,d8,d9;
  296. uint64_t m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,c;
  297. uint32_t p;
  298. r0 = in[0];
  299. r1 = in[1];
  300. r2 = in[2];
  301. r3 = in[3];
  302. r4 = in[4];
  303. r5 = in[5];
  304. r6 = in[6];
  305. r7 = in[7];
  306. r8 = in[8];
  307. r9 = in[9];
  308. do {
  309. m0 = mul32x32_64(r0, r0);
  310. r0 *= 2;
  311. m1 = mul32x32_64(r0, r1);
  312. m2 = mul32x32_64(r0, r2) + mul32x32_64(r1, r1 * 2);
  313. r1 *= 2;
  314. m3 = mul32x32_64(r0, r3) + mul32x32_64(r1, r2 );
  315. m4 = mul32x32_64(r0, r4) + mul32x32_64(r1, r3 * 2) + mul32x32_64(r2, r2);
  316. r2 *= 2;
  317. m5 = mul32x32_64(r0, r5) + mul32x32_64(r1, r4 ) + mul32x32_64(r2, r3);
  318. m6 = mul32x32_64(r0, r6) + mul32x32_64(r1, r5 * 2) + mul32x32_64(r2, r4) + mul32x32_64(r3, r3 * 2);
  319. r3 *= 2;
  320. m7 = mul32x32_64(r0, r7) + mul32x32_64(r1, r6 ) + mul32x32_64(r2, r5) + mul32x32_64(r3, r4 );
  321. m8 = mul32x32_64(r0, r8) + mul32x32_64(r1, r7 * 2) + mul32x32_64(r2, r6) + mul32x32_64(r3, r5 * 2) + mul32x32_64(r4, r4 );
  322. m9 = mul32x32_64(r0, r9) + mul32x32_64(r1, r8 ) + mul32x32_64(r2, r7) + mul32x32_64(r3, r6 ) + mul32x32_64(r4, r5 * 2);
  323. d6 = r6 * 19;
  324. d7 = r7 * 2 * 19;
  325. d8 = r8 * 19;
  326. d9 = r9 * 2 * 19;
  327. m0 += (mul32x32_64(d9, r1 ) + mul32x32_64(d8, r2 ) + mul32x32_64(d7, r3 ) + mul32x32_64(d6, r4 * 2) + mul32x32_64(r5, r5 * 2 * 19));
  328. m1 += (mul32x32_64(d9, r2 / 2) + mul32x32_64(d8, r3 ) + mul32x32_64(d7, r4 ) + mul32x32_64(d6, r5 * 2));
  329. m2 += (mul32x32_64(d9, r3 ) + mul32x32_64(d8, r4 * 2) + mul32x32_64(d7, r5 * 2) + mul32x32_64(d6, r6 ));
  330. m3 += (mul32x32_64(d9, r4 ) + mul32x32_64(d8, r5 * 2) + mul32x32_64(d7, r6 ));
  331. m4 += (mul32x32_64(d9, r5 * 2) + mul32x32_64(d8, r6 * 2) + mul32x32_64(d7, r7 ));
  332. m5 += (mul32x32_64(d9, r6 ) + mul32x32_64(d8, r7 * 2));
  333. m6 += (mul32x32_64(d9, r7 * 2) + mul32x32_64(d8, r8 ));
  334. m7 += (mul32x32_64(d9, r8 ));
  335. m8 += (mul32x32_64(d9, r9 ));
  336. r0 = (uint32_t)m0 & reduce_mask_26; c = (m0 >> 26);
  337. m1 += c; r1 = (uint32_t)m1 & reduce_mask_25; c = (m1 >> 25);
  338. m2 += c; r2 = (uint32_t)m2 & reduce_mask_26; c = (m2 >> 26);
  339. m3 += c; r3 = (uint32_t)m3 & reduce_mask_25; c = (m3 >> 25);
  340. m4 += c; r4 = (uint32_t)m4 & reduce_mask_26; c = (m4 >> 26);
  341. m5 += c; r5 = (uint32_t)m5 & reduce_mask_25; c = (m5 >> 25);
  342. m6 += c; r6 = (uint32_t)m6 & reduce_mask_26; c = (m6 >> 26);
  343. m7 += c; r7 = (uint32_t)m7 & reduce_mask_25; c = (m7 >> 25);
  344. m8 += c; r8 = (uint32_t)m8 & reduce_mask_26; c = (m8 >> 26);
  345. m9 += c; r9 = (uint32_t)m9 & reduce_mask_25; p = (uint32_t)(m9 >> 25);
  346. m0 = r0 + mul32x32_64(p,19); r0 = (uint32_t)m0 & reduce_mask_26; p = (uint32_t)(m0 >> 26);
  347. r1 += p;
  348. } while (--count);
  349. out[0] = r0;
  350. out[1] = r1;
  351. out[2] = r2;
  352. out[3] = r3;
  353. out[4] = r4;
  354. out[5] = r5;
  355. out[6] = r6;
  356. out[7] = r7;
  357. out[8] = r8;
  358. out[9] = r9;
  359. }
  360. /* Take a little-endian, 32-byte number and expand it into polynomial form */
  361. static void
  362. curve25519_expand(bignum25519 out, const unsigned char in[32]) {
  363. static const union { uint8_t b[2]; uint16_t s; } endian_check = {{1,0}};
  364. uint32_t x0,x1,x2,x3,x4,x5,x6,x7;
  365. if (endian_check.s == 1) {
  366. x0 = *(uint32_t *)(in + 0);
  367. x1 = *(uint32_t *)(in + 4);
  368. x2 = *(uint32_t *)(in + 8);
  369. x3 = *(uint32_t *)(in + 12);
  370. x4 = *(uint32_t *)(in + 16);
  371. x5 = *(uint32_t *)(in + 20);
  372. x6 = *(uint32_t *)(in + 24);
  373. x7 = *(uint32_t *)(in + 28);
  374. } else {
  375. #define F(s) \
  376. ((((uint32_t)in[s + 0]) ) | \
  377. (((uint32_t)in[s + 1]) << 8) | \
  378. (((uint32_t)in[s + 2]) << 16) | \
  379. (((uint32_t)in[s + 3]) << 24))
  380. x0 = F(0);
  381. x1 = F(4);
  382. x2 = F(8);
  383. x3 = F(12);
  384. x4 = F(16);
  385. x5 = F(20);
  386. x6 = F(24);
  387. x7 = F(28);
  388. #undef F
  389. }
  390. out[0] = ( x0 ) & 0x3ffffff;
  391. out[1] = ((((uint64_t)x1 << 32) | x0) >> 26) & 0x1ffffff;
  392. out[2] = ((((uint64_t)x2 << 32) | x1) >> 19) & 0x3ffffff;
  393. out[3] = ((((uint64_t)x3 << 32) | x2) >> 13) & 0x1ffffff;
  394. out[4] = (( x3) >> 6) & 0x3ffffff;
  395. out[5] = ( x4 ) & 0x1ffffff;
  396. out[6] = ((((uint64_t)x5 << 32) | x4) >> 25) & 0x3ffffff;
  397. out[7] = ((((uint64_t)x6 << 32) | x5) >> 19) & 0x1ffffff;
  398. out[8] = ((((uint64_t)x7 << 32) | x6) >> 12) & 0x3ffffff;
  399. out[9] = (( x7) >> 6) & 0x1ffffff;
  400. }
  401. /* Take a fully reduced polynomial form number and contract it into a
  402. * little-endian, 32-byte array
  403. */
  404. static void
  405. curve25519_contract(unsigned char out[32], const bignum25519 in) {
  406. bignum25519 f;
  407. curve25519_copy(f, in);
  408. #define carry_pass() \
  409. f[1] += f[0] >> 26; f[0] &= reduce_mask_26; \
  410. f[2] += f[1] >> 25; f[1] &= reduce_mask_25; \
  411. f[3] += f[2] >> 26; f[2] &= reduce_mask_26; \
  412. f[4] += f[3] >> 25; f[3] &= reduce_mask_25; \
  413. f[5] += f[4] >> 26; f[4] &= reduce_mask_26; \
  414. f[6] += f[5] >> 25; f[5] &= reduce_mask_25; \
  415. f[7] += f[6] >> 26; f[6] &= reduce_mask_26; \
  416. f[8] += f[7] >> 25; f[7] &= reduce_mask_25; \
  417. f[9] += f[8] >> 26; f[8] &= reduce_mask_26;
  418. #define carry_pass_full() \
  419. carry_pass() \
  420. f[0] += 19 * (f[9] >> 25); f[9] &= reduce_mask_25;
  421. #define carry_pass_final() \
  422. carry_pass() \
  423. f[9] &= reduce_mask_25;
  424. carry_pass_full()
  425. carry_pass_full()
  426. /* now t is between 0 and 2^255-1, properly carried. */
  427. /* case 1: between 0 and 2^255-20. case 2: between 2^255-19 and 2^255-1. */
  428. f[0] += 19;
  429. carry_pass_full()
  430. /* now between 19 and 2^255-1 in both cases, and offset by 19. */
  431. f[0] += (reduce_mask_26 + 1) - 19;
  432. f[1] += (reduce_mask_25 + 1) - 1;
  433. f[2] += (reduce_mask_26 + 1) - 1;
  434. f[3] += (reduce_mask_25 + 1) - 1;
  435. f[4] += (reduce_mask_26 + 1) - 1;
  436. f[5] += (reduce_mask_25 + 1) - 1;
  437. f[6] += (reduce_mask_26 + 1) - 1;
  438. f[7] += (reduce_mask_25 + 1) - 1;
  439. f[8] += (reduce_mask_26 + 1) - 1;
  440. f[9] += (reduce_mask_25 + 1) - 1;
  441. /* now between 2^255 and 2^256-20, and offset by 2^255. */
  442. carry_pass_final()
  443. #undef carry_pass
  444. #undef carry_full
  445. #undef carry_final
  446. f[1] <<= 2;
  447. f[2] <<= 3;
  448. f[3] <<= 5;
  449. f[4] <<= 6;
  450. f[6] <<= 1;
  451. f[7] <<= 3;
  452. f[8] <<= 4;
  453. f[9] <<= 6;
  454. #define F(i, s) \
  455. out[s+0] |= (unsigned char )(f[i] & 0xff); \
  456. out[s+1] = (unsigned char )((f[i] >> 8) & 0xff); \
  457. out[s+2] = (unsigned char )((f[i] >> 16) & 0xff); \
  458. out[s+3] = (unsigned char )((f[i] >> 24) & 0xff);
  459. out[0] = 0;
  460. out[16] = 0;
  461. F(0,0);
  462. F(1,3);
  463. F(2,6);
  464. F(3,9);
  465. F(4,12);
  466. F(5,16);
  467. F(6,19);
  468. F(7,22);
  469. F(8,25);
  470. F(9,28);
  471. #undef F
  472. }
  473. /* out = (flag) ? in : out */
  474. DONNA_INLINE static void
  475. curve25519_move_conditional_bytes(uint8_t out[96], const uint8_t in[96], uint32_t flag) {
  476. const uint32_t nb = flag - 1, b = ~nb;
  477. const uint32_t *inl = (const uint32_t *)in;
  478. uint32_t *outl = (uint32_t *)out;
  479. outl[0] = (outl[0] & nb) | (inl[0] & b);
  480. outl[1] = (outl[1] & nb) | (inl[1] & b);
  481. outl[2] = (outl[2] & nb) | (inl[2] & b);
  482. outl[3] = (outl[3] & nb) | (inl[3] & b);
  483. outl[4] = (outl[4] & nb) | (inl[4] & b);
  484. outl[5] = (outl[5] & nb) | (inl[5] & b);
  485. outl[6] = (outl[6] & nb) | (inl[6] & b);
  486. outl[7] = (outl[7] & nb) | (inl[7] & b);
  487. outl[8] = (outl[8] & nb) | (inl[8] & b);
  488. outl[9] = (outl[9] & nb) | (inl[9] & b);
  489. outl[10] = (outl[10] & nb) | (inl[10] & b);
  490. outl[11] = (outl[11] & nb) | (inl[11] & b);
  491. outl[12] = (outl[12] & nb) | (inl[12] & b);
  492. outl[13] = (outl[13] & nb) | (inl[13] & b);
  493. outl[14] = (outl[14] & nb) | (inl[14] & b);
  494. outl[15] = (outl[15] & nb) | (inl[15] & b);
  495. outl[16] = (outl[16] & nb) | (inl[16] & b);
  496. outl[17] = (outl[17] & nb) | (inl[17] & b);
  497. outl[18] = (outl[18] & nb) | (inl[18] & b);
  498. outl[19] = (outl[19] & nb) | (inl[19] & b);
  499. outl[20] = (outl[20] & nb) | (inl[20] & b);
  500. outl[21] = (outl[21] & nb) | (inl[21] & b);
  501. outl[22] = (outl[22] & nb) | (inl[22] & b);
  502. outl[23] = (outl[23] & nb) | (inl[23] & b);
  503. }
  504. /* if (iswap) swap(a, b) */
  505. DONNA_INLINE static void
  506. curve25519_swap_conditional(bignum25519 a, bignum25519 b, uint32_t iswap) {
  507. const uint32_t swap = (uint32_t)(-(int32_t)iswap);
  508. uint32_t x0,x1,x2,x3,x4,x5,x6,x7,x8,x9;
  509. x0 = swap & (a[0] ^ b[0]); a[0] ^= x0; b[0] ^= x0;
  510. x1 = swap & (a[1] ^ b[1]); a[1] ^= x1; b[1] ^= x1;
  511. x2 = swap & (a[2] ^ b[2]); a[2] ^= x2; b[2] ^= x2;
  512. x3 = swap & (a[3] ^ b[3]); a[3] ^= x3; b[3] ^= x3;
  513. x4 = swap & (a[4] ^ b[4]); a[4] ^= x4; b[4] ^= x4;
  514. x5 = swap & (a[5] ^ b[5]); a[5] ^= x5; b[5] ^= x5;
  515. x6 = swap & (a[6] ^ b[6]); a[6] ^= x6; b[6] ^= x6;
  516. x7 = swap & (a[7] ^ b[7]); a[7] ^= x7; b[7] ^= x7;
  517. x8 = swap & (a[8] ^ b[8]); a[8] ^= x8; b[8] ^= x8;
  518. x9 = swap & (a[9] ^ b[9]); a[9] ^= x9; b[9] ^= x9;
  519. }