mont.c 40 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411
  1. /* ===================================================================
  2. *
  3. * Copyright (c) 2018, Helder Eijs <helderijs@gmail.com>
  4. * All rights reserved.
  5. *
  6. * Redistribution and use in source and binary forms, with or without
  7. * modification, are permitted provided that the following conditions
  8. * are met:
  9. *
  10. * 1. Redistributions of source code must retain the above copyright
  11. * notice, this list of conditions and the following disclaimer.
  12. * 2. Redistributions in binary form must reproduce the above copyright
  13. * notice, this list of conditions and the following disclaimer in
  14. * the documentation and/or other materials provided with the
  15. * distribution.
  16. *
  17. * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  18. * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  19. * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  20. * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  21. * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  22. * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
  23. * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  24. * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  25. * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  26. * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
  27. * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28. * POSSIBILITY OF SUCH DAMAGE.
  29. * ===================================================================
  30. */
  31. #include <assert.h>
  32. #include "common.h"
  33. #include "endianess.h"
  34. #include "multiply.h"
  35. #include "modexp_utils.h"
  36. #include "mont.h"
  37. #if SYS_BITS == 32
  38. #include "multiply_32.c"
  39. #else
  40. #if SYS_BITS == 64
  41. #include "multiply_64.c"
  42. #else
  43. #error You must define the macro SYS_BITS
  44. #endif
  45. #endif
  46. #if defined(USE_SSE2)
  47. #if defined(HAVE_INTRIN_H)
  48. #include <intrin.h>
  49. #elif defined(HAVE_X86INTRIN_H)
  50. #include <x86intrin.h>
  51. #elif defined(HAVE_EMMINTRIN_H)
  52. #include <xmmintrin.h>
  53. #include <emmintrin.h>
  54. #endif
  55. #endif
  56. static inline unsigned is_odd(uint64_t x)
  57. {
  58. return 1 == (x & 1);
  59. }
  60. static inline unsigned is_even(uint64_t x)
  61. {
  62. return !is_odd(x);
  63. }
  64. /**
  65. * Compute the inverse modulo 2⁶⁴ of a 64-bit odd integer.
  66. *
  67. * See https://crypto.stackexchange.com/questions/47493/how-to-determine-the-multiplicative-inverse-modulo-64-or-other-power-of-two
  68. */
  69. STATIC uint64_t inverse64(uint64_t a)
  70. {
  71. uint64_t x;
  72. assert(1 & a);
  73. x = ((a << 1 ^ a) & 4) << 1 ^ a;
  74. x += x - a*x*x;
  75. x += x - a*x*x;
  76. x += x - a*x*x;
  77. x += x - a*x*x;
  78. assert((x*a & 0xFFFFFFFFFFFFFFFFULL) == 1);
  79. return x;
  80. }
  81. /**
  82. * Check if a multi-word integer x is greater than or equal to y.
  83. *
  84. * @param x The first term
  85. * @param y The second term
  86. * @param nw The number of words that make up x and y
  87. * @return 1 if x>=y, 0 if x<y
  88. */
  89. STATIC int ge(const uint64_t *x, const uint64_t *y, size_t nw)
  90. {
  91. unsigned mask = (unsigned)-1;
  92. unsigned result = 0;
  93. size_t i, j;
  94. i = nw - 1;
  95. for (j=0; j<nw; j++, i--) {
  96. unsigned greater, lower;
  97. greater = x[i] > y[i];
  98. lower = x[i] < y[i];
  99. result |= mask & (greater | (lower << 1));
  100. mask &= (greater ^ lower) - 1;
  101. }
  102. return result<2;
  103. }
  104. /*
  105. * Subtract a multi-word integer b from a.
  106. *
  107. * @param out The location where the multi-word result is stored
  108. * @param a Number to subtract from
  109. * @param b Number to subtract
  110. * @param nw The number of words of both a and b
  111. * @result 0 if there is no borrow, 1 otherwise
  112. */
  113. STATIC unsigned sub(uint64_t *out, const uint64_t *a, const uint64_t *b, size_t nw)
  114. {
  115. size_t i;
  116. unsigned borrow1 , borrow2;
  117. borrow2 = 0;
  118. for (i=0; i<nw; i++) {
  119. borrow1 = b[i] > a[i];
  120. out[i] = a[i] - b[i];
  121. borrow1 |= borrow2 > out[i];
  122. out[i] -= borrow2;
  123. borrow2 = borrow1;
  124. }
  125. return borrow2;
  126. }
  127. /*
  128. * Compute R² mod N, where R is the smallest power of 2⁶⁴ larger than N.
  129. *
  130. * @param r2_mod_n The location where the result is stored at
  131. * @param n The modulus N
  132. * @param nw The number of 64-bit words of both r2_mod_n and n
  133. */
  134. STATIC void rsquare(uint64_t *r2_mod_n, uint64_t *n, size_t nw)
  135. {
  136. size_t i;
  137. size_t R_bits;
  138. memset(r2_mod_n, 0, sizeof(uint64_t)*nw);
  139. /**
  140. * Start with R2=1, double 2*bitlen(R) times,
  141. * and reduce it as soon as it exceeds n
  142. */
  143. r2_mod_n[0] = 1;
  144. R_bits = nw * sizeof(uint64_t) * 8;
  145. for (i=0; i<R_bits*2; i++) {
  146. unsigned overflow;
  147. size_t j;
  148. /** Double, by shifting left by one bit **/
  149. overflow = (unsigned)(r2_mod_n[nw-1] >> 63);
  150. for (j=nw-1; j>0; j--) {
  151. r2_mod_n[j] = (r2_mod_n[j] << 1) + (r2_mod_n[j-1] >> 63);
  152. }
  153. /** Fill-in with zeroes **/
  154. r2_mod_n[0] <<= 1;
  155. /** Subtract n if the result exceeds it **/
  156. while (overflow || ge(r2_mod_n, n, nw)) {
  157. sub(r2_mod_n, r2_mod_n, n, nw);
  158. overflow = 0;
  159. }
  160. }
  161. }
  162. /*
  163. * Multiply a multi-word integer a by a 64-bit scalar k and
  164. * then add the result to the multi-word integer t.
  165. *
  166. * @param t The multi-word integer accumulator
  167. * @param tw The number of words of t
  168. * @param a The multi-word integer to multiply with the scalar
  169. * @param aw The number of words of a
  170. * @param k The 64-bit scalar multiplier
  171. */
  172. STATIC void addmul(uint64_t *t, size_t tw, const uint64_t *a, size_t aw, uint64_t k)
  173. {
  174. size_t i;
  175. uint64_t carry;
  176. carry = 0;
  177. for (i=0; i<aw; i++) {
  178. uint64_t prod_lo, prod_hi;
  179. DP_MULT(a[i], k, prod_lo, prod_hi);
  180. prod_lo += carry;
  181. prod_hi += prod_lo < carry;
  182. t[i] += prod_lo;
  183. prod_hi += t[i] < prod_lo;
  184. carry = prod_hi;
  185. }
  186. for (; carry; i++) {
  187. t[i] += carry;
  188. carry = t[i] < carry;
  189. }
  190. assert(i <= tw);
  191. }
  192. /**
  193. * Multiply two multi-word integers.
  194. *
  195. * @param t The location where the result is stored. It is twice as big as
  196. * either a (or b). It is an array of 2*nw words).
  197. * @param scratchpad Temporary area. It is an array of 3*nw words.
  198. * @param a The first term, array of nw words.
  199. * @param b The second term, array of nw words.
  200. * @param nw The number of words of both a and b.
  201. *
  202. */
  203. STATIC void product(uint64_t *t, uint64_t *scratchpad, const uint64_t *a, const uint64_t *b, size_t nw)
  204. {
  205. size_t i;
  206. memset(t, 0, 2*sizeof(uint64_t)*nw);
  207. for (i=0; i<(nw ^ (nw & 1)); i+=2) {
  208. addmul128(&t[i], scratchpad, a, b[i], b[i+1], 2*nw-i, nw);
  209. }
  210. if (is_odd(nw)) {
  211. addmul(&t[nw-1], nw+2, a, nw, b[nw-1]);
  212. }
  213. }
  214. /*
  215. * Select a number out of two, in constant time.
  216. *
  217. * @param out The location where the multi-word result is stored
  218. * @param a The first choice, selected if cond is true (non-zero)
  219. * @param b The second choice, selected if cond is false (zero)
  220. * @param cond The flag that drives the selection
  221. * @param words The number of words of a, b, and out
  222. * @return 0 for success, the appropriate code otherwise.
  223. */
  224. STATIC FUNC_SSE2 int mont_select(uint64_t *out, const uint64_t *a, const uint64_t *b, unsigned cond, size_t words)
  225. {
  226. uint64_t mask;
  227. #if defined(USE_SSE2)
  228. unsigned pairs, i;
  229. __m128i r0, r1, r2, r3, r4, r5;
  230. pairs = (unsigned)words / 2;
  231. mask = (uint64_t)((cond != 0) - 1); /* 0 for a, 1s for b */
  232. #if SYSBITS == 64
  233. r0 = _mm_set1_epi64x(mask);
  234. #else
  235. r0 = _mm_loadl_epi64((__m128i*)&mask);
  236. r0 = _mm_unpacklo_epi64(r0, r0);
  237. #endif
  238. for (i=0; i<pairs; i++, a+=2, b+=2, out+=2) {
  239. r1 = _mm_loadu_si128((__m128i const*)b);
  240. r2 = _mm_loadu_si128((__m128i const*)a);
  241. r3 = _mm_and_si128(r0, r1);
  242. r4 = _mm_andnot_si128(r0, r2);
  243. r5 = _mm_or_si128(r3, r4);
  244. _mm_storeu_si128((__m128i*)out, r5);
  245. }
  246. if (words & 1) {
  247. *out = (*b & mask) ^ (*a & ~mask);
  248. }
  249. #else
  250. unsigned i;
  251. mask = (uint64_t)((cond != 0) - 1);
  252. for (i=0; i<words; i++) {
  253. *out++ = (*b++ & mask) ^ (*a++ & ~mask);
  254. }
  255. #endif
  256. return 0;
  257. }
  258. /*
  259. * Add two multi-word numbers with modulo arithmetic.
  260. *
  261. * @param out The locaton where the multi-word result (nw words) is stored
  262. * @param a The first term (nw words)
  263. * @param b The second term (nw words)
  264. * @param modulus The modulus (nw words)
  265. * @param tmp1 A temporary area (nw words)
  266. * @param tmp2 A temporary area (nw words)
  267. * @param nw The number of 64-bit words in all parameters
  268. */
  269. void add_mod(uint64_t* out, const uint64_t* a, const uint64_t* b, const uint64_t *modulus, uint64_t *tmp1, uint64_t *tmp2, size_t nw)
  270. {
  271. unsigned i;
  272. unsigned carry, borrow1, borrow2;
  273. /*
  274. * Compute sum in tmp1[], and subtract modulus[]
  275. * from tmp1[] into tmp2[].
  276. */
  277. borrow2 = 0;
  278. for (i=0, carry=0; i<nw; i++) {
  279. tmp1[i] = a[i] + carry;
  280. carry = tmp1[i] < carry;
  281. tmp1[i] += b[i];
  282. carry += tmp1[i] < b[i];
  283. borrow1 = modulus[i] > tmp1[i];
  284. tmp2[i] = tmp1[i] - modulus[i];
  285. borrow1 |= borrow2 > tmp2[i];
  286. tmp2[i] -= borrow2;
  287. borrow2 = borrow1;
  288. }
  289. /*
  290. * If there is no borrow or if there is carry,
  291. * tmp1[] is larger than modulus, so we must return tmp2[].
  292. */
  293. mont_select(out, tmp2, tmp1, carry | (borrow2 ^ 1), nw);
  294. }
  295. /*
  296. * Montgomery modular multiplication, that is a*b*R mod N.
  297. *
  298. * @param out The location where the result is stored
  299. * @param a The first term (already in Montgomery form, a*R mod N)
  300. * @param b The second term (already in Montgomery form, b*R mod N)
  301. * @param n The modulus (in normal form), such that R>N
  302. * @param m0 Least-significant word of the opposite of the inverse of n modulo R, that is, -n[0]⁻¹ mod R
  303. * @param t Temporary, internal result; it must have been created with mont_number(&p,SCRATCHPAD_NR,ctx).
  304. * @param nw Number of words making up the 3 integers: out, a, and b.
  305. * It also defines R as 2^(64*nw).
  306. *
  307. * Useful read: https://alicebob.cryptoland.net/understanding-the-montgomery-reduction-algorithm/
  308. */
  309. #if SCRATCHPAD_NR < 7
  310. #error Scratchpad is too small
  311. #endif
  312. STATIC void mont_mult_generic(uint64_t *out, const uint64_t *a, const uint64_t *b, const uint64_t *n, uint64_t m0, uint64_t *tmp, size_t nw)
  313. {
  314. size_t i;
  315. uint64_t *t, *scratchpad, *t2;
  316. unsigned cond;
  317. /*
  318. * tmp is an array of SCRATCHPAD*nw words
  319. * We carve out 3 values in it:
  320. * - 3*nw words, the value a*b + m*n (we only use 2*nw+1 words)
  321. * - 3*nw words, temporary area for computing the product
  322. * - nw words, the reduced value with a final subtraction by n
  323. */
  324. t = tmp;
  325. scratchpad = tmp + 3*nw;
  326. t2 = scratchpad + 3*nw;
  327. if (a == b) {
  328. square(t, scratchpad, a, nw);
  329. } else {
  330. product(t, scratchpad, a, b, nw);
  331. }
  332. t[2*nw] = 0; /** MSW **/
  333. /** Clear lower words (two at a time) **/
  334. for (i=0; i<(nw ^ (nw & 1)); i+=2) {
  335. uint64_t k0, k1, ti1, prod_lo, prod_hi;
  336. /** Multiplier for n that will make t[i+0] go 0 **/
  337. k0 = t[i] * m0;
  338. /** Simulate Muladd for digit 0 **/
  339. DP_MULT(k0, n[0], prod_lo, prod_hi);
  340. prod_lo += t[i];
  341. prod_hi += prod_lo < t[i];
  342. /** Expected digit 1 **/
  343. ti1 = t[i+1] + n[1]*k0 + prod_hi;
  344. /** Multiplier for n that will make t[i+1] go 0 **/
  345. k1 = ti1 * m0;
  346. addmul128(&t[i], scratchpad, n, k0, k1, 2*nw+1-i, nw);
  347. }
  348. /** One left for odd number of words **/
  349. if (is_odd(nw)) {
  350. addmul(&t[nw-1], nw+2, n, nw, t[nw-1]*m0);
  351. }
  352. assert(t[2*nw] <= 1); /** MSW **/
  353. /** t[0..nw-1] == 0 **/
  354. /** Divide by R and possibly subtract n **/
  355. sub(t2, &t[nw], n, nw);
  356. cond = (unsigned)(t[2*nw] | (uint64_t)ge(&t[nw], n, nw));
  357. mont_select(out, t2, &t[nw], cond, (unsigned)nw);
  358. }
  359. STATIC void mont_mult_p256(uint64_t *out, const uint64_t *a, const uint64_t *b, const uint64_t *n, uint64_t m0, uint64_t *tmp, size_t nw)
  360. {
  361. unsigned i;
  362. uint64_t *t, *scratchpad, *t2;
  363. unsigned cond;
  364. #define WORDS_64 4U
  365. #define PREDIV_WORDS_64 (2*WORDS_64+1) /** Size of the number to divide by R **/
  366. #define WORDS_32 (WORDS_64*2)
  367. #define PREDIV_WORDS_32 (2*PREDIV_WORDS_64)
  368. #if SYS_BITS == 32
  369. uint32_t t32[18];
  370. #endif
  371. assert(nw == 4);
  372. assert(m0 == 1);
  373. t = tmp;
  374. scratchpad = tmp + 3*nw;
  375. t2 = scratchpad + 3*nw;
  376. if (a == b) {
  377. square(t, scratchpad, a, WORDS_64);
  378. } else {
  379. product(t, scratchpad, a, b, WORDS_64);
  380. }
  381. t[PREDIV_WORDS_64-1] = 0; /** MSW **/
  382. #if SYS_BITS == 32
  383. for (i=0; i<PREDIV_WORDS_64; i++) {
  384. t32[2*i] = (uint32_t)t[i];
  385. t32[2*i+1] = (uint32_t)(t[i] >> 32);
  386. }
  387. for (i=0; i<WORDS_32; i++) {
  388. uint32_t k, carry;
  389. uint64_t prod, k2;
  390. unsigned j;
  391. k = t32[i];
  392. k2 = ((uint64_t)k<<32) - k;
  393. /* p[0] = 2³²-1 */
  394. prod = k2 + t32[i+0];
  395. t32[i+0] = (uint32_t)prod;
  396. carry = (uint32_t)(prod >> 32);
  397. /* p[1] = 2³²-1 */
  398. prod = k2 + t32[i+1] + carry;
  399. t32[i+1] = (uint32_t)prod;
  400. carry = (uint32_t)(prod >> 32);
  401. /* p[2] = 2³²-1 */
  402. prod = k2 + t32[i+2] + carry;
  403. t32[i+2] = (uint32_t)prod;
  404. carry = (uint32_t)(prod >> 32);
  405. /* p[3] = 0 */
  406. t32[i+3] += carry;
  407. carry = t32[i+3] < carry;
  408. /* p[4] = 0 */
  409. t32[i+4] += carry;
  410. carry = t32[i+4] < carry;
  411. /* p[5] = 0 */
  412. t32[i+5] += carry;
  413. carry = t32[i+5] < carry;
  414. /* p[6] = 1 */
  415. t32[i+6] += carry;
  416. carry = t32[i+6] < carry;
  417. t32[i+6] += k;
  418. carry |= t32[i+6] < k;
  419. /* p[7] = 2³²-1 */
  420. prod = k2 + t32[i+7] + carry;
  421. t32[i+7] = (uint32_t)prod;
  422. carry = (uint32_t)(prod >> 32);
  423. for (j=WORDS_32; carry; j++) {
  424. t32[i+j] += carry;
  425. carry = t32[i+j] < carry;
  426. }
  427. }
  428. for (i=0; i<PREDIV_WORDS_64; i++) {
  429. t[i] = ((uint64_t)t32[2*i+1]<<32) + t32[2*i];
  430. }
  431. #elif SYS_BITS == 64
  432. for (i=0; i<WORDS_64; i++) {
  433. unsigned j;
  434. uint64_t carry, k;
  435. uint64_t prod_lo, prod_hi;
  436. k = t[i];
  437. /* n[0] = 2⁶⁴ - 1 */
  438. prod_lo = -k;
  439. prod_hi = k - (k!=0);
  440. t[i+0] += prod_lo;
  441. prod_hi += t[i+0] < prod_lo;
  442. carry = prod_hi;
  443. /* n[1] = 2³² - 1 */
  444. DP_MULT(n[1], k, prod_lo, prod_hi);
  445. prod_lo += carry;
  446. prod_hi += prod_lo < carry;
  447. t[i+1] += prod_lo;
  448. prod_hi += t[i+1] < prod_lo;
  449. carry = prod_hi;
  450. /* n[2] = 0 */
  451. t[i+2] += carry;
  452. carry = t[i+2] < carry;
  453. /* n[3] = 2⁶⁴ - 2³² + 1 */
  454. DP_MULT(n[3], k, prod_lo, prod_hi);
  455. prod_lo += carry;
  456. prod_hi += prod_lo < carry;
  457. t[i+3] += prod_lo;
  458. prod_hi += t[i+3] < prod_lo;
  459. carry = prod_hi;
  460. for (j=WORDS_64; carry; j++) {
  461. t[i+j] += carry;
  462. carry = t[i+j] < carry;
  463. }
  464. }
  465. #else
  466. #error You must define the SYS_BITS macro
  467. #endif
  468. assert(t[PREDIV_WORDS_64-1] <= 1); /** MSW **/
  469. /** t[0..nw-1] == 0 **/
  470. /** Divide by R and possibly subtract n **/
  471. sub(t2, &t[nw], n, WORDS_64);
  472. cond = (unsigned)(t[PREDIV_WORDS_64-1] | (uint64_t)ge(&t[WORDS_64], n, WORDS_64));
  473. mont_select(out, t2, &t[WORDS_64], cond, WORDS_64);
  474. #undef WORDS_64
  475. #undef PREDIV_WORDS_64
  476. #undef WORDS_32
  477. #undef PREDIV_WORDS_32
  478. }
  479. STATIC void mont_mult_p384(uint64_t *out, const uint64_t *a, const uint64_t *b, const uint64_t *n, uint64_t m0, uint64_t *tmp, size_t nw)
  480. {
  481. size_t i;
  482. uint64_t *t, *scratchpad, *t2;
  483. unsigned cond;
  484. #define WORDS_64 6U
  485. #define PREDIV_WORDS_64 (2*WORDS_64+1) /** Size of the number to divide by R **/
  486. #define WORDS_32 (WORDS_64*2)
  487. #define PREDIV_WORDS_32 (2*PREDIV_WORDS_64)
  488. #if SYS_BITS == 32
  489. uint32_t t32[PREDIV_WORDS_32];
  490. #endif
  491. assert(nw == WORDS_64);
  492. assert(m0 == 0x0000000100000001ULL);
  493. t = tmp;
  494. scratchpad = tmp + 3*nw;
  495. t2 = scratchpad + 3*nw;
  496. if (a == b) {
  497. square(t, scratchpad, a, WORDS_64);
  498. } else {
  499. product(t, scratchpad, a, b, WORDS_64);
  500. }
  501. t[PREDIV_WORDS_64-1] = 0; /** MSW **/
  502. #if SYS_BITS == 32
  503. for (i=0; i<PREDIV_WORDS_64; i++) {
  504. t32[2*i] = (uint32_t)t[i];
  505. t32[2*i+1] = (uint32_t)(t[i] >> 32);
  506. }
  507. for (i=0; i<WORDS_32; i++) {
  508. uint32_t k, carry;
  509. uint64_t prod, k2, k3;
  510. unsigned j;
  511. k = t32[i];
  512. k2 = ((uint64_t)k<<32) - k;
  513. k3 = k2 - k;
  514. /* n32[0] = 2³² - 1 */
  515. prod = k2 + t32[i+0];
  516. t32[i+0] = (uint32_t)prod;
  517. carry = (uint32_t)(prod >> 32);
  518. /* n32[1] = 0 */
  519. prod = (uint64_t)t32[i+1] + carry;
  520. t32[i+1] = (uint32_t)prod;
  521. carry = (uint32_t)(prod >> 32);
  522. /* n32[2] = 0 */
  523. prod = (uint64_t)t32[i+2] + carry;
  524. t32[i+2] = (uint32_t)prod;
  525. carry = (uint32_t)(prod >> 32);
  526. /* n32[3] = 2³² - 1 */
  527. prod = k2 + t32[i+3] + carry;
  528. t32[i+3] = (uint32_t)prod;
  529. carry = (uint32_t)(prod >> 32);
  530. /* n32[4] = 2³² - 2 */
  531. prod = k3 + t32[i+4] + carry;
  532. t32[i+4] = (uint32_t)prod;
  533. carry = (uint32_t)(prod >> 32);
  534. /* n32[5] = 2³² - 1 */
  535. prod = k2 + t32[i+5] + carry;
  536. t32[i+5] = (uint32_t)prod;
  537. carry = (uint32_t)(prod >> 32);
  538. /* n32[6] = 2³² - 1 */
  539. prod = k2 + t32[i+6] + carry;
  540. t32[i+6] = (uint32_t)prod;
  541. carry = (uint32_t)(prod >> 32);
  542. /* n32[7] = 2³² - 1 */
  543. prod = k2 + t32[i+7] + carry;
  544. t32[i+7] = (uint32_t)prod;
  545. carry = (uint32_t)(prod >> 32);
  546. /* n32[8] = 2³² - 1 */
  547. prod = k2 + t32[i+8] + carry;
  548. t32[i+8] = (uint32_t)prod;
  549. carry = (uint32_t)(prod >> 32);
  550. /* n32[9] = 2³² - 1 */
  551. prod = k2 + t32[i+9] + carry;
  552. t32[i+9] = (uint32_t)prod;
  553. carry = (uint32_t)(prod >> 32);
  554. /* n32[10] = 2³² - 1 */
  555. prod = k2 + t32[i+10] + carry;
  556. t32[i+10] = (uint32_t)prod;
  557. carry = (uint32_t)(prod >> 32);
  558. /* n32[11] = 2³² - 1 */
  559. prod = k2 + t32[i+11] + carry;
  560. t32[i+11] = (uint32_t)prod;
  561. carry = (uint32_t)(prod >> 32);
  562. for (j=WORDS_32; carry; j++) {
  563. t32[i+j] += carry;
  564. carry = t32[i+j] < carry;
  565. }
  566. }
  567. for (i=0; i<PREDIV_WORDS_64; i++) {
  568. t[i] = ((uint64_t)t32[2*i+1]<<32) + t32[2*i];
  569. }
  570. #elif SYS_BITS == 64
  571. for (i=0; i<WORDS_64; i++) {
  572. unsigned j;
  573. uint64_t carry;
  574. uint64_t k, k2_lo, k2_hi;
  575. uint64_t prod_lo, prod_hi;
  576. k = t[i] + (t[i] << 32);
  577. k2_lo = -k;
  578. k2_hi = k - (k!=0);
  579. /* n[0] = 2³² - 1 */
  580. DP_MULT(n[0], k, prod_lo, prod_hi);
  581. t[i+0] += prod_lo;
  582. prod_hi += t[i+0] < prod_lo;
  583. carry = prod_hi;
  584. /* n[1] = 2⁶⁴ - 2³² */
  585. DP_MULT(n[1], k, prod_lo, prod_hi);
  586. prod_lo += carry;
  587. prod_hi += prod_lo < carry;
  588. t[i+1] += prod_lo;
  589. prod_hi += t[i+1] < prod_lo;
  590. carry = prod_hi;
  591. /* n[2] = 2⁶⁴ - 2 */
  592. DP_MULT(n[2], k, prod_lo, prod_hi);
  593. prod_lo += carry;
  594. prod_hi += prod_lo < carry;
  595. t[i+2] += prod_lo;
  596. prod_hi += t[i+2] < prod_lo;
  597. carry = prod_hi;
  598. /* n[3] = 2⁶⁴ - 1 */
  599. prod_lo = k2_lo;
  600. prod_hi = k2_hi;
  601. prod_lo += carry;
  602. prod_hi += prod_lo < carry;
  603. t[i+3] += prod_lo;
  604. prod_hi += t[i+3] < prod_lo;
  605. carry = prod_hi;
  606. /* n[4] = 2⁶⁴ - 1 */
  607. prod_lo = k2_lo;
  608. prod_hi = k2_hi;
  609. prod_lo += carry;
  610. prod_hi += prod_lo < carry;
  611. t[i+4] += prod_lo;
  612. prod_hi += t[i+4] < prod_lo;
  613. carry = prod_hi;
  614. /* n[5] = 2⁶⁴ - 1 */
  615. prod_lo = k2_lo;
  616. prod_hi = k2_hi;
  617. prod_lo += carry;
  618. prod_hi += prod_lo < carry;
  619. t[i+5] += prod_lo;
  620. prod_hi += t[i+5] < prod_lo;
  621. carry = prod_hi;
  622. for (j=WORDS_64; carry; j++) {
  623. t[i+j] += carry;
  624. carry = t[i+j] < carry;
  625. }
  626. }
  627. #else
  628. #error You must define the SYS_BITS macro
  629. #endif
  630. assert(t[PREDIV_WORDS_64-1] <= 1); /** MSW **/
  631. /** Words t[0..WORDS_64-1] have all been set to zero **/
  632. /** Divide by R and possibly subtract n **/
  633. sub(t2, &t[WORDS_64], n, WORDS_64);
  634. cond = (unsigned)(t[PREDIV_WORDS_64-1] | (uint64_t)ge(&t[WORDS_64], n, WORDS_64));
  635. mont_select(out, t2, &t[WORDS_64], cond, WORDS_64);
  636. #undef WORDS_64
  637. #undef PREDIV_WORDS_64
  638. #undef WORDS_32
  639. #undef PREDIV_WORDS_32
  640. }
  641. STATIC void mont_mult_p521(uint64_t *out, const uint64_t *a, const uint64_t *b, const uint64_t *n, uint64_t m0, uint64_t *tmp, size_t nw)
  642. {
  643. uint64_t *t, *scratchpad, *s, *tmp1, *tmp2;
  644. assert(nw == 9);
  645. assert(m0 == 1);
  646. /*
  647. * A number in the form:
  648. * x*2⁵²¹ + y
  649. * is congruent modulo 2⁵²¹-1 to:
  650. * x + y
  651. */
  652. /*
  653. * tmp is an array of SCRATCHPAD*nw words
  654. * We carve out 3 values in it:
  655. * - 2*nw words, the value a*b
  656. * - 3*nw words, temporary area for computing the product
  657. * - nw words, the second term of the addition
  658. */
  659. t = tmp;
  660. scratchpad = t + 2*nw;
  661. s = scratchpad + 3*nw;
  662. tmp1 = scratchpad;
  663. tmp2 = scratchpad + nw;
  664. if (a == b) {
  665. square(t, scratchpad, a, 9);
  666. } else {
  667. product(t, scratchpad, a, b, 9);
  668. }
  669. /* t is a 1042-bit number, occupying 17 words (of the total 18); the MSW (t[16]) only has 18 bits */
  670. s[0] = (t[8] >> 9) | (t[9] << 55); t[8] &= 0x1FF;
  671. s[1] = (t[9] >> 9) | (t[10] << 55);
  672. s[2] = (t[10] >> 9) | (t[11] << 55);
  673. s[3] = (t[11] >> 9) | (t[12] << 55);
  674. s[4] = (t[12] >> 9) | (t[13] << 55);
  675. s[5] = (t[13] >> 9) | (t[14] << 55);
  676. s[6] = (t[14] >> 9) | (t[15] << 55);
  677. s[7] = (t[15] >> 9) | (t[16] << 55);
  678. s[8] = t[16] >> 9;
  679. add_mod(out, t, s, n, tmp1, tmp2, nw);
  680. }
  681. /* ---- PUBLIC FUNCTIONS ---- */
  682. void mont_context_free(MontContext *ctx)
  683. {
  684. if (NULL == ctx)
  685. return;
  686. free(ctx->one);
  687. free(ctx->r2_mod_n);
  688. free(ctx->r_mod_n);
  689. free(ctx->modulus);
  690. free(ctx->modulus_min_2);
  691. free(ctx);
  692. }
  693. /*
  694. * Return how many bytes a big endian multi-word number takes in memory.
  695. */
  696. size_t mont_bytes(const MontContext *ctx)
  697. {
  698. if (NULL == ctx)
  699. return 0;
  700. return ctx->bytes;
  701. }
  702. /*
  703. * Allocate memory for an array of numbers in Montgomery form
  704. * and initialize it to 0.
  705. *
  706. * @param out The location where the address of the newly allocated
  707. * array will be placed in.
  708. * The caller is responsible for deallocating the memory
  709. * using free().
  710. * @param count How many numbers the array contains.
  711. * @param ctx The Montgomery context.
  712. * @return 0 if successful, the relevant error code otherwise.
  713. *
  714. */
  715. int mont_number(uint64_t **out, unsigned count, const MontContext *ctx)
  716. {
  717. if (NULL == out || NULL == ctx)
  718. return ERR_NULL;
  719. *out = (uint64_t*)calloc(count * ctx->words, sizeof(uint64_t));
  720. if (NULL == *out)
  721. return ERR_MEMORY;
  722. return 0;
  723. }
  724. int mont_random_number(uint64_t **out, unsigned count, uint64_t seed, const MontContext *ctx)
  725. {
  726. int res;
  727. unsigned i;
  728. uint64_t *number;
  729. res = mont_number(out, count, ctx);
  730. if (res)
  731. return res;
  732. number = *out;
  733. expand_seed(seed, (uint8_t*)number, count * ctx->bytes);
  734. for (i=0; i<count; i++, number += ctx->words) {
  735. number[ctx->words-1] = 0;
  736. }
  737. return 0;
  738. }
  739. /*
  740. * Transform a big endian-encoded number into Montgomery form, by performing memory allocation.
  741. *
  742. * @param out The location where the pointer to the newly allocated memory will be put in.
  743. * The memory will contain the number encoded in Montgomery form.
  744. * The caller is responsible for deallocating the memory.
  745. * @param ctx Montgomery context, as created by mont_context_init().
  746. * @param number The big endian-encoded number to transform, strictly smaller than the modulus.
  747. * @param len The length of the big-endian number in bytes (this may be
  748. * smaller than the output of mont_bytes(ctx)).
  749. * @return 0 in case of success, the relevant error code otherwise.
  750. */
  751. int mont_from_bytes(uint64_t **out, const uint8_t *number, size_t len, const MontContext *ctx)
  752. {
  753. uint64_t *encoded = NULL;
  754. uint64_t *tmp1 = NULL;
  755. uint64_t *scratchpad = NULL;
  756. int res = 0;
  757. if (NULL == out || NULL == ctx || NULL == number)
  758. return ERR_NULL;
  759. *out = NULL;
  760. /** Removing leading zeroes but avoid a zero-length string **/
  761. if (0 == len)
  762. return ERR_NOT_ENOUGH_DATA;
  763. while (len>1 && *number==0) {
  764. len--;
  765. number++;
  766. }
  767. if (ctx->bytes < len)
  768. return ERR_VALUE;
  769. /** The caller will deallocate this memory **/
  770. *out = encoded = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  771. if (NULL == encoded)
  772. return ERR_MEMORY;
  773. /** Input number, loaded in words **/
  774. tmp1 = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  775. if (NULL == tmp1) {
  776. res = ERR_MEMORY;
  777. goto cleanup;
  778. }
  779. bytes_to_words(tmp1, ctx->words, number, len);
  780. /** Make sure number<modulus **/
  781. if (ge(tmp1, ctx->modulus, ctx->words)) {
  782. res = ERR_VALUE;
  783. goto cleanup;
  784. }
  785. /** Scratchpad **/
  786. scratchpad = (uint64_t*)calloc(SCRATCHPAD_NR, ctx->words*sizeof(uint64_t));
  787. if (NULL == scratchpad) {
  788. res = ERR_MEMORY;
  789. goto cleanup;
  790. }
  791. if (ctx->modulus_type != ModulusP521)
  792. mont_mult_generic(encoded, tmp1, ctx->r2_mod_n, ctx->modulus, ctx->m0, scratchpad, ctx->words);
  793. else
  794. mont_copy(encoded, tmp1, ctx);
  795. res = 0;
  796. cleanup:
  797. free(scratchpad);
  798. free(tmp1);
  799. if (res != 0) {
  800. free(encoded);
  801. *out = NULL;
  802. }
  803. return res;
  804. }
  805. /*
  806. * Transform a number from Montgomery representation to big endian-encoding.
  807. *
  808. * @param number The location where the number will be put in, encoded
  809. * in big-endian form and with zero padding on the left.
  810. * @param len Space allocate at number, at least ctx->modulus_len bytes.
  811. * @param ctx The address of the Montgomery context.
  812. * @param mont_number The number in Montgomery form to transform.
  813. * @return 0 if successful, the relevant error code otherwise.
  814. */
  815. int mont_to_bytes(uint8_t *number, size_t len, const uint64_t* mont_number, const MontContext *ctx)
  816. {
  817. uint64_t *tmp1 = NULL;
  818. uint64_t *scratchpad = NULL;
  819. int res;
  820. if (NULL == number || NULL == ctx || NULL == mont_number)
  821. return ERR_NULL;
  822. if (len < ctx->modulus_len)
  823. return ERR_NOT_ENOUGH_DATA;
  824. /** Number in normal form, but still in words **/
  825. tmp1 = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  826. if (NULL == tmp1)
  827. return ERR_MEMORY;
  828. /** Scratchpad **/
  829. scratchpad = (uint64_t*)calloc(SCRATCHPAD_NR, ctx->words*sizeof(uint64_t));
  830. if (NULL == scratchpad) {
  831. free(tmp1);
  832. return ERR_MEMORY;
  833. }
  834. if (ctx->modulus_type != ModulusP521)
  835. mont_mult_generic(tmp1, mont_number, ctx->one, ctx->modulus, ctx->m0, scratchpad, ctx->words);
  836. else
  837. mont_copy(tmp1, mont_number, ctx);
  838. res = words_to_bytes(number, len, tmp1, ctx->words);
  839. free(scratchpad);
  840. free(tmp1);
  841. return res;
  842. }
  843. /*
  844. * Add two numbers in Montgomery representation.
  845. *
  846. * @param out The location where the result will be stored; it must have been created with mont_number(&p,1,ctx).
  847. * @param a The first term.
  848. * @param b The second term.
  849. * @param tmp Temporary, internal result; it must have been created with mont_number(&p,SCRATCHPAD_NR,ctx).
  850. * @param ctx The Montgomery context.
  851. * @return 0 for success, the relevant error code otherwise.
  852. */
  853. int mont_add(uint64_t* out, const uint64_t* a, const uint64_t* b, uint64_t *tmp, const MontContext *ctx)
  854. {
  855. if (NULL == out || NULL == a || NULL == b || NULL == tmp || NULL == ctx)
  856. return ERR_NULL;
  857. add_mod(out, a, b, ctx->modulus, tmp, tmp + ctx->words, ctx->words);
  858. return 0;
  859. }
  860. /*
  861. * Multiply two numbers in Montgomery representation.
  862. *
  863. * @param out The location where the result will be stored at; it must have been created with mont_number(&p,1,ctx)
  864. * @param a The first term.
  865. * @param b The second term.
  866. * @param tmp Temporary, internal result; it must have been created with mont_number(&p,SCRATCHPAD_NR,ctx).
  867. * @param ctx The Montgomery context.
  868. * @return 0 for success, the relevant error code otherwise.
  869. */
  870. int mont_mult(uint64_t* out, const uint64_t* a, const uint64_t *b, uint64_t *tmp, const MontContext *ctx)
  871. {
  872. if (NULL == out || NULL == a || NULL == b || NULL == tmp || NULL == ctx)
  873. return ERR_NULL;
  874. switch (ctx->modulus_type) {
  875. case ModulusP256:
  876. mont_mult_p256(out, a, b, ctx->modulus, ctx->m0, tmp, ctx->words);
  877. break;
  878. case ModulusP384:
  879. mont_mult_p384(out, a, b, ctx->modulus, ctx->m0, tmp, ctx->words);
  880. break;
  881. case ModulusP521:
  882. mont_mult_p521(out, a, b, ctx->modulus, ctx->m0, tmp, ctx->words);
  883. break;
  884. case ModulusGeneric:
  885. mont_mult_generic(out, a, b, ctx->modulus, ctx->m0, tmp, ctx->words);
  886. break;
  887. }
  888. return 0;
  889. }
  890. /*
  891. * Subtract integer b from a.
  892. *
  893. * @param out The location where the result is stored at; it must have been created with mont_number(&p,1,ctx).
  894. * It can be the same as either a or b.
  895. * @param a The number it will be subtracted from.
  896. * @param b The number to subtract.
  897. * @param tmp Temporary, internal result; it must have been created with mont_number(&p,2,ctx).
  898. * @param ctx The Montgomery context.
  899. * @return 0 for success, the relevant error code otherwise.
  900. */
  901. int mont_sub(uint64_t *out, const uint64_t *a, const uint64_t *b, uint64_t *tmp, const MontContext *ctx)
  902. {
  903. unsigned i;
  904. unsigned carry, borrow1 , borrow2;
  905. uint64_t *scratchpad;
  906. if (NULL == out || NULL == a || NULL == b || NULL == tmp || NULL == ctx)
  907. return ERR_NULL;
  908. scratchpad = tmp + ctx->words;
  909. /*
  910. * Compute difference in tmp[], and add modulus[]
  911. * to tmp[] into scratchpad[].
  912. */
  913. borrow2 = 0;
  914. carry = 0;
  915. for (i=0; i<ctx->words; i++) {
  916. borrow1 = b[i] > a[i];
  917. tmp[i] = a[i] - b[i];
  918. borrow1 |= borrow2 > tmp[i];
  919. tmp[i] -= borrow2;
  920. borrow2 = borrow1;
  921. scratchpad[i] = tmp[i] + carry;
  922. carry = scratchpad[i] < carry;
  923. scratchpad[i] += ctx->modulus[i];
  924. carry += scratchpad[i] < ctx->modulus[i];
  925. }
  926. /*
  927. * If there is no borrow, tmp[] is smaller than modulus.
  928. */
  929. mont_select(out, scratchpad, tmp, borrow2, ctx->words);
  930. return 0;
  931. }
  932. /*
  933. * Compute the modular inverse of an integer in Montgomery form.
  934. *
  935. * Condition: the modulus defining the Montgomery context MUST BE a non-secret prime number.
  936. *
  937. * @param out The location where the result will be stored at; it must have
  938. * been allocated with mont_number(&p, 1, ctx).
  939. * @param a The number to compute the modular inverse of, already in Montgomery form.
  940. * @param ctx The Montgomery context.
  941. * @return 0 for success, the relevant error code otherwise.
  942. */
  943. int mont_inv_prime(uint64_t *out, uint64_t *a, const MontContext *ctx)
  944. {
  945. unsigned idx_word;
  946. uint64_t bit;
  947. uint64_t *tmp1 = NULL;
  948. uint64_t *scratchpad = NULL;
  949. uint64_t *exponent = NULL;
  950. int res;
  951. if (NULL == out || NULL == a || NULL == ctx)
  952. return ERR_NULL;
  953. tmp1 = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  954. if (NULL == tmp1)
  955. return ERR_MEMORY;
  956. scratchpad = (uint64_t*)calloc(SCRATCHPAD_NR, ctx->words*sizeof(uint64_t));
  957. if (NULL == scratchpad) {
  958. res = ERR_MEMORY;
  959. goto cleanup;
  960. }
  961. /** Exponent is guaranteed to be >0 **/
  962. exponent = ctx->modulus_min_2;
  963. /* Find most significant bit */
  964. idx_word = ctx->words-1;
  965. for (;;) {
  966. if (exponent[idx_word] != 0)
  967. break;
  968. if (idx_word-- == 0)
  969. break;
  970. }
  971. for (bit = (uint64_t)1U << 63; 0 == (exponent[idx_word] & bit); bit>>=1);
  972. /* Start from 1 (in Montgomery form, which is R mod N) */
  973. memcpy(out, ctx->r_mod_n, ctx->bytes);
  974. /** Left-to-right exponentiation **/
  975. for (;;) {
  976. while (bit > 0) {
  977. mont_mult(tmp1, out, out, scratchpad, ctx);
  978. if (exponent[idx_word] & bit) {
  979. mont_mult(out, tmp1, a, scratchpad, ctx);
  980. } else {
  981. memcpy(out, tmp1, ctx->bytes);
  982. }
  983. bit >>= 1;
  984. }
  985. if (idx_word-- == 0)
  986. break;
  987. bit = (uint64_t)1 << 63;
  988. }
  989. res = 0;
  990. cleanup:
  991. free(tmp1);
  992. free(scratchpad);
  993. return res;
  994. }
  995. /*
  996. * Assign a value to a number in Montgomer form.
  997. *
  998. * @param out The location where the result is stored at; it must have been created with mont_number(&p,1,ctx).
  999. * @param x The value to set.
  1000. * @param ctx The Montgomery context.
  1001. * @return 0 for success, the relevant error code otherwise.
  1002. */
  1003. int mont_set(uint64_t *out, uint64_t x, const MontContext *ctx)
  1004. {
  1005. uint64_t *tmp, *scratchpad;
  1006. if (NULL == out || NULL == ctx)
  1007. return ERR_NULL;
  1008. if (x == 0) {
  1009. memset(out, 0, ctx->bytes);
  1010. return 0;
  1011. }
  1012. if (x == 1) {
  1013. mont_copy(out, ctx->r_mod_n, ctx);
  1014. return 0;
  1015. }
  1016. tmp = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1017. if (NULL == tmp)
  1018. return ERR_MEMORY;
  1019. tmp[0] = x;
  1020. scratchpad = (uint64_t*)calloc(SCRATCHPAD_NR, ctx->words*sizeof(uint64_t));
  1021. if (NULL == scratchpad) {
  1022. free(tmp);
  1023. return ERR_MEMORY;
  1024. }
  1025. if (ctx->modulus_type != ModulusP521)
  1026. mont_mult_generic(out, tmp, ctx->r2_mod_n, ctx->modulus, ctx->m0, scratchpad, ctx->words);
  1027. else
  1028. mont_copy(out, tmp, ctx);
  1029. free(tmp);
  1030. free(scratchpad);
  1031. return 0;
  1032. }
  1033. static int cmp_modulus(const uint8_t *mod1, size_t mod1_len, const uint8_t *mod2, size_t mod2_len)
  1034. {
  1035. size_t diff;
  1036. if (mod1_len > mod2_len) {
  1037. diff = mod1_len - mod2_len;
  1038. if (0 != memcmp(mod1+diff, mod2, mod2_len))
  1039. return -1;
  1040. if (NULL != memchr_not(mod1, 0, diff))
  1041. return -1;
  1042. } else {
  1043. diff = mod2_len - mod1_len;
  1044. if (0 != memcmp(mod2+diff, mod1, mod1_len))
  1045. return -1;
  1046. if (NULL != memchr_not(mod2, 0, diff))
  1047. return -1;
  1048. }
  1049. return 0;
  1050. }
  1051. /*
  1052. * Create a new context for the Montgomery and the given odd modulus.
  1053. *
  1054. * @param out The locate where the pointer to the newly allocated data will be stored at.
  1055. * The memory will contain the new Montgomery context.
  1056. * @param modulus The modulus encoded in big endian form.
  1057. * @param mod_len The length of the modulus in bytes.
  1058. * @return 0 for success, the appropriate code otherwise.
  1059. */
  1060. int mont_context_init(MontContext **out, const uint8_t *modulus, size_t mod_len)
  1061. {
  1062. const uint8_t p256_mod[32] = "\xff\xff\xff\xff\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff";
  1063. const uint8_t p384_mod[48] = "\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xfe\xff\xff\xff\xff\x00\x00\x00\x00\x00\x00\x00\x00\xff\xff\xff\xff";
  1064. const uint8_t p521_mod[66] = "\x01\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff";
  1065. uint64_t *scratchpad = NULL;
  1066. MontContext *ctx;
  1067. int res;
  1068. if (NULL == out || NULL == modulus)
  1069. return ERR_NULL;
  1070. /** Consume leading zeros **/
  1071. while (mod_len>0 && *modulus==0) {
  1072. modulus++;
  1073. mod_len--;
  1074. }
  1075. if (0 == mod_len)
  1076. return ERR_MODULUS;
  1077. /** Ensure modulus is odd and at least 3, otherwise we can't compute its inverse over B **/
  1078. if (is_even(modulus[mod_len-1]))
  1079. return ERR_MODULUS;
  1080. if (mod_len==1 && modulus[0]==1)
  1081. return ERR_MODULUS;
  1082. *out = ctx = (MontContext*)calloc(1, sizeof(MontContext));
  1083. if (NULL == ctx)
  1084. return ERR_MEMORY;
  1085. /* Check if the modulus has a special form */
  1086. /* For P-521, modulo reduction is very simple so the Montgomery
  1087. * representation is not actually used.
  1088. */
  1089. ctx->modulus_type = ModulusGeneric;
  1090. switch (mod_len) {
  1091. case sizeof(p256_mod):
  1092. if (0 == cmp_modulus(modulus, mod_len, p256_mod, sizeof(p256_mod))) {
  1093. ctx->modulus_type = ModulusP256;
  1094. }
  1095. break;
  1096. case sizeof(p384_mod):
  1097. if (0 == cmp_modulus(modulus, mod_len, p384_mod, sizeof(p384_mod))) {
  1098. ctx->modulus_type = ModulusP384;
  1099. }
  1100. break;
  1101. case sizeof(p521_mod):
  1102. if (0 == cmp_modulus(modulus, mod_len, p521_mod, sizeof(p521_mod))) {
  1103. ctx->modulus_type = ModulusP521;
  1104. }
  1105. break;
  1106. }
  1107. ctx->words = ((unsigned)mod_len + 7) / 8;
  1108. ctx->bytes = (unsigned)(ctx->words * sizeof(uint64_t));
  1109. ctx->modulus_len = (unsigned)mod_len;
  1110. /** Load modulus N **/
  1111. ctx->modulus = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1112. if (0 == ctx->modulus) {
  1113. res = ERR_MEMORY;
  1114. goto cleanup;
  1115. }
  1116. bytes_to_words(ctx->modulus, ctx->words, modulus, mod_len);
  1117. /** Prepare 1 **/
  1118. ctx->one = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1119. if (NULL == ctx->one) {
  1120. res = ERR_MEMORY;
  1121. goto cleanup;
  1122. }
  1123. ctx->one[0] = 1;
  1124. /** Pre-compute R² mod N **/
  1125. /** Pre-compute -n[0]⁻¹ mod R **/
  1126. ctx->r2_mod_n = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1127. if (0 == ctx->r2_mod_n) {
  1128. res = ERR_MEMORY;
  1129. goto cleanup;
  1130. }
  1131. if (ctx->modulus_type != ModulusP521) {
  1132. rsquare(ctx->r2_mod_n, ctx->modulus, ctx->words);
  1133. ctx->m0 = inverse64(~ctx->modulus[0]+1);
  1134. } else {
  1135. memcpy(ctx->r2_mod_n, ctx->one, ctx->words * sizeof(uint64_t));
  1136. ctx->m0 = 1U;
  1137. }
  1138. /** Pre-compute R mod N **/
  1139. ctx->r_mod_n = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1140. if (NULL == ctx->r_mod_n) {
  1141. res = ERR_MEMORY;
  1142. goto cleanup;
  1143. }
  1144. scratchpad = (uint64_t*)calloc(SCRATCHPAD_NR, ctx->words*sizeof(uint64_t));
  1145. if (NULL == scratchpad) {
  1146. res = ERR_MEMORY;
  1147. goto cleanup;
  1148. }
  1149. if (ctx->modulus_type != ModulusP521)
  1150. mont_mult_generic(ctx->r_mod_n, ctx->one, ctx->r2_mod_n, ctx->modulus, ctx->m0, scratchpad, ctx->words);
  1151. else
  1152. memcpy(ctx->r_mod_n, ctx->one, ctx->words * sizeof(uint64_t));
  1153. /** Pre-compute modulus - 2 **/
  1154. /** Modulus is guaranteed to be >= 3 **/
  1155. ctx->modulus_min_2 = (uint64_t*)calloc(ctx->words, sizeof(uint64_t));
  1156. if (NULL == ctx->modulus_min_2) {
  1157. res = ERR_MEMORY;
  1158. goto cleanup;
  1159. }
  1160. sub(ctx->modulus_min_2, ctx->modulus, ctx->one, ctx->words);
  1161. sub(ctx->modulus_min_2, ctx->modulus_min_2, ctx->one, ctx->words);
  1162. res = 0;
  1163. cleanup:
  1164. free(scratchpad);
  1165. if (res != 0) {
  1166. mont_context_free(ctx);
  1167. }
  1168. return res;
  1169. }
  1170. int mont_is_zero(const uint64_t *a, const MontContext *ctx)
  1171. {
  1172. unsigned i;
  1173. uint64_t sum = 0;
  1174. if (NULL == a || NULL == ctx)
  1175. return -1;
  1176. for (i=0; i<ctx->words; i++) {
  1177. sum |= *a++;
  1178. }
  1179. return (sum == 0);
  1180. }
  1181. int mont_is_one(const uint64_t *a, const MontContext *ctx)
  1182. {
  1183. unsigned i;
  1184. uint64_t sum = 0;
  1185. if (NULL == a || NULL == ctx)
  1186. return -1;
  1187. for (i=0; i<ctx->words; i++) {
  1188. sum |= a[i] ^ ctx->r_mod_n[i];
  1189. }
  1190. return (sum == 0);
  1191. }
  1192. int mont_is_equal(const uint64_t *a, const uint64_t *b, const MontContext *ctx)
  1193. {
  1194. unsigned i;
  1195. uint64_t result = 0;
  1196. if (NULL == a || NULL == b || NULL == ctx)
  1197. return -1;
  1198. for (i=0; i<ctx->words; i++) {
  1199. result |= *a++ ^ *b++;
  1200. }
  1201. return (result == 0);
  1202. }
  1203. int mont_copy(uint64_t *out, const uint64_t *a, const MontContext *ctx)
  1204. {
  1205. unsigned i;
  1206. if (NULL == out || NULL == a || NULL == ctx)
  1207. return ERR_NULL;
  1208. for (i=0; i<ctx->words; i++) {
  1209. *out++ = *a++;
  1210. }
  1211. return 0;
  1212. }