mpih-div.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517
  1. // SPDX-License-Identifier: GPL-2.0-or-later
  2. /* mpihelp-div.c - MPI helper functions
  3. * Copyright (C) 1994, 1996 Free Software Foundation, Inc.
  4. * Copyright (C) 1998, 1999 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 "mpi-internal.h"
  17. #include "longlong.h"
  18. #ifndef UMUL_TIME
  19. #define UMUL_TIME 1
  20. #endif
  21. #ifndef UDIV_TIME
  22. #define UDIV_TIME UMUL_TIME
  23. #endif
  24. mpi_limb_t
  25. mpihelp_mod_1(mpi_ptr_t dividend_ptr, mpi_size_t dividend_size,
  26. mpi_limb_t divisor_limb)
  27. {
  28. mpi_size_t i;
  29. mpi_limb_t n1, n0, r;
  30. mpi_limb_t dummy __maybe_unused;
  31. /* Botch: Should this be handled at all? Rely on callers? */
  32. if (!dividend_size)
  33. return 0;
  34. /* If multiplication is much faster than division, and the
  35. * dividend is large, pre-invert the divisor, and use
  36. * only multiplications in the inner loop.
  37. *
  38. * This test should be read:
  39. * Does it ever help to use udiv_qrnnd_preinv?
  40. * && Does what we save compensate for the inversion overhead?
  41. */
  42. if (UDIV_TIME > (2 * UMUL_TIME + 6)
  43. && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) {
  44. int normalization_steps;
  45. normalization_steps = count_leading_zeros(divisor_limb);
  46. if (normalization_steps) {
  47. mpi_limb_t divisor_limb_inverted;
  48. divisor_limb <<= normalization_steps;
  49. /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
  50. * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
  51. * most significant bit (with weight 2**N) implicit.
  52. *
  53. * Special case for DIVISOR_LIMB == 100...000.
  54. */
  55. if (!(divisor_limb << 1))
  56. divisor_limb_inverted = ~(mpi_limb_t)0;
  57. else
  58. udiv_qrnnd(divisor_limb_inverted, dummy,
  59. -divisor_limb, 0, divisor_limb);
  60. n1 = dividend_ptr[dividend_size - 1];
  61. r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
  62. /* Possible optimization:
  63. * if (r == 0
  64. * && divisor_limb > ((n1 << normalization_steps)
  65. * | (dividend_ptr[dividend_size - 2] >> ...)))
  66. * ...one division less...
  67. */
  68. for (i = dividend_size - 2; i >= 0; i--) {
  69. n0 = dividend_ptr[i];
  70. UDIV_QRNND_PREINV(dummy, r, r,
  71. ((n1 << normalization_steps)
  72. | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
  73. divisor_limb, divisor_limb_inverted);
  74. n1 = n0;
  75. }
  76. UDIV_QRNND_PREINV(dummy, r, r,
  77. n1 << normalization_steps,
  78. divisor_limb, divisor_limb_inverted);
  79. return r >> normalization_steps;
  80. } else {
  81. mpi_limb_t divisor_limb_inverted;
  82. /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
  83. * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
  84. * most significant bit (with weight 2**N) implicit.
  85. *
  86. * Special case for DIVISOR_LIMB == 100...000.
  87. */
  88. if (!(divisor_limb << 1))
  89. divisor_limb_inverted = ~(mpi_limb_t)0;
  90. else
  91. udiv_qrnnd(divisor_limb_inverted, dummy,
  92. -divisor_limb, 0, divisor_limb);
  93. i = dividend_size - 1;
  94. r = dividend_ptr[i];
  95. if (r >= divisor_limb)
  96. r = 0;
  97. else
  98. i--;
  99. for ( ; i >= 0; i--) {
  100. n0 = dividend_ptr[i];
  101. UDIV_QRNND_PREINV(dummy, r, r,
  102. n0, divisor_limb, divisor_limb_inverted);
  103. }
  104. return r;
  105. }
  106. } else {
  107. if (UDIV_NEEDS_NORMALIZATION) {
  108. int normalization_steps;
  109. normalization_steps = count_leading_zeros(divisor_limb);
  110. if (normalization_steps) {
  111. divisor_limb <<= normalization_steps;
  112. n1 = dividend_ptr[dividend_size - 1];
  113. r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
  114. /* Possible optimization:
  115. * if (r == 0
  116. * && divisor_limb > ((n1 << normalization_steps)
  117. * | (dividend_ptr[dividend_size - 2] >> ...)))
  118. * ...one division less...
  119. */
  120. for (i = dividend_size - 2; i >= 0; i--) {
  121. n0 = dividend_ptr[i];
  122. udiv_qrnnd(dummy, r, r,
  123. ((n1 << normalization_steps)
  124. | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
  125. divisor_limb);
  126. n1 = n0;
  127. }
  128. udiv_qrnnd(dummy, r, r,
  129. n1 << normalization_steps,
  130. divisor_limb);
  131. return r >> normalization_steps;
  132. }
  133. }
  134. /* No normalization needed, either because udiv_qrnnd doesn't require
  135. * it, or because DIVISOR_LIMB is already normalized.
  136. */
  137. i = dividend_size - 1;
  138. r = dividend_ptr[i];
  139. if (r >= divisor_limb)
  140. r = 0;
  141. else
  142. i--;
  143. for (; i >= 0; i--) {
  144. n0 = dividend_ptr[i];
  145. udiv_qrnnd(dummy, r, r, n0, divisor_limb);
  146. }
  147. return r;
  148. }
  149. }
  150. /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
  151. * the NSIZE-DSIZE least significant quotient limbs at QP
  152. * and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
  153. * non-zero, generate that many fraction bits and append them after the
  154. * other quotient limbs.
  155. * Return the most significant limb of the quotient, this is always 0 or 1.
  156. *
  157. * Preconditions:
  158. * 0. NSIZE >= DSIZE.
  159. * 1. The most significant bit of the divisor must be set.
  160. * 2. QP must either not overlap with the input operands at all, or
  161. * QP + DSIZE >= NP must hold true. (This means that it's
  162. * possible to put the quotient in the high part of NUM, right after the
  163. * remainder in NUM.
  164. * 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero.
  165. */
  166. mpi_limb_t
  167. mpihelp_divrem(mpi_ptr_t qp, mpi_size_t qextra_limbs,
  168. mpi_ptr_t np, mpi_size_t nsize, mpi_ptr_t dp, mpi_size_t dsize)
  169. {
  170. mpi_limb_t most_significant_q_limb = 0;
  171. switch (dsize) {
  172. case 0:
  173. /* We are asked to divide by zero, so go ahead and do it! (To make
  174. the compiler not remove this statement, return the value.) */
  175. /*
  176. * existing clients of this function have been modified
  177. * not to call it with dsize == 0, so this should not happen
  178. */
  179. return 1 / dsize;
  180. case 1:
  181. {
  182. mpi_size_t i;
  183. mpi_limb_t n1;
  184. mpi_limb_t d;
  185. d = dp[0];
  186. n1 = np[nsize - 1];
  187. if (n1 >= d) {
  188. n1 -= d;
  189. most_significant_q_limb = 1;
  190. }
  191. qp += qextra_limbs;
  192. for (i = nsize - 2; i >= 0; i--)
  193. udiv_qrnnd(qp[i], n1, n1, np[i], d);
  194. qp -= qextra_limbs;
  195. for (i = qextra_limbs - 1; i >= 0; i--)
  196. udiv_qrnnd(qp[i], n1, n1, 0, d);
  197. np[0] = n1;
  198. }
  199. break;
  200. case 2:
  201. {
  202. mpi_size_t i;
  203. mpi_limb_t n1, n0, n2;
  204. mpi_limb_t d1, d0;
  205. np += nsize - 2;
  206. d1 = dp[1];
  207. d0 = dp[0];
  208. n1 = np[1];
  209. n0 = np[0];
  210. if (n1 >= d1 && (n1 > d1 || n0 >= d0)) {
  211. sub_ddmmss(n1, n0, n1, n0, d1, d0);
  212. most_significant_q_limb = 1;
  213. }
  214. for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--) {
  215. mpi_limb_t q;
  216. mpi_limb_t r;
  217. if (i >= qextra_limbs)
  218. np--;
  219. else
  220. np[0] = 0;
  221. if (n1 == d1) {
  222. /* Q should be either 111..111 or 111..110. Need special
  223. * treatment of this rare case as normal division would
  224. * give overflow. */
  225. q = ~(mpi_limb_t) 0;
  226. r = n0 + d1;
  227. if (r < d1) { /* Carry in the addition? */
  228. add_ssaaaa(n1, n0, r - d0,
  229. np[0], 0, d0);
  230. qp[i] = q;
  231. continue;
  232. }
  233. n1 = d0 - (d0 != 0 ? 1 : 0);
  234. n0 = -d0;
  235. } else {
  236. udiv_qrnnd(q, r, n1, n0, d1);
  237. umul_ppmm(n1, n0, d0, q);
  238. }
  239. n2 = np[0];
  240. q_test:
  241. if (n1 > r || (n1 == r && n0 > n2)) {
  242. /* The estimated Q was too large. */
  243. q--;
  244. sub_ddmmss(n1, n0, n1, n0, 0, d0);
  245. r += d1;
  246. if (r >= d1) /* If not carry, test Q again. */
  247. goto q_test;
  248. }
  249. qp[i] = q;
  250. sub_ddmmss(n1, n0, r, n2, n1, n0);
  251. }
  252. np[1] = n1;
  253. np[0] = n0;
  254. }
  255. break;
  256. default:
  257. {
  258. mpi_size_t i;
  259. mpi_limb_t dX, d1, n0;
  260. np += nsize - dsize;
  261. dX = dp[dsize - 1];
  262. d1 = dp[dsize - 2];
  263. n0 = np[dsize - 1];
  264. if (n0 >= dX) {
  265. if (n0 > dX
  266. || mpihelp_cmp(np, dp, dsize - 1) >= 0) {
  267. mpihelp_sub_n(np, np, dp, dsize);
  268. n0 = np[dsize - 1];
  269. most_significant_q_limb = 1;
  270. }
  271. }
  272. for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--) {
  273. mpi_limb_t q;
  274. mpi_limb_t n1, n2;
  275. mpi_limb_t cy_limb;
  276. if (i >= qextra_limbs) {
  277. np--;
  278. n2 = np[dsize];
  279. } else {
  280. n2 = np[dsize - 1];
  281. MPN_COPY_DECR(np + 1, np, dsize - 1);
  282. np[0] = 0;
  283. }
  284. if (n0 == dX) {
  285. /* This might over-estimate q, but it's probably not worth
  286. * the extra code here to find out. */
  287. q = ~(mpi_limb_t) 0;
  288. } else {
  289. mpi_limb_t r;
  290. udiv_qrnnd(q, r, n0, np[dsize - 1], dX);
  291. umul_ppmm(n1, n0, d1, q);
  292. while (n1 > r
  293. || (n1 == r
  294. && n0 > np[dsize - 2])) {
  295. q--;
  296. r += dX;
  297. if (r < dX) /* I.e. "carry in previous addition?" */
  298. break;
  299. n1 -= n0 < d1;
  300. n0 -= d1;
  301. }
  302. }
  303. /* Possible optimization: We already have (q * n0) and (1 * n1)
  304. * after the calculation of q. Taking advantage of that, we
  305. * could make this loop make two iterations less. */
  306. cy_limb = mpihelp_submul_1(np, dp, dsize, q);
  307. if (n2 != cy_limb) {
  308. mpihelp_add_n(np, np, dp, dsize);
  309. q--;
  310. }
  311. qp[i] = q;
  312. n0 = np[dsize - 1];
  313. }
  314. }
  315. }
  316. return most_significant_q_limb;
  317. }
  318. /****************
  319. * Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
  320. * Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
  321. * Return the single-limb remainder.
  322. * There are no constraints on the value of the divisor.
  323. *
  324. * QUOT_PTR and DIVIDEND_PTR might point to the same limb.
  325. */
  326. mpi_limb_t
  327. mpihelp_divmod_1(mpi_ptr_t quot_ptr,
  328. mpi_ptr_t dividend_ptr, mpi_size_t dividend_size,
  329. mpi_limb_t divisor_limb)
  330. {
  331. mpi_size_t i;
  332. mpi_limb_t n1, n0, r;
  333. mpi_limb_t dummy __maybe_unused;
  334. if (!dividend_size)
  335. return 0;
  336. /* If multiplication is much faster than division, and the
  337. * dividend is large, pre-invert the divisor, and use
  338. * only multiplications in the inner loop.
  339. *
  340. * This test should be read:
  341. * Does it ever help to use udiv_qrnnd_preinv?
  342. * && Does what we save compensate for the inversion overhead?
  343. */
  344. if (UDIV_TIME > (2 * UMUL_TIME + 6)
  345. && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME) {
  346. int normalization_steps;
  347. normalization_steps = count_leading_zeros(divisor_limb);
  348. if (normalization_steps) {
  349. mpi_limb_t divisor_limb_inverted;
  350. divisor_limb <<= normalization_steps;
  351. /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
  352. * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
  353. * most significant bit (with weight 2**N) implicit.
  354. */
  355. /* Special case for DIVISOR_LIMB == 100...000. */
  356. if (!(divisor_limb << 1))
  357. divisor_limb_inverted = ~(mpi_limb_t)0;
  358. else
  359. udiv_qrnnd(divisor_limb_inverted, dummy,
  360. -divisor_limb, 0, divisor_limb);
  361. n1 = dividend_ptr[dividend_size - 1];
  362. r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
  363. /* Possible optimization:
  364. * if (r == 0
  365. * && divisor_limb > ((n1 << normalization_steps)
  366. * | (dividend_ptr[dividend_size - 2] >> ...)))
  367. * ...one division less...
  368. */
  369. for (i = dividend_size - 2; i >= 0; i--) {
  370. n0 = dividend_ptr[i];
  371. UDIV_QRNND_PREINV(quot_ptr[i + 1], r, r,
  372. ((n1 << normalization_steps)
  373. | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
  374. divisor_limb, divisor_limb_inverted);
  375. n1 = n0;
  376. }
  377. UDIV_QRNND_PREINV(quot_ptr[0], r, r,
  378. n1 << normalization_steps,
  379. divisor_limb, divisor_limb_inverted);
  380. return r >> normalization_steps;
  381. } else {
  382. mpi_limb_t divisor_limb_inverted;
  383. /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
  384. * result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
  385. * most significant bit (with weight 2**N) implicit.
  386. */
  387. /* Special case for DIVISOR_LIMB == 100...000. */
  388. if (!(divisor_limb << 1))
  389. divisor_limb_inverted = ~(mpi_limb_t) 0;
  390. else
  391. udiv_qrnnd(divisor_limb_inverted, dummy,
  392. -divisor_limb, 0, divisor_limb);
  393. i = dividend_size - 1;
  394. r = dividend_ptr[i];
  395. if (r >= divisor_limb)
  396. r = 0;
  397. else
  398. quot_ptr[i--] = 0;
  399. for ( ; i >= 0; i--) {
  400. n0 = dividend_ptr[i];
  401. UDIV_QRNND_PREINV(quot_ptr[i], r, r,
  402. n0, divisor_limb, divisor_limb_inverted);
  403. }
  404. return r;
  405. }
  406. } else {
  407. if (UDIV_NEEDS_NORMALIZATION) {
  408. int normalization_steps;
  409. normalization_steps = count_leading_zeros(divisor_limb);
  410. if (normalization_steps) {
  411. divisor_limb <<= normalization_steps;
  412. n1 = dividend_ptr[dividend_size - 1];
  413. r = n1 >> (BITS_PER_MPI_LIMB - normalization_steps);
  414. /* Possible optimization:
  415. * if (r == 0
  416. * && divisor_limb > ((n1 << normalization_steps)
  417. * | (dividend_ptr[dividend_size - 2] >> ...)))
  418. * ...one division less...
  419. */
  420. for (i = dividend_size - 2; i >= 0; i--) {
  421. n0 = dividend_ptr[i];
  422. udiv_qrnnd(quot_ptr[i + 1], r, r,
  423. ((n1 << normalization_steps)
  424. | (n0 >> (BITS_PER_MPI_LIMB - normalization_steps))),
  425. divisor_limb);
  426. n1 = n0;
  427. }
  428. udiv_qrnnd(quot_ptr[0], r, r,
  429. n1 << normalization_steps,
  430. divisor_limb);
  431. return r >> normalization_steps;
  432. }
  433. }
  434. /* No normalization needed, either because udiv_qrnnd doesn't require
  435. * it, or because DIVISOR_LIMB is already normalized.
  436. */
  437. i = dividend_size - 1;
  438. r = dividend_ptr[i];
  439. if (r >= divisor_limb)
  440. r = 0;
  441. else
  442. quot_ptr[i--] = 0;
  443. for (; i >= 0; i--) {
  444. n0 = dividend_ptr[i];
  445. udiv_qrnnd(quot_ptr[i], r, r, n0, divisor_limb);
  446. }
  447. return r;
  448. }
  449. }