bn_div.c 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460
  1. /*
  2. * Copyright 1995-2022 The OpenSSL Project Authors. All Rights Reserved.
  3. *
  4. * Licensed under the Apache License 2.0 (the "License"). You may not use
  5. * this file except in compliance with the License. You can obtain a copy
  6. * in the file LICENSE in the source distribution or at
  7. * https://www.openssl.org/source/license.html
  8. */
  9. #include <assert.h>
  10. #include <openssl/bn.h>
  11. #include "internal/cryptlib.h"
  12. #include "bn_local.h"
  13. /* The old slow way */
  14. #if 0
  15. int BN_div(BIGNUM *dv, BIGNUM *rem, const BIGNUM *m, const BIGNUM *d,
  16. BN_CTX *ctx)
  17. {
  18. int i, nm, nd;
  19. int ret = 0;
  20. BIGNUM *D;
  21. bn_check_top(m);
  22. bn_check_top(d);
  23. if (BN_is_zero(d)) {
  24. ERR_raise(ERR_LIB_BN, BN_R_DIV_BY_ZERO);
  25. return 0;
  26. }
  27. if (BN_ucmp(m, d) < 0) {
  28. if (rem != NULL) {
  29. if (BN_copy(rem, m) == NULL)
  30. return 0;
  31. }
  32. if (dv != NULL)
  33. BN_zero(dv);
  34. return 1;
  35. }
  36. BN_CTX_start(ctx);
  37. D = BN_CTX_get(ctx);
  38. if (dv == NULL)
  39. dv = BN_CTX_get(ctx);
  40. if (rem == NULL)
  41. rem = BN_CTX_get(ctx);
  42. if (D == NULL || dv == NULL || rem == NULL)
  43. goto end;
  44. nd = BN_num_bits(d);
  45. nm = BN_num_bits(m);
  46. if (BN_copy(D, d) == NULL)
  47. goto end;
  48. if (BN_copy(rem, m) == NULL)
  49. goto end;
  50. /*
  51. * The next 2 are needed so we can do a dv->d[0]|=1 later since
  52. * BN_lshift1 will only work once there is a value :-)
  53. */
  54. BN_zero(dv);
  55. if (bn_wexpand(dv, 1) == NULL)
  56. goto end;
  57. dv->top = 1;
  58. if (!BN_lshift(D, D, nm - nd))
  59. goto end;
  60. for (i = nm - nd; i >= 0; i--) {
  61. if (!BN_lshift1(dv, dv))
  62. goto end;
  63. if (BN_ucmp(rem, D) >= 0) {
  64. dv->d[0] |= 1;
  65. if (!BN_usub(rem, rem, D))
  66. goto end;
  67. }
  68. /* CAN IMPROVE (and have now :=) */
  69. if (!BN_rshift1(D, D))
  70. goto end;
  71. }
  72. rem->neg = BN_is_zero(rem) ? 0 : m->neg;
  73. dv->neg = m->neg ^ d->neg;
  74. ret = 1;
  75. end:
  76. BN_CTX_end(ctx);
  77. return ret;
  78. }
  79. #else
  80. # if defined(BN_DIV3W)
  81. BN_ULONG bn_div_3_words(const BN_ULONG *m, BN_ULONG d1, BN_ULONG d0);
  82. # elif 0
  83. /*
  84. * This is #if-ed away, because it's a reference for assembly implementations,
  85. * where it can and should be made constant-time. But if you want to test it,
  86. * just replace 0 with 1.
  87. */
  88. # if BN_BITS2 == 64 && defined(__SIZEOF_INT128__) && __SIZEOF_INT128__==16
  89. # undef BN_ULLONG
  90. # define BN_ULLONG uint128_t
  91. # define BN_LLONG
  92. # endif
  93. # ifdef BN_LLONG
  94. # define BN_DIV3W
  95. /*
  96. * Interface is somewhat quirky, |m| is pointer to most significant limb,
  97. * and less significant limb is referred at |m[-1]|. This means that caller
  98. * is responsible for ensuring that |m[-1]| is valid. Second condition that
  99. * has to be met is that |d0|'s most significant bit has to be set. Or in
  100. * other words divisor has to be "bit-aligned to the left." bn_div_fixed_top
  101. * does all this. The subroutine considers four limbs, two of which are
  102. * "overlapping," hence the name...
  103. */
  104. static BN_ULONG bn_div_3_words(const BN_ULONG *m, BN_ULONG d1, BN_ULONG d0)
  105. {
  106. BN_ULLONG R = ((BN_ULLONG)m[0] << BN_BITS2) | m[-1];
  107. BN_ULLONG D = ((BN_ULLONG)d0 << BN_BITS2) | d1;
  108. BN_ULONG Q = 0, mask;
  109. int i;
  110. for (i = 0; i < BN_BITS2; i++) {
  111. Q <<= 1;
  112. if (R >= D) {
  113. Q |= 1;
  114. R -= D;
  115. }
  116. D >>= 1;
  117. }
  118. mask = 0 - (Q >> (BN_BITS2 - 1)); /* does it overflow? */
  119. Q <<= 1;
  120. Q |= (R >= D);
  121. return (Q | mask) & BN_MASK2;
  122. }
  123. # endif
  124. # endif
  125. static int bn_left_align(BIGNUM *num)
  126. {
  127. BN_ULONG *d = num->d, n, m, rmask;
  128. int top = num->top;
  129. int rshift = BN_num_bits_word(d[top - 1]), lshift, i;
  130. lshift = BN_BITS2 - rshift;
  131. rshift %= BN_BITS2; /* say no to undefined behaviour */
  132. rmask = (BN_ULONG)0 - rshift; /* rmask = 0 - (rshift != 0) */
  133. rmask |= rmask >> 8;
  134. for (i = 0, m = 0; i < top; i++) {
  135. n = d[i];
  136. d[i] = ((n << lshift) | m) & BN_MASK2;
  137. m = (n >> rshift) & rmask;
  138. }
  139. return lshift;
  140. }
  141. # if !defined(OPENSSL_NO_ASM) && !defined(OPENSSL_NO_INLINE_ASM) \
  142. && !defined(PEDANTIC) && !defined(BN_DIV3W)
  143. # if defined(__GNUC__) && __GNUC__>=2
  144. # if defined(__i386) || defined (__i386__)
  145. /*-
  146. * There were two reasons for implementing this template:
  147. * - GNU C generates a call to a function (__udivdi3 to be exact)
  148. * in reply to ((((BN_ULLONG)n0)<<BN_BITS2)|n1)/d0 (I fail to
  149. * understand why...);
  150. * - divl doesn't only calculate quotient, but also leaves
  151. * remainder in %edx which we can definitely use here:-)
  152. */
  153. # undef bn_div_words
  154. # define bn_div_words(n0,n1,d0) \
  155. ({ asm volatile ( \
  156. "divl %4" \
  157. : "=a"(q), "=d"(rem) \
  158. : "a"(n1), "d"(n0), "r"(d0) \
  159. : "cc"); \
  160. q; \
  161. })
  162. # define REMAINDER_IS_ALREADY_CALCULATED
  163. # elif defined(__x86_64) && defined(SIXTY_FOUR_BIT_LONG)
  164. /*
  165. * Same story here, but it's 128-bit by 64-bit division. Wow!
  166. */
  167. # undef bn_div_words
  168. # define bn_div_words(n0,n1,d0) \
  169. ({ asm volatile ( \
  170. "divq %4" \
  171. : "=a"(q), "=d"(rem) \
  172. : "a"(n1), "d"(n0), "r"(d0) \
  173. : "cc"); \
  174. q; \
  175. })
  176. # define REMAINDER_IS_ALREADY_CALCULATED
  177. # endif /* __<cpu> */
  178. # endif /* __GNUC__ */
  179. # endif /* OPENSSL_NO_ASM */
  180. /*-
  181. * BN_div computes dv := num / divisor, rounding towards
  182. * zero, and sets up rm such that dv*divisor + rm = num holds.
  183. * Thus:
  184. * dv->neg == num->neg ^ divisor->neg (unless the result is zero)
  185. * rm->neg == num->neg (unless the remainder is zero)
  186. * If 'dv' or 'rm' is NULL, the respective value is not returned.
  187. */
  188. int BN_div(BIGNUM *dv, BIGNUM *rm, const BIGNUM *num, const BIGNUM *divisor,
  189. BN_CTX *ctx)
  190. {
  191. int ret;
  192. if (BN_is_zero(divisor)) {
  193. ERR_raise(ERR_LIB_BN, BN_R_DIV_BY_ZERO);
  194. return 0;
  195. }
  196. /*
  197. * Invalid zero-padding would have particularly bad consequences so don't
  198. * just rely on bn_check_top() here (bn_check_top() works only for
  199. * BN_DEBUG builds)
  200. */
  201. if (divisor->d[divisor->top - 1] == 0) {
  202. ERR_raise(ERR_LIB_BN, BN_R_NOT_INITIALIZED);
  203. return 0;
  204. }
  205. ret = bn_div_fixed_top(dv, rm, num, divisor, ctx);
  206. if (ret) {
  207. if (dv != NULL)
  208. bn_correct_top(dv);
  209. if (rm != NULL)
  210. bn_correct_top(rm);
  211. }
  212. return ret;
  213. }
  214. /*
  215. * It's argued that *length* of *significant* part of divisor is public.
  216. * Even if it's private modulus that is. Again, *length* is assumed
  217. * public, but not *value*. Former is likely to be pre-defined by
  218. * algorithm with bit granularity, though below subroutine is invariant
  219. * of limb length. Thanks to this assumption we can require that |divisor|
  220. * may not be zero-padded, yet claim this subroutine "constant-time"(*).
  221. * This is because zero-padded dividend, |num|, is tolerated, so that
  222. * caller can pass dividend of public length(*), but with smaller amount
  223. * of significant limbs. This naturally means that quotient, |dv|, would
  224. * contain correspongly less significant limbs as well, and will be zero-
  225. * padded accordingly. Returned remainder, |rm|, will have same bit length
  226. * as divisor, also zero-padded if needed. These actually leave sign bits
  227. * in ambiguous state. In sense that we try to avoid negative zeros, while
  228. * zero-padded zeros would retain sign.
  229. *
  230. * (*) "Constant-time-ness" has two pre-conditions:
  231. *
  232. * - availability of constant-time bn_div_3_words;
  233. * - dividend is at least as "wide" as divisor, limb-wise, zero-padded
  234. * if so required, which shouldn't be a privacy problem, because
  235. * divisor's length is considered public;
  236. */
  237. int bn_div_fixed_top(BIGNUM *dv, BIGNUM *rm, const BIGNUM *num,
  238. const BIGNUM *divisor, BN_CTX *ctx)
  239. {
  240. int norm_shift, i, j, loop;
  241. BIGNUM *tmp, *snum, *sdiv, *res;
  242. BN_ULONG *resp, *wnum, *wnumtop;
  243. BN_ULONG d0, d1;
  244. int num_n, div_n, num_neg;
  245. assert(divisor->top > 0 && divisor->d[divisor->top - 1] != 0);
  246. bn_check_top(num);
  247. bn_check_top(divisor);
  248. bn_check_top(dv);
  249. bn_check_top(rm);
  250. BN_CTX_start(ctx);
  251. res = (dv == NULL) ? BN_CTX_get(ctx) : dv;
  252. tmp = BN_CTX_get(ctx);
  253. snum = BN_CTX_get(ctx);
  254. sdiv = BN_CTX_get(ctx);
  255. if (sdiv == NULL)
  256. goto err;
  257. /* First we normalise the numbers */
  258. if (!BN_copy(sdiv, divisor))
  259. goto err;
  260. norm_shift = bn_left_align(sdiv);
  261. sdiv->neg = 0;
  262. /*
  263. * Note that bn_lshift_fixed_top's output is always one limb longer
  264. * than input, even when norm_shift is zero. This means that amount of
  265. * inner loop iterations is invariant of dividend value, and that one
  266. * doesn't need to compare dividend and divisor if they were originally
  267. * of the same bit length.
  268. */
  269. if (!(bn_lshift_fixed_top(snum, num, norm_shift)))
  270. goto err;
  271. div_n = sdiv->top;
  272. num_n = snum->top;
  273. if (num_n <= div_n) {
  274. /* caller didn't pad dividend -> no constant-time guarantee... */
  275. if (bn_wexpand(snum, div_n + 1) == NULL)
  276. goto err;
  277. memset(&(snum->d[num_n]), 0, (div_n - num_n + 1) * sizeof(BN_ULONG));
  278. snum->top = num_n = div_n + 1;
  279. }
  280. loop = num_n - div_n;
  281. /*
  282. * Lets setup a 'window' into snum This is the part that corresponds to
  283. * the current 'area' being divided
  284. */
  285. wnum = &(snum->d[loop]);
  286. wnumtop = &(snum->d[num_n - 1]);
  287. /* Get the top 2 words of sdiv */
  288. d0 = sdiv->d[div_n - 1];
  289. d1 = (div_n == 1) ? 0 : sdiv->d[div_n - 2];
  290. /* Setup quotient */
  291. if (!bn_wexpand(res, loop))
  292. goto err;
  293. num_neg = num->neg;
  294. res->neg = (num_neg ^ divisor->neg);
  295. res->top = loop;
  296. res->flags |= BN_FLG_FIXED_TOP;
  297. resp = &(res->d[loop]);
  298. /* space for temp */
  299. if (!bn_wexpand(tmp, (div_n + 1)))
  300. goto err;
  301. for (i = 0; i < loop; i++, wnumtop--) {
  302. BN_ULONG q, l0;
  303. /*
  304. * the first part of the loop uses the top two words of snum and sdiv
  305. * to calculate a BN_ULONG q such that | wnum - sdiv * q | < sdiv
  306. */
  307. # if defined(BN_DIV3W)
  308. q = bn_div_3_words(wnumtop, d1, d0);
  309. # else
  310. BN_ULONG n0, n1, rem = 0;
  311. n0 = wnumtop[0];
  312. n1 = wnumtop[-1];
  313. if (n0 == d0)
  314. q = BN_MASK2;
  315. else { /* n0 < d0 */
  316. BN_ULONG n2 = (wnumtop == wnum) ? 0 : wnumtop[-2];
  317. # ifdef BN_LLONG
  318. BN_ULLONG t2;
  319. # if defined(BN_LLONG) && defined(BN_DIV2W) && !defined(bn_div_words)
  320. q = (BN_ULONG)(((((BN_ULLONG) n0) << BN_BITS2) | n1) / d0);
  321. # else
  322. q = bn_div_words(n0, n1, d0);
  323. # endif
  324. # ifndef REMAINDER_IS_ALREADY_CALCULATED
  325. /*
  326. * rem doesn't have to be BN_ULLONG. The least we
  327. * know it's less that d0, isn't it?
  328. */
  329. rem = (n1 - q * d0) & BN_MASK2;
  330. # endif
  331. t2 = (BN_ULLONG) d1 *q;
  332. for (;;) {
  333. if (t2 <= ((((BN_ULLONG) rem) << BN_BITS2) | n2))
  334. break;
  335. q--;
  336. rem += d0;
  337. if (rem < d0)
  338. break; /* don't let rem overflow */
  339. t2 -= d1;
  340. }
  341. # else /* !BN_LLONG */
  342. BN_ULONG t2l, t2h;
  343. q = bn_div_words(n0, n1, d0);
  344. # ifndef REMAINDER_IS_ALREADY_CALCULATED
  345. rem = (n1 - q * d0) & BN_MASK2;
  346. # endif
  347. # if defined(BN_UMULT_LOHI)
  348. BN_UMULT_LOHI(t2l, t2h, d1, q);
  349. # elif defined(BN_UMULT_HIGH)
  350. t2l = d1 * q;
  351. t2h = BN_UMULT_HIGH(d1, q);
  352. # else
  353. {
  354. BN_ULONG ql, qh;
  355. t2l = LBITS(d1);
  356. t2h = HBITS(d1);
  357. ql = LBITS(q);
  358. qh = HBITS(q);
  359. mul64(t2l, t2h, ql, qh); /* t2=(BN_ULLONG)d1*q; */
  360. }
  361. # endif
  362. for (;;) {
  363. if ((t2h < rem) || ((t2h == rem) && (t2l <= n2)))
  364. break;
  365. q--;
  366. rem += d0;
  367. if (rem < d0)
  368. break; /* don't let rem overflow */
  369. if (t2l < d1)
  370. t2h--;
  371. t2l -= d1;
  372. }
  373. # endif /* !BN_LLONG */
  374. }
  375. # endif /* !BN_DIV3W */
  376. l0 = bn_mul_words(tmp->d, sdiv->d, div_n, q);
  377. tmp->d[div_n] = l0;
  378. wnum--;
  379. /*
  380. * ignore top values of the bignums just sub the two BN_ULONG arrays
  381. * with bn_sub_words
  382. */
  383. l0 = bn_sub_words(wnum, wnum, tmp->d, div_n + 1);
  384. q -= l0;
  385. /*
  386. * Note: As we have considered only the leading two BN_ULONGs in
  387. * the calculation of q, sdiv * q might be greater than wnum (but
  388. * then (q-1) * sdiv is less or equal than wnum)
  389. */
  390. for (l0 = 0 - l0, j = 0; j < div_n; j++)
  391. tmp->d[j] = sdiv->d[j] & l0;
  392. l0 = bn_add_words(wnum, wnum, tmp->d, div_n);
  393. (*wnumtop) += l0;
  394. assert((*wnumtop) == 0);
  395. /* store part of the result */
  396. *--resp = q;
  397. }
  398. /* snum holds remainder, it's as wide as divisor */
  399. snum->neg = num_neg;
  400. snum->top = div_n;
  401. snum->flags |= BN_FLG_FIXED_TOP;
  402. if (rm != NULL && bn_rshift_fixed_top(rm, snum, norm_shift) == 0)
  403. goto err;
  404. BN_CTX_end(ctx);
  405. return 1;
  406. err:
  407. bn_check_top(rm);
  408. BN_CTX_end(ctx);
  409. return 0;
  410. }
  411. #endif