pierrelf.com
we (web engine): Experimental web browser project to understand the limits of Claude
1//! Arbitrary-precision unsigned integer arithmetic for RSA.
2//!
3//! Provides a `BigUint` type backed by little-endian `u64` limbs, with
4//! modular exponentiation (square-and-multiply) for up to 4096-bit keys.
5
6use core::cmp::Ordering;
7
8/// Arbitrary-precision unsigned integer stored as little-endian `u64` limbs.
9#[derive(Clone, Debug)]
10pub struct BigUint {
11 /// Limbs in little-endian order (limbs[0] is the least significant).
12 limbs: Vec<u64>,
13}
14
15impl BigUint {
16 /// The number zero.
17 pub fn zero() -> Self {
18 Self { limbs: vec![0] }
19 }
20
21 /// The number one.
22 pub fn one() -> Self {
23 Self { limbs: vec![1] }
24 }
25
26 /// Create from a single u64.
27 pub fn from_u64(v: u64) -> Self {
28 Self { limbs: vec![v] }
29 }
30
31 /// Create from big-endian bytes (as used in DER INTEGER encoding).
32 pub fn from_be_bytes(bytes: &[u8]) -> Self {
33 if bytes.is_empty() {
34 return Self::zero();
35 }
36
37 // Convert big-endian bytes to little-endian u64 limbs.
38 let mut limbs = Vec::new();
39 let mut i = bytes.len();
40 while i > 0 {
41 let start = i.saturating_sub(8);
42 let chunk = &bytes[start..i];
43 let mut val = 0u64;
44 for &b in chunk {
45 val = (val << 8) | b as u64;
46 }
47 limbs.push(val);
48 i = start;
49 }
50
51 let mut result = Self { limbs };
52 result.normalize();
53 result
54 }
55
56 /// Export to big-endian bytes.
57 pub fn to_be_bytes(&self) -> Vec<u8> {
58 if self.is_zero() {
59 return vec![0];
60 }
61
62 let mut bytes = Vec::new();
63 // Start from the most significant limb.
64 let mut started = false;
65 for &limb in self.limbs.iter().rev() {
66 for shift in (0..8).rev() {
67 let byte = (limb >> (shift * 8)) as u8;
68 if !started && byte == 0 {
69 continue;
70 }
71 started = true;
72 bytes.push(byte);
73 }
74 }
75
76 if bytes.is_empty() {
77 bytes.push(0);
78 }
79 bytes
80 }
81
82 /// Export to big-endian bytes, zero-padded to exactly `len` bytes.
83 pub fn to_be_bytes_padded(&self, len: usize) -> Vec<u8> {
84 let raw = self.to_be_bytes();
85 if raw.len() >= len {
86 // Take the last `len` bytes (truncate leading bytes).
87 return raw[raw.len() - len..].to_vec();
88 }
89 let mut out = vec![0u8; len - raw.len()];
90 out.extend_from_slice(&raw);
91 out
92 }
93
94 /// True if the value is zero.
95 pub fn is_zero(&self) -> bool {
96 self.limbs.iter().all(|&l| l == 0)
97 }
98
99 /// Number of significant bits.
100 pub fn bit_len(&self) -> usize {
101 if self.is_zero() {
102 return 0;
103 }
104 let top = self.limbs.len() - 1;
105 let top_bits = 64 - self.limbs[top].leading_zeros() as usize;
106 top * 64 + top_bits
107 }
108
109 /// Get bit at position `i` (0-indexed from LSB).
110 pub fn bit(&self, i: usize) -> bool {
111 let limb_idx = i / 64;
112 let bit_idx = i % 64;
113 if limb_idx >= self.limbs.len() {
114 return false;
115 }
116 (self.limbs[limb_idx] >> bit_idx) & 1 == 1
117 }
118
119 /// Remove trailing zero limbs (keep at least one limb).
120 fn normalize(&mut self) {
121 while self.limbs.len() > 1 && *self.limbs.last().unwrap() == 0 {
122 self.limbs.pop();
123 }
124 }
125
126 /// Addition: self + other.
127 pub fn add(&self, other: &Self) -> Self {
128 let max_len = self.limbs.len().max(other.limbs.len());
129 let mut result = Vec::with_capacity(max_len + 1);
130 let mut carry = 0u64;
131
132 for i in 0..max_len {
133 let a = if i < self.limbs.len() {
134 self.limbs[i]
135 } else {
136 0
137 };
138 let b = if i < other.limbs.len() {
139 other.limbs[i]
140 } else {
141 0
142 };
143 let (sum1, c1) = a.overflowing_add(b);
144 let (sum2, c2) = sum1.overflowing_add(carry);
145 result.push(sum2);
146 carry = (c1 as u64) + (c2 as u64);
147 }
148
149 if carry > 0 {
150 result.push(carry);
151 }
152
153 let mut r = Self { limbs: result };
154 r.normalize();
155 r
156 }
157
158 /// Subtraction: self - other. Panics if other > self.
159 pub fn sub(&self, other: &Self) -> Self {
160 debug_assert!(self.cmp(other) != Ordering::Less);
161 let mut result = Vec::with_capacity(self.limbs.len());
162 let mut borrow = 0u64;
163
164 for i in 0..self.limbs.len() {
165 let a = self.limbs[i];
166 let b = if i < other.limbs.len() {
167 other.limbs[i]
168 } else {
169 0
170 };
171 let (diff1, b1) = a.overflowing_sub(b);
172 let (diff2, b2) = diff1.overflowing_sub(borrow);
173 result.push(diff2);
174 borrow = (b1 as u64) + (b2 as u64);
175 }
176
177 let mut r = Self { limbs: result };
178 r.normalize();
179 r
180 }
181
182 /// Multiplication: self * other.
183 pub fn mul(&self, other: &Self) -> Self {
184 let n = self.limbs.len();
185 let m = other.limbs.len();
186 let mut result = vec![0u64; n + m];
187
188 for i in 0..n {
189 let mut carry = 0u128;
190 for j in 0..m {
191 let prod = (self.limbs[i] as u128) * (other.limbs[j] as u128)
192 + result[i + j] as u128
193 + carry;
194 result[i + j] = prod as u64;
195 carry = prod >> 64;
196 }
197 result[i + m] = carry as u64;
198 }
199
200 let mut r = Self { limbs: result };
201 r.normalize();
202 r
203 }
204
205 /// Division with remainder: returns (quotient, remainder).
206 pub fn div_rem(&self, divisor: &Self) -> (Self, Self) {
207 assert!(!divisor.is_zero(), "division by zero");
208
209 if self.cmp(divisor) == Ordering::Less {
210 return (Self::zero(), self.clone());
211 }
212
213 if divisor.limbs.len() == 1 {
214 return self.div_rem_single(divisor.limbs[0]);
215 }
216
217 self.div_rem_knuth(divisor)
218 }
219
220 /// Single-limb division.
221 fn div_rem_single(&self, d: u64) -> (Self, Self) {
222 let mut quotient = vec![0u64; self.limbs.len()];
223 let mut rem = 0u128;
224
225 for i in (0..self.limbs.len()).rev() {
226 rem = (rem << 64) | self.limbs[i] as u128;
227 quotient[i] = (rem / d as u128) as u64;
228 rem %= d as u128;
229 }
230
231 let mut q = Self { limbs: quotient };
232 q.normalize();
233 (q, Self::from_u64(rem as u64))
234 }
235
236 /// Multi-limb division using Knuth's Algorithm D.
237 fn div_rem_knuth(&self, divisor: &Self) -> (Self, Self) {
238 let n = divisor.limbs.len();
239 let m = self.limbs.len() - n;
240
241 // Normalize: shift so that the top bit of divisor is set.
242 let shift = divisor.limbs[n - 1].leading_zeros();
243 let u = self.shl_bits(shift);
244 let v = divisor.shl_bits(shift);
245
246 let mut u_limbs = u.limbs.clone();
247 // Ensure u has m + n + 1 limbs.
248 while u_limbs.len() <= m + n {
249 u_limbs.push(0);
250 }
251
252 let mut q = vec![0u64; m + 1];
253
254 for j in (0..=m).rev() {
255 // Estimate q_hat.
256 let u_hi = ((u_limbs[j + n] as u128) << 64) | u_limbs[j + n - 1] as u128;
257 let mut q_hat = u_hi / v.limbs[n - 1] as u128;
258 let mut r_hat = u_hi % v.limbs[n - 1] as u128;
259
260 // Refine estimate.
261 loop {
262 if q_hat >= (1u128 << 64)
263 || (q_hat * v.limbs[n - 2] as u128
264 > ((r_hat << 64) | u_limbs[j + n - 2] as u128))
265 {
266 q_hat -= 1;
267 r_hat += v.limbs[n - 1] as u128;
268 if r_hat < (1u128 << 64) {
269 continue;
270 }
271 }
272 break;
273 }
274
275 // Multiply and subtract.
276 let mut borrow: i128 = 0;
277 for i in 0..n {
278 let prod = q_hat * v.limbs[i] as u128;
279 let diff = u_limbs[j + i] as i128 - borrow - (prod as u64) as i128;
280 u_limbs[j + i] = diff as u64;
281 borrow = (prod >> 64) as i128 - (diff >> 64);
282 }
283 let diff = u_limbs[j + n] as i128 - borrow;
284 u_limbs[j + n] = diff as u64;
285
286 q[j] = q_hat as u64;
287
288 // If we subtracted too much, add back.
289 if diff < 0 {
290 q[j] -= 1;
291 let mut carry = 0u64;
292 for i in 0..n {
293 let sum = u_limbs[j + i] as u128 + v.limbs[i] as u128 + carry as u128;
294 u_limbs[j + i] = sum as u64;
295 carry = (sum >> 64) as u64;
296 }
297 u_limbs[j + n] = u_limbs[j + n].wrapping_add(carry);
298 }
299 }
300
301 // Remainder: unnormalize.
302 u_limbs.truncate(n);
303 let rem = Self { limbs: u_limbs }.shr_bits(shift);
304
305 let mut quotient = Self { limbs: q };
306 quotient.normalize();
307
308 (quotient, rem)
309 }
310
311 /// Left shift by `bits` positions (bits < 64).
312 fn shl_bits(&self, bits: u32) -> Self {
313 if bits == 0 {
314 return self.clone();
315 }
316 let mut result = Vec::with_capacity(self.limbs.len() + 1);
317 let mut carry = 0u64;
318 for &limb in &self.limbs {
319 result.push((limb << bits) | carry);
320 carry = limb >> (64 - bits);
321 }
322 if carry > 0 {
323 result.push(carry);
324 }
325 let mut r = Self { limbs: result };
326 r.normalize();
327 r
328 }
329
330 /// Right shift by `bits` positions (bits < 64).
331 fn shr_bits(&self, bits: u32) -> Self {
332 if bits == 0 {
333 return self.clone();
334 }
335 let mut result = Vec::with_capacity(self.limbs.len());
336 let mut carry = 0u64;
337 for &limb in self.limbs.iter().rev() {
338 result.push((limb >> bits) | carry);
339 carry = limb << (64 - bits);
340 }
341 result.reverse();
342 let mut r = Self { limbs: result };
343 r.normalize();
344 r
345 }
346
347 /// Modular exponentiation: self^exp mod modulus.
348 /// Uses left-to-right binary method (square-and-multiply).
349 pub fn modpow(&self, exp: &Self, modulus: &Self) -> Self {
350 assert!(!modulus.is_zero(), "modulus must be non-zero");
351
352 if modulus.limbs == [1] {
353 return Self::zero();
354 }
355
356 // Montgomery multiplication requires an odd modulus.
357 // RSA moduli are always odd, but handle even case with simple modpow.
358 if modulus.limbs[0] & 1 == 0 {
359 return simple_modpow(self, exp, modulus);
360 }
361
362 montgomery_modpow(self, exp, modulus)
363 }
364
365 /// Simple modular reduction: self mod modulus.
366 pub fn modulo(&self, modulus: &Self) -> Self {
367 self.div_rem(modulus).1
368 }
369
370 /// Number of limbs.
371 pub fn num_limbs(&self) -> usize {
372 self.limbs.len()
373 }
374}
375
376impl Ord for BigUint {
377 fn cmp(&self, other: &Self) -> Ordering {
378 let a_len = self.limbs.len();
379 let b_len = other.limbs.len();
380 let max_len = a_len.max(b_len);
381 for i in (0..max_len).rev() {
382 let a = if i < a_len { self.limbs[i] } else { 0 };
383 let b = if i < b_len { other.limbs[i] } else { 0 };
384 match a.cmp(&b) {
385 Ordering::Equal => continue,
386 ord => return ord,
387 }
388 }
389 Ordering::Equal
390 }
391}
392
393impl PartialOrd for BigUint {
394 fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
395 Some(self.cmp(other))
396 }
397}
398
399impl PartialEq for BigUint {
400 fn eq(&self, other: &Self) -> bool {
401 self.cmp(other) == Ordering::Equal
402 }
403}
404
405impl Eq for BigUint {}
406
407// ---------------------------------------------------------------------------
408// Montgomery multiplication for efficient modular exponentiation
409// ---------------------------------------------------------------------------
410
411/// Montgomery context for a given modulus.
412struct Montgomery {
413 /// The modulus N.
414 n: BigUint,
415 /// Number of limbs in N.
416 num_limbs: usize,
417 /// R = 2^(64*num_limbs) (not stored, implicit).
418 /// R^2 mod N — used to convert into Montgomery form.
419 r_squared: BigUint,
420 /// n0_inv = -N^(-1) mod 2^64.
421 n0_inv: u64,
422}
423
424impl Montgomery {
425 fn new(n: &BigUint) -> Self {
426 let num_limbs = n.limbs.len();
427
428 // Compute n0_inv = -N^(-1) mod 2^64 using Newton's method.
429 // Start with x = 1 (since n is odd for RSA moduli).
430 let n0 = n.limbs[0];
431 let mut inv = 1u64;
432 for _ in 0..63 {
433 inv = inv.wrapping_mul(2u64.wrapping_sub(n0.wrapping_mul(inv)));
434 }
435 let n0_inv = inv.wrapping_neg(); // -N^(-1) mod 2^64
436
437 // Compute R^2 mod N where R = 2^(64*num_limbs).
438 // Start with R mod N, then square it.
439 let r_mod_n = {
440 // R = 2^(64*num_limbs). We compute R mod N by creating R and reducing.
441 let mut r_limbs = vec![0u64; num_limbs + 1];
442 r_limbs[num_limbs] = 1;
443 let r = BigUint { limbs: r_limbs };
444 r.div_rem(n).1
445 };
446 let r_squared = r_mod_n.mul(&r_mod_n).div_rem(n).1;
447
448 Self {
449 n: n.clone(),
450 num_limbs,
451 r_squared,
452 n0_inv,
453 }
454 }
455
456 /// Convert a value into Montgomery form: aR mod N.
457 fn to_montgomery(&self, a: &BigUint) -> BigUint {
458 self.montgomery_mul(a, &self.r_squared)
459 }
460
461 /// Convert from Montgomery form back to normal: a * R^(-1) mod N.
462 fn reduce(&self, a: &BigUint) -> BigUint {
463 self.montgomery_mul(a, &BigUint::one())
464 }
465
466 /// Montgomery multiplication: computes a * b * R^(-1) mod N.
467 fn montgomery_mul(&self, a: &BigUint, b: &BigUint) -> BigUint {
468 let n = self.num_limbs;
469 let mut t = vec![0u64; 2 * n + 2];
470
471 for i in 0..n {
472 // t = t + a[i] * b
473 let ai = if i < a.limbs.len() { a.limbs[i] } else { 0 };
474 let mut carry = 0u128;
475 for j in 0..n {
476 let bj = if j < b.limbs.len() { b.limbs[j] } else { 0 };
477 let sum = t[i + j] as u128 + (ai as u128) * (bj as u128) + carry;
478 t[i + j] = sum as u64;
479 carry = sum >> 64;
480 }
481 // Propagate carry.
482 let mut k = i + n;
483 while carry > 0 {
484 let sum = t[k] as u128 + carry;
485 t[k] = sum as u64;
486 carry = sum >> 64;
487 k += 1;
488 }
489
490 // Montgomery reduction step.
491 let m = t[i].wrapping_mul(self.n0_inv);
492 carry = 0u128;
493 for j in 0..n {
494 let sum = t[i + j] as u128 + (m as u128) * (self.n.limbs[j] as u128) + carry;
495 t[i + j] = sum as u64;
496 carry = sum >> 64;
497 }
498 let mut k = i + n;
499 while carry > 0 {
500 let sum = t[k] as u128 + carry;
501 t[k] = sum as u64;
502 carry = sum >> 64;
503 k += 1;
504 }
505 }
506
507 // Result is t[n..2n].
508 let result_limbs: Vec<u64> = t[n..2 * n + 1].to_vec();
509 let mut result = BigUint {
510 limbs: result_limbs,
511 };
512 result.normalize();
513
514 // Final subtraction if result >= N.
515 if result.cmp(&self.n) != Ordering::Less {
516 result = result.sub(&self.n);
517 }
518
519 result
520 }
521}
522
523/// Simple modular exponentiation (for even moduli where Montgomery doesn't work).
524fn simple_modpow(base: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint {
525 let base_reduced = base.modulo(modulus);
526 let mut result = BigUint::one();
527 let bits = exp.bit_len();
528 for i in (0..bits).rev() {
529 result = result.mul(&result).modulo(modulus);
530 if exp.bit(i) {
531 result = result.mul(&base_reduced).modulo(modulus);
532 }
533 }
534 result
535}
536
537/// Modular exponentiation using Montgomery multiplication.
538fn montgomery_modpow(base: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint {
539 let mont = Montgomery::new(modulus);
540
541 // Convert base to Montgomery form.
542 let base_reduced = base.modulo(modulus);
543 let mut result = mont.to_montgomery(&BigUint::one()); // 1 in Montgomery form = R mod N
544 let base_mont = mont.to_montgomery(&base_reduced);
545
546 // Left-to-right binary method.
547 let bits = exp.bit_len();
548 for i in (0..bits).rev() {
549 result = mont.montgomery_mul(&result, &result); // square
550 if exp.bit(i) {
551 result = mont.montgomery_mul(&result, &base_mont); // multiply
552 }
553 }
554
555 mont.reduce(&result)
556}
557
558// ---------------------------------------------------------------------------
559// Tests
560// ---------------------------------------------------------------------------
561
562#[cfg(test)]
563mod tests {
564 use super::*;
565
566 fn hex(bytes: &[u8]) -> String {
567 bytes.iter().map(|b| format!("{b:02x}")).collect()
568 }
569
570 fn from_hex(s: &str) -> Vec<u8> {
571 let s = s.replace(' ', "");
572 (0..s.len())
573 .step_by(2)
574 .map(|i| u8::from_str_radix(&s[i..i + 2], 16).unwrap())
575 .collect()
576 }
577
578 #[test]
579 fn zero_and_one() {
580 assert!(BigUint::zero().is_zero());
581 assert!(!BigUint::one().is_zero());
582 assert_eq!(BigUint::zero().bit_len(), 0);
583 assert_eq!(BigUint::one().bit_len(), 1);
584 }
585
586 #[test]
587 fn from_be_bytes_roundtrip() {
588 let bytes = from_hex("0123456789abcdef");
589 let n = BigUint::from_be_bytes(&bytes);
590 assert_eq!(hex(&n.to_be_bytes()), "0123456789abcdef");
591 }
592
593 #[test]
594 fn from_be_bytes_large() {
595 // 128-bit number.
596 let bytes = from_hex("ffffffffffffffffffffffffffffffff");
597 let n = BigUint::from_be_bytes(&bytes);
598 assert_eq!(hex(&n.to_be_bytes()), "ffffffffffffffffffffffffffffffff");
599 }
600
601 #[test]
602 fn addition() {
603 let a = BigUint::from_be_bytes(&from_hex("ffffffffffffffff"));
604 let b = BigUint::from_u64(1);
605 let c = a.add(&b);
606 assert_eq!(hex(&c.to_be_bytes()), "010000000000000000");
607 }
608
609 #[test]
610 fn subtraction() {
611 let a = BigUint::from_be_bytes(&from_hex("010000000000000000"));
612 let b = BigUint::from_u64(1);
613 let c = a.sub(&b);
614 assert_eq!(hex(&c.to_be_bytes()), "ffffffffffffffff");
615 }
616
617 #[test]
618 fn multiplication() {
619 let a = BigUint::from_u64(0xFFFFFFFF);
620 let b = BigUint::from_u64(0xFFFFFFFF);
621 let c = a.mul(&b);
622 // 0xFFFFFFFF * 0xFFFFFFFF = 0xFFFFFFFE00000001
623 assert_eq!(hex(&c.to_be_bytes()), "fffffffe00000001");
624 }
625
626 #[test]
627 fn division() {
628 let a = BigUint::from_be_bytes(&from_hex("fffffffe00000001"));
629 let b = BigUint::from_u64(0xFFFFFFFF);
630 let (q, r) = a.div_rem(&b);
631 assert_eq!(hex(&q.to_be_bytes()), "ffffffff");
632 assert!(r.is_zero());
633 }
634
635 #[test]
636 fn division_with_remainder() {
637 let a = BigUint::from_u64(17);
638 let b = BigUint::from_u64(5);
639 let (q, r) = a.div_rem(&b);
640 assert_eq!(q, BigUint::from_u64(3));
641 assert_eq!(r, BigUint::from_u64(2));
642 }
643
644 #[test]
645 fn modpow_small() {
646 // 2^10 mod 1000 = 1024 mod 1000 = 24
647 let base = BigUint::from_u64(2);
648 let exp = BigUint::from_u64(10);
649 let modulus = BigUint::from_u64(1000);
650 let result = base.modpow(&exp, &modulus);
651 assert_eq!(result, BigUint::from_u64(24));
652 }
653
654 #[test]
655 fn modpow_fermat() {
656 // Fermat's little theorem: a^(p-1) ≡ 1 (mod p) for prime p, gcd(a,p)=1.
657 // p = 65537 (prime), a = 12345.
658 let a = BigUint::from_u64(12345);
659 let p = BigUint::from_u64(65537);
660 let exp = BigUint::from_u64(65536); // p-1
661 let result = a.modpow(&exp, &p);
662 assert_eq!(result, BigUint::one());
663 }
664
665 #[test]
666 fn modpow_rsa_sized() {
667 // Test with a larger modulus to verify it works at scale.
668 // 3^17 mod 2^128+1 — verify with known computation.
669 let base = BigUint::from_u64(3);
670 let exp = BigUint::from_u64(17);
671 // 2^128 + 1
672 let mut mod_bytes = vec![0u8; 17];
673 mod_bytes[0] = 1;
674 mod_bytes[16] = 1;
675 let modulus = BigUint::from_be_bytes(&mod_bytes);
676 let result = base.modpow(&exp, &modulus);
677 // 3^17 = 129140163. This is < 2^128+1, so result = 129140163.
678 assert_eq!(result, BigUint::from_u64(129140163));
679 }
680
681 #[test]
682 fn to_be_bytes_padded() {
683 let n = BigUint::from_u64(0xFF);
684 let padded = n.to_be_bytes_padded(4);
685 assert_eq!(padded, vec![0, 0, 0, 0xFF]);
686 }
687
688 #[test]
689 fn comparison() {
690 let a = BigUint::from_u64(100);
691 let b = BigUint::from_u64(200);
692 assert_eq!(a.cmp(&b), Ordering::Less);
693 assert_eq!(b.cmp(&a), Ordering::Greater);
694 assert_eq!(a.cmp(&a), Ordering::Equal);
695 }
696
697 #[test]
698 fn multi_limb_division() {
699 // (2^128 - 1) / (2^64 - 1) = 2^64 + 1
700 let a = BigUint::from_be_bytes(&from_hex("ffffffffffffffffffffffffffffffff"));
701 let b = BigUint::from_be_bytes(&from_hex("ffffffffffffffff"));
702 let (q, r) = a.div_rem(&b);
703 assert_eq!(hex(&q.to_be_bytes()), "010000000000000001");
704 assert!(r.is_zero());
705 }
706
707 #[test]
708 fn modpow_256bit() {
709 // Verify modpow with 256-bit numbers.
710 // base = 2, exp = 255, mod = 2^256 - 189 (a prime near 2^256)
711 // 2^255 mod (2^256 - 189)
712 let base = BigUint::from_u64(2);
713 let exp = BigUint::from_u64(255);
714
715 // mod = 2^256 - 189
716 let mut mod_bytes = vec![0xFF; 32];
717 // 2^256 - 189 = 0xFFFFFF...FF43
718 mod_bytes[31] = 0x43; // 0xFF - 188 = 0x43
719
720 let modulus = BigUint::from_be_bytes(&mod_bytes);
721 let result = base.modpow(&exp, &modulus);
722
723 // Verify: result^2 * 2 should equal 2^256 mod (2^256 - 189)
724 // 2^256 mod (2^256 - 189) = 189
725 // So result * 2 mod m should equal... let's just verify result < modulus.
726 assert_eq!(result.cmp(&modulus), Ordering::Less);
727
728 // Double-check: result should be 2^255.
729 // 2^255 < 2^256 - 189, so result = 2^255 exactly.
730 let expected = {
731 let mut b = vec![0u8; 32];
732 b[0] = 0x80;
733 BigUint::from_be_bytes(&b)
734 };
735 assert_eq!(result, expected);
736 }
737}