bn_mp_prime_next_prime.c 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132
  1. #include "tommath_private.h"
  2. #ifdef BN_MP_PRIME_NEXT_PRIME_C
  3. /* LibTomMath, multiple-precision integer library -- Tom St Denis */
  4. /* SPDX-License-Identifier: Unlicense */
  5. /* finds the next prime after the number "a" using "t" trials
  6. * of Miller-Rabin.
  7. *
  8. * bbs_style = 1 means the prime must be congruent to 3 mod 4
  9. */
  10. mp_err mp_prime_next_prime(mp_int *a, int t, int bbs_style)
  11. {
  12. int x, y;
  13. mp_ord cmp;
  14. mp_err err;
  15. mp_bool res = MP_NO;
  16. mp_digit res_tab[PRIVATE_MP_PRIME_TAB_SIZE], step, kstep;
  17. mp_int b;
  18. /* force positive */
  19. a->sign = MP_ZPOS;
  20. /* simple algo if a is less than the largest prime in the table */
  21. if (mp_cmp_d(a, s_mp_prime_tab[PRIVATE_MP_PRIME_TAB_SIZE-1]) == MP_LT) {
  22. /* find which prime it is bigger than "a" */
  23. for (x = 0; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
  24. cmp = mp_cmp_d(a, s_mp_prime_tab[x]);
  25. if (cmp == MP_EQ) {
  26. continue;
  27. }
  28. if (cmp != MP_GT) {
  29. if ((bbs_style == 1) && ((s_mp_prime_tab[x] & 3u) != 3u)) {
  30. /* try again until we get a prime congruent to 3 mod 4 */
  31. continue;
  32. } else {
  33. mp_set(a, s_mp_prime_tab[x]);
  34. return MP_OKAY;
  35. }
  36. }
  37. }
  38. /* fall through to the sieve */
  39. }
  40. /* generate a prime congruent to 3 mod 4 or 1/3 mod 4? */
  41. if (bbs_style == 1) {
  42. kstep = 4;
  43. } else {
  44. kstep = 2;
  45. }
  46. /* at this point we will use a combination of a sieve and Miller-Rabin */
  47. if (bbs_style == 1) {
  48. /* if a mod 4 != 3 subtract the correct value to make it so */
  49. if ((a->dp[0] & 3u) != 3u) {
  50. if ((err = mp_sub_d(a, (a->dp[0] & 3u) + 1u, a)) != MP_OKAY) {
  51. return err;
  52. }
  53. }
  54. } else {
  55. if (MP_IS_EVEN(a)) {
  56. /* force odd */
  57. if ((err = mp_sub_d(a, 1uL, a)) != MP_OKAY) {
  58. return err;
  59. }
  60. }
  61. }
  62. /* generate the restable */
  63. for (x = 1; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
  64. if ((err = mp_mod_d(a, s_mp_prime_tab[x], res_tab + x)) != MP_OKAY) {
  65. return err;
  66. }
  67. }
  68. /* init temp used for Miller-Rabin Testing */
  69. if ((err = mp_init(&b)) != MP_OKAY) {
  70. return err;
  71. }
  72. for (;;) {
  73. /* skip to the next non-trivially divisible candidate */
  74. step = 0;
  75. do {
  76. /* y == 1 if any residue was zero [e.g. cannot be prime] */
  77. y = 0;
  78. /* increase step to next candidate */
  79. step += kstep;
  80. /* compute the new residue without using division */
  81. for (x = 1; x < PRIVATE_MP_PRIME_TAB_SIZE; x++) {
  82. /* add the step to each residue */
  83. res_tab[x] += kstep;
  84. /* subtract the modulus [instead of using division] */
  85. if (res_tab[x] >= s_mp_prime_tab[x]) {
  86. res_tab[x] -= s_mp_prime_tab[x];
  87. }
  88. /* set flag if zero */
  89. if (res_tab[x] == 0u) {
  90. y = 1;
  91. }
  92. }
  93. } while ((y == 1) && (step < (((mp_digit)1 << MP_DIGIT_BIT) - kstep)));
  94. /* add the step */
  95. if ((err = mp_add_d(a, step, a)) != MP_OKAY) {
  96. goto LBL_ERR;
  97. }
  98. /* if didn't pass sieve and step == MP_MAX then skip test */
  99. if ((y == 1) && (step >= (((mp_digit)1 << MP_DIGIT_BIT) - kstep))) {
  100. continue;
  101. }
  102. if ((err = mp_prime_is_prime(a, t, &res)) != MP_OKAY) {
  103. goto LBL_ERR;
  104. }
  105. if (res == MP_YES) {
  106. break;
  107. }
  108. }
  109. err = MP_OKAY;
  110. LBL_ERR:
  111. mp_clear(&b);
  112. return err;
  113. }
  114. #endif