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

randombit / botan / 20279912191

16 Dec 2025 07:18PM UTC coverage: 90.518% (+0.2%) from 90.36%
20279912191

Pull #5167

github

web-flow
Merge b45d77a93 into 3d96b675e
Pull Request #5167: Changes to reduce unnecessary inclusions

101149 of 111745 relevant lines covered (90.52%)

12673176.65 hits per line

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

96.69
/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

20
namespace Botan {
21

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

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

35
   BOTAN_ARG_CHECK(p.is_odd(), "invalid prime");
2,922✔
36

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

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

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

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

53
   q -= 1;
359✔
54
   q >>= 1;
359✔
55

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

60
   if(n == 1) {
359✔
61
      return r;
172✔
62
   }
63

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

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

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

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

85
   while(n > 1) {
592✔
86
      q = n;
412✔
87

88
      size_t i = 0;
412✔
89
      while(q != 1) {
12,152✔
90
         q = mod_p.square(q);
11,746✔
91
         ++i;
11,746✔
92

93
         if(i >= s) {
11,746✔
94
            return BigInt::from_s32(-1);
6✔
95
         }
96
      }
97

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

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

109
   return r;
180✔
110
}
3,045✔
111

112
/*
113
* Calculate the Jacobi symbol
114
*
115
* See Algorithm 2.149 in Handbook of Applied Cryptography
116
*/
117
int32_t jacobi(BigInt a, BigInt n) {
101,696✔
118
   BOTAN_ARG_CHECK(n.is_odd() && n >= 3, "Argument n must be an odd integer >= 3");
305,088✔
119

120
   if(a < 0 || a >= n) {
172,287✔
121
      a %= n;
31,377✔
122
   }
123

124
   if(a == 0) {
101,696✔
125
      return 0;
126
   }
127
   if(a == 1) {
101,673✔
128
      return 1;
129
   }
130

131
   int32_t s = 1;
132

133
   for(;;) {
388,300✔
134
      const size_t e = low_zero_bits(a);
388,300✔
135
      a >>= e;
388,300✔
136
      const word n_mod_8 = n.word_at(0) % 8;
388,300✔
137
      const word n_mod_4 = n_mod_8 % 4;
388,300✔
138

139
      if(e % 2 == 1 && (n_mod_8 == 3 || n_mod_8 == 5)) {
388,300✔
140
         s = -s;
66,857✔
141
      }
142

143
      if(n_mod_4 == 3 && a % 4 == 3) {
388,300✔
144
         s = -s;
68,028✔
145
      }
146

147
      /*
148
      * The HAC presentation of the algorithm uses recursion, which is not
149
      * desirable or necessary.
150
      *
151
      * Instead we loop accumulating the product of the various jacobi()
152
      * subcomputations into s, until we reach algorithm termination, which
153
      * occurs in one of two ways.
154
      *
155
      * If a == 1 then the recursion has completed; we can return the value of s.
156
      *
157
      * Otherwise, after swapping and reducing, check for a == 0 [this value is
158
      * called `n1` in HAC's presentation]. This would imply that jacobi(n1,a1)
159
      * would have the value 0, due to Line 1 in HAC 2.149, in which case the
160
      * entire product is zero, and we can immediately return that result.
161
      */
162

163
      if(a == 1) {
388,300✔
164
         return s;
85,047✔
165
      }
166

167
      std::swap(a, n);
303,253✔
168

169
      BOTAN_ASSERT_NOMSG(n.is_odd());
606,506✔
170

171
      a %= n;
303,253✔
172

173
      if(a == 0) {
303,253✔
174
         return 0;
175
      }
176
   }
177
}
178

179
/*
180
* Square a BigInt
181
*/
182
BigInt square(const BigInt& x) {
1,021✔
183
   BigInt z = x;
1,021✔
184
   secure_vector<word> ws;
1,021✔
185
   z.square(ws);
1,021✔
186
   return z;
1,021✔
187
}
1,021✔
188

189
/*
190
* Return the number of 0 bits at the end of n
191
*/
192
size_t low_zero_bits(const BigInt& n) {
627,096✔
193
   size_t low_zero = 0;
627,096✔
194

195
   auto seen_nonempty_word = CT::Mask<word>::cleared();
627,096✔
196

197
   for(size_t i = 0; i != n.size(); ++i) {
6,642,880✔
198
      const word x = n.word_at(i);
6,015,784✔
199

200
      // ctz(0) will return sizeof(word)
201
      const size_t tz_x = ctz(x);
6,015,784✔
202

203
      // if x > 0 we want to count tz_x in total but not any
204
      // further words, so set the mask after the addition
205
      low_zero += seen_nonempty_word.if_not_set_return(tz_x);
6,015,784✔
206

207
      seen_nonempty_word |= CT::Mask<word>::expand(x);
6,015,784✔
208
   }
209

210
   // if we saw no words with x > 0 then n == 0 and the value we have
211
   // computed is meaningless. Instead return BigInt::zero() in that case.
212
   return static_cast<size_t>(seen_nonempty_word.if_set_return(low_zero));
627,096✔
213
}
214

215
/*
216
* Calculate the GCD in constant time
217
*/
218
BigInt gcd(const BigInt& a, const BigInt& b) {
26,615✔
219
   if(a.is_zero()) {
28,764✔
220
      return abs(b);
3✔
221
   }
222
   if(b.is_zero()) {
27,052✔
223
      return abs(a);
26,615✔
224
   }
225

226
   const size_t sz = std::max(a.sig_words(), b.sig_words());
26,610✔
227
   auto u = BigInt::with_capacity(sz);
26,610✔
228
   auto v = BigInt::with_capacity(sz);
26,610✔
229
   u += a;
26,610✔
230
   v += b;
26,610✔
231

232
   CT::poison_all(u, v);
26,610✔
233

234
   u.set_sign(BigInt::Positive);
26,610✔
235
   v.set_sign(BigInt::Positive);
26,610✔
236

237
   // In the worst case we have two fully populated big ints. After right
238
   // shifting so many times, we'll have reached the result for sure.
239
   const size_t loop_cnt = u.bits() + v.bits();
26,610✔
240

241
   // This temporary is big enough to hold all intermediate results of the
242
   // algorithm. No reallocation will happen during the loop.
243
   // Note however, that `ct_cond_assign()` will invalidate the 'sig_words'
244
   // cache, which _does not_ shrink the capacity of the underlying buffer.
245
   auto tmp = BigInt::with_capacity(sz);
26,610✔
246
   secure_vector<word> ws(sz * 2);
26,610✔
247
   size_t factors_of_two = 0;
26,610✔
248
   for(size_t i = 0; i != loop_cnt; ++i) {
5,705,179✔
249
      auto both_odd = CT::Mask<word>::expand_bool(u.is_odd()) & CT::Mask<word>::expand_bool(v.is_odd());
17,035,707✔
250

251
      // Subtract the smaller from the larger if both are odd
252
      auto u_gt_v = CT::Mask<word>::expand_bool(bigint_cmp(u._data(), u.size(), v._data(), v.size()) > 0);
5,678,569✔
253
      bigint_sub_abs(tmp.mutable_data(), u._data(), v._data(), sz, ws.data());
5,678,569✔
254
      u.ct_cond_assign((u_gt_v & both_odd).as_bool(), tmp);
5,678,569✔
255
      v.ct_cond_assign((~u_gt_v & both_odd).as_bool(), tmp);
5,678,569✔
256

257
      const auto u_is_even = CT::Mask<word>::expand_bool(u.is_even());
11,357,138✔
258
      const auto v_is_even = CT::Mask<word>::expand_bool(v.is_even());
11,357,138✔
259
      BOTAN_DEBUG_ASSERT((u_is_even | v_is_even).as_bool());
5,678,569✔
260

261
      // When both are even, we're going to eliminate a factor of 2.
262
      // We have to reapply this factor to the final result.
263
      factors_of_two += (u_is_even & v_is_even).if_set_return(1);
5,678,569✔
264

265
      // remove one factor of 2, if u is even
266
      bigint_shr2(tmp.mutable_data(), u._data(), sz, 1);
5,678,569✔
267
      u.ct_cond_assign(u_is_even.as_bool(), tmp);
5,678,569✔
268

269
      // remove one factor of 2, if v is even
270
      bigint_shr2(tmp.mutable_data(), v._data(), sz, 1);
5,678,569✔
271
      v.ct_cond_assign(v_is_even.as_bool(), tmp);
5,678,569✔
272
   }
273

274
   // The GCD (without factors of two) is either in u or v, the other one is
275
   // zero. The non-zero variable _must_ be odd, because all factors of two were
276
   // removed in the loop iterations above.
277
   BOTAN_DEBUG_ASSERT(u.is_zero() || v.is_zero());
26,610✔
278
   BOTAN_DEBUG_ASSERT(u.is_odd() || v.is_odd());
26,610✔
279

280
   // make sure that the GCD (without factors of two) is in u
281
   u.ct_cond_assign(u.is_even() /* .is_zero() would not be constant time */, v);
53,220✔
282

283
   // re-apply the factors of two
284
   u.ct_shift_left(factors_of_two);
26,610✔
285

286
   CT::unpoison_all(u, v);
26,610✔
287

288
   return u;
26,610✔
289
}
26,610✔
290

291
/*
292
* Calculate the LCM
293
*/
294
BigInt lcm(const BigInt& a, const BigInt& b) {
4,628✔
295
   if(a == b) {
4,628✔
296
      return a;
×
297
   }
298

299
   auto ab = a * b;
4,628✔
300
   ab.set_sign(BigInt::Positive);  // ignore the signs of a & b
4,628✔
301
   const auto g = gcd(a, b);
4,628✔
302
   return ct_divide(ab, g);
4,628✔
303
}
4,628✔
304

305
/*
306
* Modular Exponentiation
307
*/
308
BigInt power_mod(const BigInt& base, const BigInt& exp, const BigInt& mod) {
68✔
309
   if(mod.is_negative() || mod == 1) {
136✔
310
      return BigInt::zero();
×
311
   }
312

313
   if(base.is_zero() || mod.is_zero()) {
197✔
314
      if(exp.is_zero()) {
2✔
315
         return BigInt::one();
1✔
316
      }
317
      return BigInt::zero();
×
318
   }
319

320
   auto reduce_mod = Barrett_Reduction::for_secret_modulus(mod);
67✔
321

322
   const size_t exp_bits = exp.bits();
67✔
323

324
   if(mod.is_odd()) {
67✔
325
      const Montgomery_Params monty_params(mod, reduce_mod);
53✔
326
      return monty_exp(monty_params, ct_modulo(base, mod), exp, exp_bits).value();
53✔
327
   }
53✔
328

329
   /*
330
   Support for even modulus is just a convenience and not considered
331
   cryptographically important, so this implementation is slow ...
332
   */
333
   BigInt accum = BigInt::one();
14✔
334
   BigInt g = ct_modulo(base, mod);
14✔
335
   BigInt t;
14✔
336

337
   for(size_t i = 0; i != exp_bits; ++i) {
2,530✔
338
      t = reduce_mod.multiply(g, accum);
2,516✔
339
      g = reduce_mod.square(g);
2,516✔
340
      accum.ct_cond_assign(exp.get_bit(i), t);
5,032✔
341
   }
342
   return accum;
14✔
343
}
81✔
344

345
BigInt is_perfect_square(const BigInt& C) {
708✔
346
   if(C < 1) {
708✔
347
      throw Invalid_Argument("is_perfect_square requires C >= 1");
×
348
   }
349
   if(C == 1) {
708✔
350
      return BigInt::one();
×
351
   }
352

353
   const size_t n = C.bits();
708✔
354
   const size_t m = (n + 1) / 2;
708✔
355
   const BigInt B = C + BigInt::power_of_2(m);
708✔
356

357
   BigInt X = BigInt::power_of_2(m) - 1;
1,416✔
358
   BigInt X2 = (X * X);
708✔
359

360
   for(;;) {
1,677✔
361
      X = (X2 + C) / (2 * X);
1,677✔
362
      X2 = (X * X);
1,677✔
363

364
      if(X2 < B) {
1,677✔
365
         break;
366
      }
367
   }
368

369
   if(X2 == C) {
708✔
370
      return X;
63✔
371
   } else {
372
      return BigInt::zero();
645✔
373
   }
374
}
708✔
375

376
/*
377
* Test for primality using Miller-Rabin
378
*/
379
bool is_prime(const BigInt& n, RandomNumberGenerator& rng, size_t prob, bool is_random) {
53,777✔
380
   if(n == 2) {
53,777✔
381
      return true;
382
   }
383
   if(n <= 1 || n.is_even()) {
107,548✔
384
      return false;
385
   }
386

387
   const size_t n_bits = n.bits();
53,706✔
388

389
   // Fast path testing for small numbers (<= 65521)
390
   if(n_bits <= 16) {
53,706✔
391
      const uint16_t num = static_cast<uint16_t>(n.word_at(0));
32,802✔
392

393
      return std::binary_search(PRIMES, PRIMES + PRIME_TABLE_SIZE, num);
32,802✔
394
   }
395

396
   auto mod_n = Barrett_Reduction::for_secret_modulus(n);
20,904✔
397

398
   if(rng.is_seeded()) {
20,904✔
399
      const size_t t = miller_rabin_test_iterations(n_bits, prob, is_random);
20,904✔
400

401
      if(!is_miller_rabin_probable_prime(n, mod_n, rng, t)) {
20,904✔
402
         return false;
403
      }
404

405
      if(is_random) {
3,363✔
406
         return true;
407
      } else {
408
         return is_lucas_probable_prime(n, mod_n);
3,164✔
409
      }
410
   } else {
411
      return is_bailie_psw_probable_prime(n, mod_n);
×
412
   }
413
}
20,904✔
414

415
}  // 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