mpih-mul.c 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509
  1. // SPDX-License-Identifier: GPL-2.0-or-later
  2. /* mpihelp-mul.c - MPI helper functions
  3. * Copyright (C) 1994, 1996, 1998, 1999,
  4. * 2000 Free Software Foundation, Inc.
  5. *
  6. * This file is part of GnuPG.
  7. *
  8. * Note: This code is heavily based on the GNU MP Library.
  9. * Actually it's the same code with only minor changes in the
  10. * way the data is stored; this is to support the abstraction
  11. * of an optional secure memory allocation which may be used
  12. * to avoid revealing of sensitive data due to paging etc.
  13. * The GNU MP Library itself is published under the LGPL;
  14. * however I decided to publish this code under the plain GPL.
  15. */
  16. #include <linux/string.h>
  17. #include "mpi-internal.h"
  18. #include "longlong.h"
  19. #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
  20. do { \
  21. if ((size) < KARATSUBA_THRESHOLD) \
  22. mul_n_basecase(prodp, up, vp, size); \
  23. else \
  24. mul_n(prodp, up, vp, size, tspace); \
  25. } while (0);
  26. #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
  27. do { \
  28. if ((size) < KARATSUBA_THRESHOLD) \
  29. mpih_sqr_n_basecase(prodp, up, size); \
  30. else \
  31. mpih_sqr_n(prodp, up, size, tspace); \
  32. } while (0);
  33. /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
  34. * both with SIZE limbs, and store the result at PRODP. 2 * SIZE limbs are
  35. * always stored. Return the most significant limb.
  36. *
  37. * Argument constraints:
  38. * 1. PRODP != UP and PRODP != VP, i.e. the destination
  39. * must be distinct from the multiplier and the multiplicand.
  40. *
  41. *
  42. * Handle simple cases with traditional multiplication.
  43. *
  44. * This is the most critical code of multiplication. All multiplies rely
  45. * on this, both small and huge. Small ones arrive here immediately. Huge
  46. * ones arrive here as this is the base case for Karatsuba's recursive
  47. * algorithm below.
  48. */
  49. static mpi_limb_t
  50. mul_n_basecase(mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
  51. {
  52. mpi_size_t i;
  53. mpi_limb_t cy;
  54. mpi_limb_t v_limb;
  55. /* Multiply by the first limb in V separately, as the result can be
  56. * stored (not added) to PROD. We also avoid a loop for zeroing. */
  57. v_limb = vp[0];
  58. if (v_limb <= 1) {
  59. if (v_limb == 1)
  60. MPN_COPY(prodp, up, size);
  61. else
  62. MPN_ZERO(prodp, size);
  63. cy = 0;
  64. } else
  65. cy = mpihelp_mul_1(prodp, up, size, v_limb);
  66. prodp[size] = cy;
  67. prodp++;
  68. /* For each iteration in the outer loop, multiply one limb from
  69. * U with one limb from V, and add it to PROD. */
  70. for (i = 1; i < size; i++) {
  71. v_limb = vp[i];
  72. if (v_limb <= 1) {
  73. cy = 0;
  74. if (v_limb == 1)
  75. cy = mpihelp_add_n(prodp, prodp, up, size);
  76. } else
  77. cy = mpihelp_addmul_1(prodp, up, size, v_limb);
  78. prodp[size] = cy;
  79. prodp++;
  80. }
  81. return cy;
  82. }
  83. static void
  84. mul_n(mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp,
  85. mpi_size_t size, mpi_ptr_t tspace)
  86. {
  87. if (size & 1) {
  88. /* The size is odd, and the code below doesn't handle that.
  89. * Multiply the least significant (size - 1) limbs with a recursive
  90. * call, and handle the most significant limb of S1 and S2
  91. * separately.
  92. * A slightly faster way to do this would be to make the Karatsuba
  93. * code below behave as if the size were even, and let it check for
  94. * odd size in the end. I.e., in essence move this code to the end.
  95. * Doing so would save us a recursive call, and potentially make the
  96. * stack grow a lot less.
  97. */
  98. mpi_size_t esize = size - 1; /* even size */
  99. mpi_limb_t cy_limb;
  100. MPN_MUL_N_RECURSE(prodp, up, vp, esize, tspace);
  101. cy_limb = mpihelp_addmul_1(prodp + esize, up, esize, vp[esize]);
  102. prodp[esize + esize] = cy_limb;
  103. cy_limb = mpihelp_addmul_1(prodp + esize, vp, size, up[esize]);
  104. prodp[esize + size] = cy_limb;
  105. } else {
  106. /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
  107. *
  108. * Split U in two pieces, U1 and U0, such that
  109. * U = U0 + U1*(B**n),
  110. * and V in V1 and V0, such that
  111. * V = V0 + V1*(B**n).
  112. *
  113. * UV is then computed recursively using the identity
  114. *
  115. * 2n n n n
  116. * UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
  117. * 1 1 1 0 0 1 0 0
  118. *
  119. * Where B = 2**BITS_PER_MP_LIMB.
  120. */
  121. mpi_size_t hsize = size >> 1;
  122. mpi_limb_t cy;
  123. int negflg;
  124. /* Product H. ________________ ________________
  125. * |_____U1 x V1____||____U0 x V0_____|
  126. * Put result in upper part of PROD and pass low part of TSPACE
  127. * as new TSPACE.
  128. */
  129. MPN_MUL_N_RECURSE(prodp + size, up + hsize, vp + hsize, hsize,
  130. tspace);
  131. /* Product M. ________________
  132. * |_(U1-U0)(V0-V1)_|
  133. */
  134. if (mpihelp_cmp(up + hsize, up, hsize) >= 0) {
  135. mpihelp_sub_n(prodp, up + hsize, up, hsize);
  136. negflg = 0;
  137. } else {
  138. mpihelp_sub_n(prodp, up, up + hsize, hsize);
  139. negflg = 1;
  140. }
  141. if (mpihelp_cmp(vp + hsize, vp, hsize) >= 0) {
  142. mpihelp_sub_n(prodp + hsize, vp + hsize, vp, hsize);
  143. negflg ^= 1;
  144. } else {
  145. mpihelp_sub_n(prodp + hsize, vp, vp + hsize, hsize);
  146. /* No change of NEGFLG. */
  147. }
  148. /* Read temporary operands from low part of PROD.
  149. * Put result in low part of TSPACE using upper part of TSPACE
  150. * as new TSPACE.
  151. */
  152. MPN_MUL_N_RECURSE(tspace, prodp, prodp + hsize, hsize,
  153. tspace + size);
  154. /* Add/copy product H. */
  155. MPN_COPY(prodp + hsize, prodp + size, hsize);
  156. cy = mpihelp_add_n(prodp + size, prodp + size,
  157. prodp + size + hsize, hsize);
  158. /* Add product M (if NEGFLG M is a negative number) */
  159. if (negflg)
  160. cy -=
  161. mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace,
  162. size);
  163. else
  164. cy +=
  165. mpihelp_add_n(prodp + hsize, prodp + hsize, tspace,
  166. size);
  167. /* Product L. ________________ ________________
  168. * |________________||____U0 x V0_____|
  169. * Read temporary operands from low part of PROD.
  170. * Put result in low part of TSPACE using upper part of TSPACE
  171. * as new TSPACE.
  172. */
  173. MPN_MUL_N_RECURSE(tspace, up, vp, hsize, tspace + size);
  174. /* Add/copy Product L (twice) */
  175. cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
  176. if (cy)
  177. mpihelp_add_1(prodp + hsize + size,
  178. prodp + hsize + size, hsize, cy);
  179. MPN_COPY(prodp, tspace, hsize);
  180. cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize,
  181. hsize);
  182. if (cy)
  183. mpihelp_add_1(prodp + size, prodp + size, size, 1);
  184. }
  185. }
  186. void mpih_sqr_n_basecase(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size)
  187. {
  188. mpi_size_t i;
  189. mpi_limb_t cy_limb;
  190. mpi_limb_t v_limb;
  191. /* Multiply by the first limb in V separately, as the result can be
  192. * stored (not added) to PROD. We also avoid a loop for zeroing. */
  193. v_limb = up[0];
  194. if (v_limb <= 1) {
  195. if (v_limb == 1)
  196. MPN_COPY(prodp, up, size);
  197. else
  198. MPN_ZERO(prodp, size);
  199. cy_limb = 0;
  200. } else
  201. cy_limb = mpihelp_mul_1(prodp, up, size, v_limb);
  202. prodp[size] = cy_limb;
  203. prodp++;
  204. /* For each iteration in the outer loop, multiply one limb from
  205. * U with one limb from V, and add it to PROD. */
  206. for (i = 1; i < size; i++) {
  207. v_limb = up[i];
  208. if (v_limb <= 1) {
  209. cy_limb = 0;
  210. if (v_limb == 1)
  211. cy_limb = mpihelp_add_n(prodp, prodp, up, size);
  212. } else
  213. cy_limb = mpihelp_addmul_1(prodp, up, size, v_limb);
  214. prodp[size] = cy_limb;
  215. prodp++;
  216. }
  217. }
  218. void
  219. mpih_sqr_n(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size, mpi_ptr_t tspace)
  220. {
  221. if (size & 1) {
  222. /* The size is odd, and the code below doesn't handle that.
  223. * Multiply the least significant (size - 1) limbs with a recursive
  224. * call, and handle the most significant limb of S1 and S2
  225. * separately.
  226. * A slightly faster way to do this would be to make the Karatsuba
  227. * code below behave as if the size were even, and let it check for
  228. * odd size in the end. I.e., in essence move this code to the end.
  229. * Doing so would save us a recursive call, and potentially make the
  230. * stack grow a lot less.
  231. */
  232. mpi_size_t esize = size - 1; /* even size */
  233. mpi_limb_t cy_limb;
  234. MPN_SQR_N_RECURSE(prodp, up, esize, tspace);
  235. cy_limb = mpihelp_addmul_1(prodp + esize, up, esize, up[esize]);
  236. prodp[esize + esize] = cy_limb;
  237. cy_limb = mpihelp_addmul_1(prodp + esize, up, size, up[esize]);
  238. prodp[esize + size] = cy_limb;
  239. } else {
  240. mpi_size_t hsize = size >> 1;
  241. mpi_limb_t cy;
  242. /* Product H. ________________ ________________
  243. * |_____U1 x U1____||____U0 x U0_____|
  244. * Put result in upper part of PROD and pass low part of TSPACE
  245. * as new TSPACE.
  246. */
  247. MPN_SQR_N_RECURSE(prodp + size, up + hsize, hsize, tspace);
  248. /* Product M. ________________
  249. * |_(U1-U0)(U0-U1)_|
  250. */
  251. if (mpihelp_cmp(up + hsize, up, hsize) >= 0)
  252. mpihelp_sub_n(prodp, up + hsize, up, hsize);
  253. else
  254. mpihelp_sub_n(prodp, up, up + hsize, hsize);
  255. /* Read temporary operands from low part of PROD.
  256. * Put result in low part of TSPACE using upper part of TSPACE
  257. * as new TSPACE. */
  258. MPN_SQR_N_RECURSE(tspace, prodp, hsize, tspace + size);
  259. /* Add/copy product H */
  260. MPN_COPY(prodp + hsize, prodp + size, hsize);
  261. cy = mpihelp_add_n(prodp + size, prodp + size,
  262. prodp + size + hsize, hsize);
  263. /* Add product M (if NEGFLG M is a negative number). */
  264. cy -= mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace, size);
  265. /* Product L. ________________ ________________
  266. * |________________||____U0 x U0_____|
  267. * Read temporary operands from low part of PROD.
  268. * Put result in low part of TSPACE using upper part of TSPACE
  269. * as new TSPACE. */
  270. MPN_SQR_N_RECURSE(tspace, up, hsize, tspace + size);
  271. /* Add/copy Product L (twice). */
  272. cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
  273. if (cy)
  274. mpihelp_add_1(prodp + hsize + size,
  275. prodp + hsize + size, hsize, cy);
  276. MPN_COPY(prodp, tspace, hsize);
  277. cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize,
  278. hsize);
  279. if (cy)
  280. mpihelp_add_1(prodp + size, prodp + size, size, 1);
  281. }
  282. }
  283. void mpihelp_mul_n(mpi_ptr_t prodp,
  284. mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
  285. {
  286. if (up == vp) {
  287. if (size < KARATSUBA_THRESHOLD)
  288. mpih_sqr_n_basecase(prodp, up, size);
  289. else {
  290. mpi_ptr_t tspace;
  291. tspace = mpi_alloc_limb_space(2 * size);
  292. mpih_sqr_n(prodp, up, size, tspace);
  293. mpi_free_limb_space(tspace);
  294. }
  295. } else {
  296. if (size < KARATSUBA_THRESHOLD)
  297. mul_n_basecase(prodp, up, vp, size);
  298. else {
  299. mpi_ptr_t tspace;
  300. tspace = mpi_alloc_limb_space(2 * size);
  301. mul_n(prodp, up, vp, size, tspace);
  302. mpi_free_limb_space(tspace);
  303. }
  304. }
  305. }
  306. int
  307. mpihelp_mul_karatsuba_case(mpi_ptr_t prodp,
  308. mpi_ptr_t up, mpi_size_t usize,
  309. mpi_ptr_t vp, mpi_size_t vsize,
  310. struct karatsuba_ctx *ctx)
  311. {
  312. mpi_limb_t cy;
  313. if (!ctx->tspace || ctx->tspace_size < vsize) {
  314. if (ctx->tspace)
  315. mpi_free_limb_space(ctx->tspace);
  316. ctx->tspace = mpi_alloc_limb_space(2 * vsize);
  317. if (!ctx->tspace)
  318. return -ENOMEM;
  319. ctx->tspace_size = vsize;
  320. }
  321. MPN_MUL_N_RECURSE(prodp, up, vp, vsize, ctx->tspace);
  322. prodp += vsize;
  323. up += vsize;
  324. usize -= vsize;
  325. if (usize >= vsize) {
  326. if (!ctx->tp || ctx->tp_size < vsize) {
  327. if (ctx->tp)
  328. mpi_free_limb_space(ctx->tp);
  329. ctx->tp = mpi_alloc_limb_space(2 * vsize);
  330. if (!ctx->tp) {
  331. if (ctx->tspace)
  332. mpi_free_limb_space(ctx->tspace);
  333. ctx->tspace = NULL;
  334. return -ENOMEM;
  335. }
  336. ctx->tp_size = vsize;
  337. }
  338. do {
  339. MPN_MUL_N_RECURSE(ctx->tp, up, vp, vsize, ctx->tspace);
  340. cy = mpihelp_add_n(prodp, prodp, ctx->tp, vsize);
  341. mpihelp_add_1(prodp + vsize, ctx->tp + vsize, vsize,
  342. cy);
  343. prodp += vsize;
  344. up += vsize;
  345. usize -= vsize;
  346. } while (usize >= vsize);
  347. }
  348. if (usize) {
  349. if (usize < KARATSUBA_THRESHOLD) {
  350. mpi_limb_t tmp;
  351. if (mpihelp_mul(ctx->tspace, vp, vsize, up, usize, &tmp)
  352. < 0)
  353. return -ENOMEM;
  354. } else {
  355. if (!ctx->next) {
  356. ctx->next = kzalloc(sizeof *ctx, GFP_KERNEL);
  357. if (!ctx->next)
  358. return -ENOMEM;
  359. }
  360. if (mpihelp_mul_karatsuba_case(ctx->tspace,
  361. vp, vsize,
  362. up, usize,
  363. ctx->next) < 0)
  364. return -ENOMEM;
  365. }
  366. cy = mpihelp_add_n(prodp, prodp, ctx->tspace, vsize);
  367. mpihelp_add_1(prodp + vsize, ctx->tspace + vsize, usize, cy);
  368. }
  369. return 0;
  370. }
  371. void mpihelp_release_karatsuba_ctx(struct karatsuba_ctx *ctx)
  372. {
  373. struct karatsuba_ctx *ctx2;
  374. if (ctx->tp)
  375. mpi_free_limb_space(ctx->tp);
  376. if (ctx->tspace)
  377. mpi_free_limb_space(ctx->tspace);
  378. for (ctx = ctx->next; ctx; ctx = ctx2) {
  379. ctx2 = ctx->next;
  380. if (ctx->tp)
  381. mpi_free_limb_space(ctx->tp);
  382. if (ctx->tspace)
  383. mpi_free_limb_space(ctx->tspace);
  384. kfree(ctx);
  385. }
  386. }
  387. /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
  388. * and v (pointed to by VP, with VSIZE limbs), and store the result at
  389. * PRODP. USIZE + VSIZE limbs are always stored, but if the input
  390. * operands are normalized. Return the most significant limb of the
  391. * result.
  392. *
  393. * NOTE: The space pointed to by PRODP is overwritten before finished
  394. * with U and V, so overlap is an error.
  395. *
  396. * Argument constraints:
  397. * 1. USIZE >= VSIZE.
  398. * 2. PRODP != UP and PRODP != VP, i.e. the destination
  399. * must be distinct from the multiplier and the multiplicand.
  400. */
  401. int
  402. mpihelp_mul(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t usize,
  403. mpi_ptr_t vp, mpi_size_t vsize, mpi_limb_t *_result)
  404. {
  405. mpi_ptr_t prod_endp = prodp + usize + vsize - 1;
  406. mpi_limb_t cy;
  407. struct karatsuba_ctx ctx;
  408. if (vsize < KARATSUBA_THRESHOLD) {
  409. mpi_size_t i;
  410. mpi_limb_t v_limb;
  411. if (!vsize) {
  412. *_result = 0;
  413. return 0;
  414. }
  415. /* Multiply by the first limb in V separately, as the result can be
  416. * stored (not added) to PROD. We also avoid a loop for zeroing. */
  417. v_limb = vp[0];
  418. if (v_limb <= 1) {
  419. if (v_limb == 1)
  420. MPN_COPY(prodp, up, usize);
  421. else
  422. MPN_ZERO(prodp, usize);
  423. cy = 0;
  424. } else
  425. cy = mpihelp_mul_1(prodp, up, usize, v_limb);
  426. prodp[usize] = cy;
  427. prodp++;
  428. /* For each iteration in the outer loop, multiply one limb from
  429. * U with one limb from V, and add it to PROD. */
  430. for (i = 1; i < vsize; i++) {
  431. v_limb = vp[i];
  432. if (v_limb <= 1) {
  433. cy = 0;
  434. if (v_limb == 1)
  435. cy = mpihelp_add_n(prodp, prodp, up,
  436. usize);
  437. } else
  438. cy = mpihelp_addmul_1(prodp, up, usize, v_limb);
  439. prodp[usize] = cy;
  440. prodp++;
  441. }
  442. *_result = cy;
  443. return 0;
  444. }
  445. memset(&ctx, 0, sizeof ctx);
  446. if (mpihelp_mul_karatsuba_case(prodp, up, usize, vp, vsize, &ctx) < 0)
  447. return -ENOMEM;
  448. mpihelp_release_karatsuba_ctx(&ctx);
  449. *_result = *prod_endp;
  450. return 0;
  451. }