• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

randombit / botan / 16511306696

25 Jul 2025 01:07AM UTC coverage: 90.698% (+0.005%) from 90.693%
16511306696

push

github

web-flow
Merge pull request #5017 from randombit/jack/mp-cleanup

Clean up low level mp interfaces

99956 of 110208 relevant lines covered (90.7%)

12233642.11 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.17
/src/lib/math/numbertheory/numthry.cpp
1
/*
2
* Number Theory Functions
3
* (C) 1999-2011,2016,2018,2019 Jack Lloyd
4
* (C) 2007,2008 Falko Strenzke, FlexSecure GmbH
5
*
6
* Botan is released under the Simplified BSD License (see license.txt)
7
*/
8

9
#include <botan/numthry.h>
10

11
#include <botan/rng.h>
12
#include <botan/internal/barrett.h>
13
#include <botan/internal/ct_utils.h>
14
#include <botan/internal/divide.h>
15
#include <botan/internal/monty.h>
16
#include <botan/internal/monty_exp.h>
17
#include <botan/internal/mp_core.h>
18
#include <botan/internal/primality.h>
19
#include <algorithm>
20

21
namespace Botan {
22

23
/*
24
* Tonelli-Shanks algorithm
25
*/
26
BigInt sqrt_modulo_prime(const BigInt& a, const BigInt& p) {
1,450✔
27
   BOTAN_ARG_CHECK(p > 1, "invalid prime");
1,450✔
28
   BOTAN_ARG_CHECK(a < p, "value to solve for must be less than p");
1,450✔
29
   BOTAN_ARG_CHECK(a >= 0, "value to solve for must not be negative");
1,450✔
30

31
   // some very easy cases
32
   if(p == 2 || a <= 1) {
2,898✔
33
      return a;
3✔
34
   }
35

36
   BOTAN_ARG_CHECK(p.is_odd(), "invalid prime");
2,894✔
37

38
   if(jacobi(a, p) != 1) {  // not a quadratic residue
1,447✔
39
      return BigInt::from_s32(-1);
110✔
40
   }
41

42
   auto mod_p = Barrett_Reduction::for_public_modulus(p);
1,337✔
43
   const Montgomery_Params monty_p(p, mod_p);
1,337✔
44

45
   // If p == 3 (mod 4) there is a simple solution
46
   if(p % 4 == 3) {
1,337✔
47
      return monty_exp_vartime(monty_p, a, ((p + 1) >> 2)).value();
984✔
48
   }
49

50
   // Otherwise we have to use Shanks-Tonelli
51
   size_t s = low_zero_bits(p - 1);
353✔
52
   BigInt q = p >> s;
353✔
53

54
   q -= 1;
353✔
55
   q >>= 1;
353✔
56

57
   BigInt r = monty_exp_vartime(monty_p, a, q).value();
353✔
58
   BigInt n = mod_p.multiply(a, mod_p.square(r));
353✔
59
   r = mod_p.multiply(r, a);
353✔
60

61
   if(n == 1) {
353✔
62
      return r;
161✔
63
   }
64

65
   // find random quadratic nonresidue z
66
   word z = 2;
67
   for(;;) {
508✔
68
      if(jacobi(BigInt::from_word(z), p) == -1) {  // found one
508✔
69
         break;
70
      }
71

72
      z += 1;  // try next z
317✔
73

74
      /*
75
      * The expected number of tests to find a non-residue modulo a
76
      * prime is 2. If we have not found one after 256 then almost
77
      * certainly we have been given a non-prime p.
78
      */
79
      if(z >= 256) {
317✔
80
         return BigInt::from_s32(-1);
1✔
81
      }
82
   }
83

84
   BigInt c = monty_exp_vartime(monty_p, BigInt::from_word(z), (q << 1) + 1).value();
382✔
85

86
   while(n > 1) {
585✔
87
      q = n;
400✔
88

89
      size_t i = 0;
400✔
90
      while(q != 1) {
11,018✔
91
         q = mod_p.square(q);
10,624✔
92
         ++i;
10,624✔
93

94
         if(i >= s) {
10,624✔
95
            return BigInt::from_s32(-1);
6✔
96
         }
97
      }
98

99
      BOTAN_ASSERT_NOMSG(s >= (i + 1));  // No underflow!
394✔
100
      c = monty_exp_vartime(monty_p, c, BigInt::power_of_2(s - i - 1)).value();
394✔
101
      r = mod_p.multiply(r, c);
394✔
102
      c = mod_p.square(c);
394✔
103
      n = mod_p.multiply(n, c);
394✔
104

105
      // s decreases as the algorithm proceeds
106
      BOTAN_ASSERT_NOMSG(s >= i);
394✔
107
      s = i;
108
   }
109

110
   return r;
185✔
111
}
3,027✔
112

113
/*
114
* Calculate the Jacobi symbol
115
*/
116
int32_t jacobi(const BigInt& a, const BigInt& n) {
101,739✔
117
   if(n.is_even() || n < 2) {
305,217✔
118
      throw Invalid_Argument("jacobi: second argument must be odd and > 1");
×
119
   }
120

121
   BigInt x = a % n;
101,739✔
122
   BigInt y = n;
101,739✔
123
   int32_t J = 1;
101,739✔
124

125
   while(y > 1) {
410,064✔
126
      x %= y;
324,972✔
127
      if(x > y / 2) {
649,944✔
128
         x = y - x;
127,198✔
129
         if(y % 4 == 3) {
127,198✔
130
            J = -J;
56,978✔
131
         }
132
      }
133
      if(x.is_zero()) {
452,170✔
134
         return 0;
135
      }
136

137
      size_t shifts = low_zero_bits(x);
308,325✔
138
      x >>= shifts;
308,325✔
139
      if(shifts % 2 == 1) {
308,325✔
140
         word y_mod_8 = y % 8;
70,156✔
141
         if(y_mod_8 == 3 || y_mod_8 == 5) {
70,156✔
142
            J = -J;
43,133✔
143
         }
144
      }
145

146
      if(x % 4 == 3 && y % 4 == 3) {
308,325✔
147
         J = -J;
51,258✔
148
      }
149
      std::swap(x, y);
308,325✔
150
   }
151
   return J;
152
}
101,739✔
153

154
/*
155
* Square a BigInt
156
*/
157
BigInt square(const BigInt& x) {
1,021✔
158
   BigInt z = x;
1,021✔
159
   secure_vector<word> ws;
1,021✔
160
   z.square(ws);
1,021✔
161
   return z;
1,021✔
162
}
1,021✔
163

164
/*
165
* Return the number of 0 bits at the end of n
166
*/
167
size_t low_zero_bits(const BigInt& n) {
542,358✔
168
   size_t low_zero = 0;
542,358✔
169

170
   auto seen_nonempty_word = CT::Mask<word>::cleared();
542,358✔
171

172
   for(size_t i = 0; i != n.size(); ++i) {
5,788,731✔
173
      const word x = n.word_at(i);
5,246,373✔
174

175
      // ctz(0) will return sizeof(word)
176
      const size_t tz_x = ctz(x);
5,246,373✔
177

178
      // if x > 0 we want to count tz_x in total but not any
179
      // further words, so set the mask after the addition
180
      low_zero += seen_nonempty_word.if_not_set_return(tz_x);
5,246,373✔
181

182
      seen_nonempty_word |= CT::Mask<word>::expand(x);
5,246,373✔
183
   }
184

185
   // if we saw no words with x > 0 then n == 0 and the value we have
186
   // computed is meaningless. Instead return BigInt::zero() in that case.
187
   return static_cast<size_t>(seen_nonempty_word.if_set_return(low_zero));
542,358✔
188
}
189

190
/*
191
* Calculate the GCD in constant time
192
*/
193
BigInt gcd(const BigInt& a, const BigInt& b) {
26,631✔
194
   if(a.is_zero()) {
28,792✔
195
      return abs(b);
3✔
196
   }
197
   if(b.is_zero()) {
27,068✔
198
      return abs(a);
26,631✔
199
   }
200

201
   const size_t sz = std::max(a.sig_words(), b.sig_words());
26,626✔
202
   auto u = BigInt::with_capacity(sz);
26,626✔
203
   auto v = BigInt::with_capacity(sz);
26,626✔
204
   u += a;
26,626✔
205
   v += b;
26,626✔
206

207
   CT::poison_all(u, v);
26,626✔
208

209
   u.set_sign(BigInt::Positive);
26,626✔
210
   v.set_sign(BigInt::Positive);
26,626✔
211

212
   // In the worst case we have two fully populated big ints. After right
213
   // shifting so many times, we'll have reached the result for sure.
214
   const size_t loop_cnt = u.bits() + v.bits();
26,626✔
215

216
   // This temporary is big enough to hold all intermediate results of the
217
   // algorithm. No reallocation will happen during the loop.
218
   // Note however, that `ct_cond_assign()` will invalidate the 'sig_words'
219
   // cache, which _does not_ shrink the capacity of the underlying buffer.
220
   auto tmp = BigInt::with_capacity(sz);
26,626✔
221
   secure_vector<word> ws(sz * 2);
26,626✔
222
   size_t factors_of_two = 0;
26,626✔
223
   for(size_t i = 0; i != loop_cnt; ++i) {
5,718,110✔
224
      auto both_odd = CT::Mask<word>::expand_bool(u.is_odd()) & CT::Mask<word>::expand_bool(v.is_odd());
17,074,452✔
225

226
      // Subtract the smaller from the larger if both are odd
227
      auto u_gt_v = CT::Mask<word>::expand(bigint_cmp(u._data(), u.size(), v._data(), v.size()) > 0);
5,691,484✔
228
      bigint_sub_abs(tmp.mutable_data(), u._data(), v._data(), sz, ws.data());
5,691,484✔
229
      u.ct_cond_assign((u_gt_v & both_odd).as_bool(), tmp);
5,691,484✔
230
      v.ct_cond_assign((~u_gt_v & both_odd).as_bool(), tmp);
5,691,484✔
231

232
      const auto u_is_even = CT::Mask<word>::expand_bool(u.is_even());
11,382,968✔
233
      const auto v_is_even = CT::Mask<word>::expand_bool(v.is_even());
11,382,968✔
234
      BOTAN_DEBUG_ASSERT((u_is_even | v_is_even).as_bool());
5,691,484✔
235

236
      // When both are even, we're going to eliminate a factor of 2.
237
      // We have to reapply this factor to the final result.
238
      factors_of_two += (u_is_even & v_is_even).if_set_return(1);
5,691,484✔
239

240
      // remove one factor of 2, if u is even
241
      bigint_shr2(tmp.mutable_data(), u._data(), sz, 1);
5,691,484✔
242
      u.ct_cond_assign(u_is_even.as_bool(), tmp);
5,691,484✔
243

244
      // remove one factor of 2, if v is even
245
      bigint_shr2(tmp.mutable_data(), v._data(), sz, 1);
5,691,484✔
246
      v.ct_cond_assign(v_is_even.as_bool(), tmp);
5,691,484✔
247
   }
248

249
   // The GCD (without factors of two) is either in u or v, the other one is
250
   // zero. The non-zero variable _must_ be odd, because all factors of two were
251
   // removed in the loop iterations above.
252
   BOTAN_DEBUG_ASSERT(u.is_zero() || v.is_zero());
26,626✔
253
   BOTAN_DEBUG_ASSERT(u.is_odd() || v.is_odd());
26,626✔
254

255
   // make sure that the GCD (without factors of two) is in u
256
   u.ct_cond_assign(u.is_even() /* .is_zero() would not be constant time */, v);
53,252✔
257

258
   // re-apply the factors of two
259
   u.ct_shift_left(factors_of_two);
26,626✔
260

261
   CT::unpoison_all(u, v);
26,626✔
262

263
   return u;
26,626✔
264
}
26,626✔
265

266
/*
267
* Calculate the LCM
268
*/
269
BigInt lcm(const BigInt& a, const BigInt& b) {
4,630✔
270
   if(a == b) {
4,630✔
271
      return a;
×
272
   }
273

274
   auto ab = a * b;
4,630✔
275
   ab.set_sign(BigInt::Positive);  // ignore the signs of a & b
4,630✔
276
   const auto g = gcd(a, b);
4,630✔
277
   return ct_divide(ab, g);
4,630✔
278
}
4,630✔
279

280
/*
281
* Modular Exponentiation
282
*/
283
BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) {
65✔
284
   if(mod.is_negative() || mod == 1) {
130✔
285
      return BigInt::zero();
×
286
   }
287

288
   if(base.is_zero() || mod.is_zero()) {
188✔
289
      if(exp.is_zero()) {
2✔
290
         return BigInt::one();
1✔
291
      }
292
      return BigInt::zero();
×
293
   }
294

295
   auto reduce_mod = Barrett_Reduction::for_secret_modulus(mod);
64✔
296

297
   const size_t exp_bits = exp.bits();
64✔
298

299
   if(mod.is_odd()) {
64✔
300
      const Montgomery_Params monty_params(mod, reduce_mod);
49✔
301
      return monty_exp(monty_params, ct_modulo(base, mod), exp, exp_bits).value();
49✔
302
   }
49✔
303

304
   /*
305
   Support for even modulus is just a convenience and not considered
306
   cryptographically important, so this implementation is slow ...
307
   */
308
   BigInt accum = BigInt::one();
15✔
309
   BigInt g = ct_modulo(base, mod);
15✔
310
   BigInt t;
15✔
311

312
   for(size_t i = 0; i != exp_bits; ++i) {
2,882✔
313
      t = reduce_mod.multiply(g, accum);
2,867✔
314
      g = reduce_mod.square(g);
2,867✔
315
      accum.ct_cond_assign(exp.get_bit(i), t);
5,734✔
316
   }
317
   return accum;
15✔
318
}
79✔
319

320
BigInt is_perfect_square(const BigInt& C) {
705✔
321
   if(C < 1) {
705✔
322
      throw Invalid_Argument("is_perfect_square requires C >= 1");
×
323
   }
324
   if(C == 1) {
705✔
325
      return BigInt::one();
×
326
   }
327

328
   const size_t n = C.bits();
705✔
329
   const size_t m = (n + 1) / 2;
705✔
330
   const BigInt B = C + BigInt::power_of_2(m);
705✔
331

332
   BigInt X = BigInt::power_of_2(m) - 1;
1,410✔
333
   BigInt X2 = (X * X);
705✔
334

335
   for(;;) {
1,659✔
336
      X = (X2 + C) / (2 * X);
1,659✔
337
      X2 = (X * X);
1,659✔
338

339
      if(X2 < B) {
1,659✔
340
         break;
341
      }
342
   }
343

344
   if(X2 == C) {
705✔
345
      return X;
63✔
346
   } else {
347
      return BigInt::zero();
642✔
348
   }
349
}
705✔
350

351
/*
352
* Test for primality using Miller-Rabin
353
*/
354
bool is_prime(const BigInt& n, RandomNumberGenerator& rng, size_t prob, bool is_random) {
52,179✔
355
   if(n == 2) {
52,179✔
356
      return true;
357
   }
358
   if(n <= 1 || n.is_even()) {
104,352✔
359
      return false;
360
   }
361

362
   const size_t n_bits = n.bits();
52,108✔
363

364
   // Fast path testing for small numbers (<= 65521)
365
   if(n_bits <= 16) {
52,108✔
366
      const uint16_t num = static_cast<uint16_t>(n.word_at(0));
32,802✔
367

368
      return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
32,802✔
369
   }
370

371
   auto mod_n = Barrett_Reduction::for_secret_modulus(n);
19,306✔
372

373
   if(rng.is_seeded()) {
19,306✔
374
      const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
19,306✔
375

376
      if(is_miller_rabin_probable_prime(n, mod_n, rng, t) == false) {
19,306✔
377
         return false;
378
      }
379

380
      if(is_random) {
3,359✔
381
         return true;
382
      } else {
383
         return is_lucas_probable_prime(n, mod_n);
3,164✔
384
      }
385
   } else {
386
      return is_bailie_psw_probable_prime(n, mod_n);
×
387
   }
388
}
19,306✔
389

390
}  // namespace Botan
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc