mpih-mul.c 15 KB

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