curve25519-donna-c64.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451
  1. /* Copyright 2008, Google Inc.
  2. * All rights reserved.
  3. *
  4. * Code released into the public domain.
  5. *
  6. * curve25519-donna: Curve25519 elliptic curve, public key function
  7. *
  8. * http://code.google.com/p/curve25519-donna/
  9. *
  10. * Adam Langley <agl@imperialviolet.org>
  11. *
  12. * Derived from public domain C code by Daniel J. Bernstein <djb@cr.yp.to>
  13. *
  14. * More information about curve25519 can be found here
  15. * http://cr.yp.to/ecdh.html
  16. *
  17. * djb's sample implementation of curve25519 is written in a special assembly
  18. * language called qhasm and uses the floating point registers.
  19. *
  20. * This is, almost, a clean room reimplementation from the curve25519 paper. It
  21. * uses many of the tricks described therein. Only the crecip function is taken
  22. * from the sample implementation.
  23. */
  24. #include "orconfig.h"
  25. #include <string.h>
  26. #include "lib/cc/torint.h"
  27. typedef uint8_t u8;
  28. typedef uint64_t limb;
  29. typedef limb felem[5];
  30. // This is a special gcc mode for 128-bit integers. It's implemented on 64-bit
  31. // platforms only as far as I know.
  32. typedef unsigned uint128_t __attribute__((mode(TI)));
  33. #undef force_inline
  34. #define force_inline __attribute__((always_inline))
  35. /* Sum two numbers: output += in */
  36. static inline void force_inline
  37. fsum(limb *output, const limb *in) {
  38. output[0] += in[0];
  39. output[1] += in[1];
  40. output[2] += in[2];
  41. output[3] += in[3];
  42. output[4] += in[4];
  43. }
  44. /* Find the difference of two numbers: output = in - output
  45. * (note the order of the arguments!)
  46. *
  47. * Assumes that out[i] < 2**52
  48. * On return, out[i] < 2**55
  49. */
  50. static inline void force_inline
  51. fdifference_backwards(felem out, const felem in) {
  52. /* 152 is 19 << 3 */
  53. static const limb two54m152 = (((limb)1) << 54) - 152;
  54. static const limb two54m8 = (((limb)1) << 54) - 8;
  55. out[0] = in[0] + two54m152 - out[0];
  56. out[1] = in[1] + two54m8 - out[1];
  57. out[2] = in[2] + two54m8 - out[2];
  58. out[3] = in[3] + two54m8 - out[3];
  59. out[4] = in[4] + two54m8 - out[4];
  60. }
  61. /* Multiply a number by a scalar: output = in * scalar */
  62. static inline void force_inline
  63. fscalar_product(felem output, const felem in, const limb scalar) {
  64. uint128_t a;
  65. a = ((uint128_t) in[0]) * scalar;
  66. output[0] = ((limb)a) & 0x7ffffffffffff;
  67. a = ((uint128_t) in[1]) * scalar + ((limb) (a >> 51));
  68. output[1] = ((limb)a) & 0x7ffffffffffff;
  69. a = ((uint128_t) in[2]) * scalar + ((limb) (a >> 51));
  70. output[2] = ((limb)a) & 0x7ffffffffffff;
  71. a = ((uint128_t) in[3]) * scalar + ((limb) (a >> 51));
  72. output[3] = ((limb)a) & 0x7ffffffffffff;
  73. a = ((uint128_t) in[4]) * scalar + ((limb) (a >> 51));
  74. output[4] = ((limb)a) & 0x7ffffffffffff;
  75. output[0] += (a >> 51) * 19;
  76. }
  77. /* Multiply two numbers: output = in2 * in
  78. *
  79. * output must be distinct to both inputs. The inputs are reduced coefficient
  80. * form, the output is not.
  81. *
  82. * Assumes that in[i] < 2**55 and likewise for in2.
  83. * On return, output[i] < 2**52
  84. */
  85. static inline void force_inline
  86. fmul(felem output, const felem in2, const felem in) {
  87. uint128_t t[5];
  88. limb r0,r1,r2,r3,r4,s0,s1,s2,s3,s4,c;
  89. r0 = in[0];
  90. r1 = in[1];
  91. r2 = in[2];
  92. r3 = in[3];
  93. r4 = in[4];
  94. s0 = in2[0];
  95. s1 = in2[1];
  96. s2 = in2[2];
  97. s3 = in2[3];
  98. s4 = in2[4];
  99. t[0] = ((uint128_t) r0) * s0;
  100. t[1] = ((uint128_t) r0) * s1 + ((uint128_t) r1) * s0;
  101. t[2] = ((uint128_t) r0) * s2 + ((uint128_t) r2) * s0 + ((uint128_t) r1) * s1;
  102. t[3] = ((uint128_t) r0) * s3 + ((uint128_t) r3) * s0 + ((uint128_t) r1) * s2 + ((uint128_t) r2) * s1;
  103. t[4] = ((uint128_t) r0) * s4 + ((uint128_t) r4) * s0 + ((uint128_t) r3) * s1 + ((uint128_t) r1) * s3 + ((uint128_t) r2) * s2;
  104. r4 *= 19;
  105. r1 *= 19;
  106. r2 *= 19;
  107. r3 *= 19;
  108. t[0] += ((uint128_t) r4) * s1 + ((uint128_t) r1) * s4 + ((uint128_t) r2) * s3 + ((uint128_t) r3) * s2;
  109. t[1] += ((uint128_t) r4) * s2 + ((uint128_t) r2) * s4 + ((uint128_t) r3) * s3;
  110. t[2] += ((uint128_t) r4) * s3 + ((uint128_t) r3) * s4;
  111. t[3] += ((uint128_t) r4) * s4;
  112. r0 = (limb)t[0] & 0x7ffffffffffff; c = (limb)(t[0] >> 51);
  113. t[1] += c; r1 = (limb)t[1] & 0x7ffffffffffff; c = (limb)(t[1] >> 51);
  114. t[2] += c; r2 = (limb)t[2] & 0x7ffffffffffff; c = (limb)(t[2] >> 51);
  115. t[3] += c; r3 = (limb)t[3] & 0x7ffffffffffff; c = (limb)(t[3] >> 51);
  116. t[4] += c; r4 = (limb)t[4] & 0x7ffffffffffff; c = (limb)(t[4] >> 51);
  117. r0 += c * 19; c = r0 >> 51; r0 = r0 & 0x7ffffffffffff;
  118. r1 += c; c = r1 >> 51; r1 = r1 & 0x7ffffffffffff;
  119. r2 += c;
  120. output[0] = r0;
  121. output[1] = r1;
  122. output[2] = r2;
  123. output[3] = r3;
  124. output[4] = r4;
  125. }
  126. static inline void force_inline
  127. fsquare_times(felem output, const felem in, limb count) {
  128. uint128_t t[5];
  129. limb r0,r1,r2,r3,r4,c;
  130. limb d0,d1,d2,d4,d419;
  131. r0 = in[0];
  132. r1 = in[1];
  133. r2 = in[2];
  134. r3 = in[3];
  135. r4 = in[4];
  136. do {
  137. d0 = r0 * 2;
  138. d1 = r1 * 2;
  139. d2 = r2 * 2 * 19;
  140. d419 = r4 * 19;
  141. d4 = d419 * 2;
  142. t[0] = ((uint128_t) r0) * r0 + ((uint128_t) d4) * r1 + (((uint128_t) d2) * (r3 ));
  143. t[1] = ((uint128_t) d0) * r1 + ((uint128_t) d4) * r2 + (((uint128_t) r3) * (r3 * 19));
  144. t[2] = ((uint128_t) d0) * r2 + ((uint128_t) r1) * r1 + (((uint128_t) d4) * (r3 ));
  145. t[3] = ((uint128_t) d0) * r3 + ((uint128_t) d1) * r2 + (((uint128_t) r4) * (d419 ));
  146. t[4] = ((uint128_t) d0) * r4 + ((uint128_t) d1) * r3 + (((uint128_t) r2) * (r2 ));
  147. r0 = (limb)t[0] & 0x7ffffffffffff; c = (limb)(t[0] >> 51);
  148. t[1] += c; r1 = (limb)t[1] & 0x7ffffffffffff; c = (limb)(t[1] >> 51);
  149. t[2] += c; r2 = (limb)t[2] & 0x7ffffffffffff; c = (limb)(t[2] >> 51);
  150. t[3] += c; r3 = (limb)t[3] & 0x7ffffffffffff; c = (limb)(t[3] >> 51);
  151. t[4] += c; r4 = (limb)t[4] & 0x7ffffffffffff; c = (limb)(t[4] >> 51);
  152. r0 += c * 19; c = r0 >> 51; r0 = r0 & 0x7ffffffffffff;
  153. r1 += c; c = r1 >> 51; r1 = r1 & 0x7ffffffffffff;
  154. r2 += c;
  155. } while(--count);
  156. output[0] = r0;
  157. output[1] = r1;
  158. output[2] = r2;
  159. output[3] = r3;
  160. output[4] = r4;
  161. }
  162. /* Load a little-endian 64-bit number */
  163. static limb
  164. load_limb(const u8 *in) {
  165. return
  166. ((limb)in[0]) |
  167. (((limb)in[1]) << 8) |
  168. (((limb)in[2]) << 16) |
  169. (((limb)in[3]) << 24) |
  170. (((limb)in[4]) << 32) |
  171. (((limb)in[5]) << 40) |
  172. (((limb)in[6]) << 48) |
  173. (((limb)in[7]) << 56);
  174. }
  175. static void
  176. store_limb(u8 *out, limb in) {
  177. out[0] = in & 0xff;
  178. out[1] = (in >> 8) & 0xff;
  179. out[2] = (in >> 16) & 0xff;
  180. out[3] = (in >> 24) & 0xff;
  181. out[4] = (in >> 32) & 0xff;
  182. out[5] = (in >> 40) & 0xff;
  183. out[6] = (in >> 48) & 0xff;
  184. out[7] = (in >> 56) & 0xff;
  185. }
  186. /* Take a little-endian, 32-byte number and expand it into polynomial form */
  187. static void
  188. fexpand(limb *output, const u8 *in) {
  189. output[0] = load_limb(in) & 0x7ffffffffffff;
  190. output[1] = (load_limb(in+6) >> 3) & 0x7ffffffffffff;
  191. output[2] = (load_limb(in+12) >> 6) & 0x7ffffffffffff;
  192. output[3] = (load_limb(in+19) >> 1) & 0x7ffffffffffff;
  193. output[4] = (load_limb(in+24) >> 12) & 0x7ffffffffffff;
  194. }
  195. /* Take a fully reduced polynomial form number and contract it into a
  196. * little-endian, 32-byte array
  197. */
  198. static void
  199. fcontract(u8 *output, const felem input) {
  200. uint128_t t[5];
  201. t[0] = input[0];
  202. t[1] = input[1];
  203. t[2] = input[2];
  204. t[3] = input[3];
  205. t[4] = input[4];
  206. t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
  207. t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
  208. t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
  209. t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
  210. t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;
  211. t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
  212. t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
  213. t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
  214. t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
  215. t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;
  216. /* now t is between 0 and 2^255-1, properly carried. */
  217. /* case 1: between 0 and 2^255-20. case 2: between 2^255-19 and 2^255-1. */
  218. t[0] += 19;
  219. t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
  220. t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
  221. t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
  222. t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
  223. t[0] += 19 * (t[4] >> 51); t[4] &= 0x7ffffffffffff;
  224. /* now between 19 and 2^255-1 in both cases, and offset by 19. */
  225. t[0] += 0x8000000000000 - 19;
  226. t[1] += 0x8000000000000 - 1;
  227. t[2] += 0x8000000000000 - 1;
  228. t[3] += 0x8000000000000 - 1;
  229. t[4] += 0x8000000000000 - 1;
  230. /* now between 2^255 and 2^256-20, and offset by 2^255. */
  231. t[1] += t[0] >> 51; t[0] &= 0x7ffffffffffff;
  232. t[2] += t[1] >> 51; t[1] &= 0x7ffffffffffff;
  233. t[3] += t[2] >> 51; t[2] &= 0x7ffffffffffff;
  234. t[4] += t[3] >> 51; t[3] &= 0x7ffffffffffff;
  235. t[4] &= 0x7ffffffffffff;
  236. store_limb(output, t[0] | (t[1] << 51));
  237. store_limb(output+8, (t[1] >> 13) | (t[2] << 38));
  238. store_limb(output+16, (t[2] >> 26) | (t[3] << 25));
  239. store_limb(output+24, (t[3] >> 39) | (t[4] << 12));
  240. }
  241. /* Input: Q, Q', Q-Q'
  242. * Output: 2Q, Q+Q'
  243. *
  244. * x2 z3: long form
  245. * x3 z3: long form
  246. * x z: short form, destroyed
  247. * xprime zprime: short form, destroyed
  248. * qmqp: short form, preserved
  249. */
  250. static void
  251. fmonty(limb *x2, limb *z2, /* output 2Q */
  252. limb *x3, limb *z3, /* output Q + Q' */
  253. limb *x, limb *z, /* input Q */
  254. limb *xprime, limb *zprime, /* input Q' */
  255. const limb *qmqp /* input Q - Q' */) {
  256. limb origx[5], origxprime[5], zzz[5], xx[5], zz[5], xxprime[5],
  257. zzprime[5], zzzprime[5];
  258. memcpy(origx, x, 5 * sizeof(limb));
  259. fsum(x, z);
  260. fdifference_backwards(z, origx); // does x - z
  261. memcpy(origxprime, xprime, sizeof(limb) * 5);
  262. fsum(xprime, zprime);
  263. fdifference_backwards(zprime, origxprime);
  264. fmul(xxprime, xprime, z);
  265. fmul(zzprime, x, zprime);
  266. memcpy(origxprime, xxprime, sizeof(limb) * 5);
  267. fsum(xxprime, zzprime);
  268. fdifference_backwards(zzprime, origxprime);
  269. fsquare_times(x3, xxprime, 1);
  270. fsquare_times(zzzprime, zzprime, 1);
  271. fmul(z3, zzzprime, qmqp);
  272. fsquare_times(xx, x, 1);
  273. fsquare_times(zz, z, 1);
  274. fmul(x2, xx, zz);
  275. fdifference_backwards(zz, xx); // does zz = xx - zz
  276. fscalar_product(zzz, zz, 121665);
  277. fsum(zzz, xx);
  278. fmul(z2, zz, zzz);
  279. }
  280. // -----------------------------------------------------------------------------
  281. // Maybe swap the contents of two limb arrays (@a and @b), each @len elements
  282. // long. Perform the swap iff @swap is non-zero.
  283. //
  284. // This function performs the swap without leaking any side-channel
  285. // information.
  286. // -----------------------------------------------------------------------------
  287. static void
  288. swap_conditional(limb a[5], limb b[5], limb iswap) {
  289. unsigned i;
  290. const limb swap = -iswap;
  291. for (i = 0; i < 5; ++i) {
  292. const limb x = swap & (a[i] ^ b[i]);
  293. a[i] ^= x;
  294. b[i] ^= x;
  295. }
  296. }
  297. /* Calculates nQ where Q is the x-coordinate of a point on the curve
  298. *
  299. * resultx/resultz: the x coordinate of the resulting curve point (short form)
  300. * n: a little endian, 32-byte number
  301. * q: a point of the curve (short form)
  302. */
  303. static void
  304. cmult(limb *resultx, limb *resultz, const u8 *n, const limb *q) {
  305. limb a[5] = {0}, b[5] = {1}, c[5] = {1}, d[5] = {0};
  306. limb *nqpqx = a, *nqpqz = b, *nqx = c, *nqz = d, *t;
  307. limb e[5] = {0}, f[5] = {1}, g[5] = {0}, h[5] = {1};
  308. limb *nqpqx2 = e, *nqpqz2 = f, *nqx2 = g, *nqz2 = h;
  309. unsigned i, j;
  310. memcpy(nqpqx, q, sizeof(limb) * 5);
  311. for (i = 0; i < 32; ++i) {
  312. u8 byte = n[31 - i];
  313. for (j = 0; j < 8; ++j) {
  314. const limb bit = byte >> 7;
  315. swap_conditional(nqx, nqpqx, bit);
  316. swap_conditional(nqz, nqpqz, bit);
  317. fmonty(nqx2, nqz2,
  318. nqpqx2, nqpqz2,
  319. nqx, nqz,
  320. nqpqx, nqpqz,
  321. q);
  322. swap_conditional(nqx2, nqpqx2, bit);
  323. swap_conditional(nqz2, nqpqz2, bit);
  324. t = nqx;
  325. nqx = nqx2;
  326. nqx2 = t;
  327. t = nqz;
  328. nqz = nqz2;
  329. nqz2 = t;
  330. t = nqpqx;
  331. nqpqx = nqpqx2;
  332. nqpqx2 = t;
  333. t = nqpqz;
  334. nqpqz = nqpqz2;
  335. nqpqz2 = t;
  336. byte <<= 1;
  337. }
  338. }
  339. memcpy(resultx, nqx, sizeof(limb) * 5);
  340. memcpy(resultz, nqz, sizeof(limb) * 5);
  341. }
  342. // -----------------------------------------------------------------------------
  343. // Shamelessly copied from djb's code, tightened a little
  344. // -----------------------------------------------------------------------------
  345. static void
  346. crecip(felem out, const felem z) {
  347. felem a,t0,b,c;
  348. /* 2 */ fsquare_times(a, z, 1); // a = 2
  349. /* 8 */ fsquare_times(t0, a, 2);
  350. /* 9 */ fmul(b, t0, z); // b = 9
  351. /* 11 */ fmul(a, b, a); // a = 11
  352. /* 22 */ fsquare_times(t0, a, 1);
  353. /* 2^5 - 2^0 = 31 */ fmul(b, t0, b);
  354. /* 2^10 - 2^5 */ fsquare_times(t0, b, 5);
  355. /* 2^10 - 2^0 */ fmul(b, t0, b);
  356. /* 2^20 - 2^10 */ fsquare_times(t0, b, 10);
  357. /* 2^20 - 2^0 */ fmul(c, t0, b);
  358. /* 2^40 - 2^20 */ fsquare_times(t0, c, 20);
  359. /* 2^40 - 2^0 */ fmul(t0, t0, c);
  360. /* 2^50 - 2^10 */ fsquare_times(t0, t0, 10);
  361. /* 2^50 - 2^0 */ fmul(b, t0, b);
  362. /* 2^100 - 2^50 */ fsquare_times(t0, b, 50);
  363. /* 2^100 - 2^0 */ fmul(c, t0, b);
  364. /* 2^200 - 2^100 */ fsquare_times(t0, c, 100);
  365. /* 2^200 - 2^0 */ fmul(t0, t0, c);
  366. /* 2^250 - 2^50 */ fsquare_times(t0, t0, 50);
  367. /* 2^250 - 2^0 */ fmul(t0, t0, b);
  368. /* 2^255 - 2^5 */ fsquare_times(t0, t0, 5);
  369. /* 2^255 - 21 */ fmul(out, t0, a);
  370. }
  371. int curve25519_donna(u8 *, const u8 *, const u8 *);
  372. int
  373. curve25519_donna(u8 *mypublic, const u8 *secret, const u8 *basepoint) {
  374. limb bp[5], x[5], z[5], zmone[5];
  375. uint8_t e[32];
  376. int i;
  377. for (i = 0;i < 32;++i) e[i] = secret[i];
  378. e[0] &= 248;
  379. e[31] &= 127;
  380. e[31] |= 64;
  381. fexpand(bp, basepoint);
  382. cmult(x, z, e, bp);
  383. crecip(zmone, z);
  384. fmul(z, x, zmone);
  385. fcontract(mypublic, z);
  386. return 0;
  387. }