jaulib v1.3.8
Jau Support Library (C++, Java, ..)
Loading...
Searching...
No Matches
big_int.hpp
Go to the documentation of this file.
1/**
2 * big_int_t (this file)
3 * (c) 2024 Gothel Software e.K.
4 *
5 * Includes code from: Botan BigInt (bigint.h, ..)
6 * (C) 1999-2011,2012,2014,2019 Jack Lloyd
7 * 2007 FlexSecure
8 *
9 * Includes code from: Botan Division Algorithms (divide.h)
10 * (C) 1999-2007,2012,2018,2021 Jack Lloyd
11 *
12 * jaulib including this code is released under the MIT License (see COPYING)
13 * Botan itself is released under the Simplified BSD License (see COPYING)
14 *
15 * SPDX-License-Identifier: MIT
16 *
17 * This Source Code Form is subject to the terms of the MIT License
18 * If a copy of the MIT was not distributed with this file,
19 * you can obtain one at https://opensource.org/license/mit/.
20 */
21
22#ifndef JAU_BIG_INT_HPP_
23#define JAU_BIG_INT_HPP_
24
25#include <cstdint>
26#include <cassert>
27
29#include <jau/byte_util.hpp>
30#include <jau/string_util.hpp>
31
32namespace jau::mp {
33
34 /** \addtogroup Integer
35 *
36 * @{
37 */
38
39 /**
40 * Arbitrary precision integer type
41 *
42 * @anchor bigint_storage_format
43 * ### Local storage format
44 * Internally the big integer is stored in an array of mp_word_t ordered little-endian alike,
45 * with the least significant word at the array-bottom and most significant word at the array-top.
46 *
47 * The mp_word_t itself is stored in jau::endian::native!
48 */
49 class BigInt {
50 public:
51 /**
52 * Sign symbol definitions for positive and negative numbers
53 */
54 enum sign_t { negative = 0, positive = 1 };
55
56 public:
57 BigInt() noexcept = default;
58
59 BigInt(const BigInt& o) noexcept = default;
60
61 ~BigInt() noexcept = default;
62
63 /**
64 * Create a 0-value big_int
65 */
66 static BigInt zero() { return BigInt(); }
67
68 /**
69 * Create a 1-value big_int
70 */
71 static BigInt one() { return BigInt::from_word(1); }
72
73 /**
74 * Create big_int from an unsigned 64 bit integer
75 * @param n initial value of this big_int
76 */
77 static BigInt from_u64(uint64_t n) {
78 BigInt bn;
79 if( 64 == mp_word_bits ) {
80 bn.set_word_at(0, n);
81 } else {
82 bn.set_word_at(1, static_cast<mp_word_t>(n >> 32));
83 bn.set_word_at(0, static_cast<mp_word_t>(n));
84 }
85 return bn;
86 }
87
88 /**
89 * Create big_int from a mp_word_t (limb)
90 * @param n initial value of this big_int
91 */
93 BigInt bn;
94 bn.set_word_at(0, n);
95 return bn;
96 }
97
98 /**
99 * Create big_int from a signed 32 bit integer
100 * @param n initial value of this big_int
101 */
102 static BigInt from_s32(int32_t n) {
103 if(n >= 0) {
104 return BigInt::from_u64(static_cast<uint64_t>(n));
105 } else {
106 return -BigInt::from_u64(static_cast<uint64_t>(-n));
107 }
108 }
109
110 /**
111 * Create big_int of specified size, all zeros
112 * @param n size of the internal register in words
113 */
114 static BigInt with_capacity(size_t n) {
115 BigInt bn;
116 bn.grow_to(n);
117 return bn;
118 }
119
120 /**
121 * Create a power of two
122 * @param n the power of two to create
123 * @return big_int_t representing 2^n
124 */
125 static BigInt power_of_2(size_t n) {
126 BigInt b;
127 b.set_bit(n);
128 return b;
129 }
130
131 /**
132 * Create big_int_t from an unsigned 64 bit integer
133 * @param n initial value of this big_int_t
134 */
135 BigInt(uint64_t n) {
136 if( 64 == mp_word_bits ) {
137 m_data.set_word_at(0, n);
138 } else {
139 m_data.set_word_at(1, static_cast<mp_word_t>(n >> 32));
140 m_data.set_word_at(0, static_cast<mp_word_t>(n));
141 }
142 }
143
144 /**
145 * Construct a big_int_t from a string encoded as hexadecimal or decimal.
146 *
147 * Both number bases may lead a `-`, denoting a negative number.
148 *
149 * Hexadecimal is detected by a leading `0x`.
150 */
151 BigInt(const std::string& str) {
152 size_t markers = 0;
153 bool is_negative = false;
154
155 if(str.length() > 0 && str[0] == '-')
156 {
157 markers += 1;
158 is_negative = true;
159 }
160
161 if(str.length() > markers + 2 && str[markers ] == '0' &&
162 str[markers + 1] == 'x')
163 {
164 markers += 2;
165 *this = hex_decode(cast_char_ptr_to_uint8(str.data()) + markers, str.length() - markers, lb_endian_t::big);
166 } else {
167 *this = dec_decode(cast_char_ptr_to_uint8(str.data()) + markers, str.length() - markers);
168 }
169
171 else set_sign(positive);
172 }
173
174 /**
175 * Create a big_int_t from an integer in a byte array with given byte_len,
176 * considering the given byte_order.
177 *
178 * The value is stored in the local storage format, see \ref bigint_storage_format
179 *
180 * @param buf the byte array holding the value
181 * @param byte_len size of buf in bytes
182 * @param littleEndian
183 */
184 BigInt(const uint8_t buf[], size_t byte_len, const lb_endian_t byte_order) {
185 binary_decode(buf, byte_len, byte_order);
186 }
187
188 BigInt(std::vector<mp_word_t>&& other_reg) noexcept {
189 this->swap_reg(other_reg);
190 }
191
192 BigInt(BigInt&& other) noexcept {
193 this->swap(other);
194 }
195
196 BigInt& operator=(const BigInt& r) = default;
197
198 BigInt& operator=(BigInt&& other) noexcept {
199 if(this != &other) {
200 this->swap(other);
201 }
202 return *this;
203 }
204
205 /**
206 * Swap this value with another
207 * @param other big_int to swap values with
208 */
209 void swap(BigInt& other) noexcept {
210 m_data.swap(other.m_data);
211 std::swap(m_signedness, other.m_signedness);
212 }
213
214 private:
215 void swap_reg(std::vector<mp_word_t>& reg) noexcept {
216 m_data.swap(reg);
217 // sign left unchanged
218 }
219
220 public:
221 /** Unary negation operator, returns new negative instance of this. */
222 BigInt operator-() const noexcept { return BigInt(*this).flip_sign(); }
223
224 /**
225 * @param n the offset to get a byte from
226 * @result byte at offset n
227 */
228 uint8_t byte_at(size_t n) const noexcept {
229 return get_byte_var_be(sizeof(mp_word_t) - (n % sizeof(mp_word_t)) - 1,
230 word_at(n / sizeof(mp_word_t)));
231 }
232
233 /**
234 * Return the mp_word_t at a specified position of the internal register
235 * @param n position in the register
236 * @return value at position n
237 */
238 mp_word_t word_at(size_t n) const noexcept {
239 return m_data.get_word_at(n);
240 }
241
242 /**
243 * Return a const pointer to the register
244 * @result a pointer to the start of the internal register
245 */
246 const mp_word_t* data() const { return m_data.const_data(); }
247
248 /**
249 * Tests if the sign of the integer is negative
250 * @result true, iff the integer has a negative sign
251 */
252 bool is_negative() const noexcept { return sign() == negative; }
253
254 /**
255 * Tests if the sign of the integer is positive
256 * @result true, iff the integer has a positive sign
257 */
258 bool is_positive() const noexcept { return sign() == positive; }
259
260 /**
261 * Return the sign of the integer
262 * @result the sign of the integer
263 */
264 sign_t sign() const noexcept { return m_signedness; }
265
266 /**
267 * @result the opposite sign of the represented integer value
268 */
269 sign_t reverse_sign() const noexcept {
270 if(sign() == positive) {
271 return negative;
272 }
273 return positive;
274 }
275
276 /**
277 * Flip the sign of this big_int
278 */
279 BigInt& flip_sign() noexcept {
280 return set_sign(reverse_sign());
281 }
282
283 /**
284 * Set sign of the integer
285 * @param sign new Sign to set
286 */
288 if(sign == negative && is_zero()) {
289 sign = positive;
290 }
291 m_signedness = sign;
292 return *this;
293 }
294
295 /** Returns absolute (positive) value of this instance */
296 BigInt abs() const noexcept {
297 return BigInt(*this).set_sign(positive);
298 }
299
300 /**
301 * Give size of internal register
302 * @result size of internal register in words
303 */
304 size_t size() const noexcept { return m_data.size(); }
305
306 /**
307 * Return how many words we need to hold this value
308 * @result significant words of the represented integer value
309 */
310 size_t sig_words() const noexcept { return m_data.sig_words(); }
311
312 /** Returns byte length of this integer */
313 size_t bytes() const { return jau::round_up(bits(), 8U) / 8U; }
314
315 /** Returns bit length of this integer */
316 size_t bits() const noexcept {
317 const size_t words = sig_words();
318 if(words == 0) {
319 return 0;
320 }
321 const size_t full_words = (words - 1) * mp_word_bits;
322 const size_t top_bits = mp_word_bits - top_bits_free();
323 return full_words + top_bits;
324 }
325
326 /**
327 * Zeroize the big_int. The size of the underlying register is not
328 * modified.
329 */
330 void clear() { m_data.set_to_zero(); m_signedness = positive; }
331
332 /**
333 * Compares this instance against other considering both sign() value
334 * Returns
335 * - -1 if this < other
336 * - 0 if this == other
337 * - 1 if this > other
338 */
339 int compare(const BigInt& b) const noexcept {
340 return cmp(b, true);
341 }
342
343 /**
344 * Test if the integer has an even value
345 * @result true if the integer is even, false otherwise
346 */
347 bool is_even() const noexcept { return get_bit(0) == 0; }
348
349 /**
350 * Test if the integer has an odd value
351 * @result true if the integer is odd, false otherwise
352 */
353 bool is_odd() const noexcept { return get_bit(0) == 1; }
354
355 /**
356 * Test if the integer is not zero
357 * @result true if the integer is non-zero, false otherwise
358 */
359 bool is_nonzero() const noexcept { return !is_zero(); }
360
361 /**
362 * Test if the integer is zero
363 * @result true if the integer is zero, false otherwise
364 */
365 bool is_zero() const noexcept {
366 return sig_words() == 0;
367 }
368
369 /**
370 * Return bit value at specified position
371 * @param n the bit offset to test
372 * @result true, if the bit at position n is set, false otherwise
373 */
374 bool get_bit(size_t n) const noexcept {
375 return (word_at(n / mp_word_bits) >> (n % mp_word_bits)) & 1;
376 }
377
378 /**
379 * Set bit at specified position
380 * @param n bit position to set
381 */
382 void set_bit(size_t n) noexcept {
383 conditionally_set_bit(n, true);
384 }
385
386 /**
387 * Conditionally set bit at specified position. Note if set_it is
388 * false, nothing happens, and if the bit is already set, it
389 * remains set.
390 *
391 * @param n bit position to set
392 * @param set_it if the bit should be set
393 */
394 void conditionally_set_bit(size_t n, bool set_it) noexcept {
395 const size_t which = n / mp_word_bits;
396 const mp_word_t mask = static_cast<mp_word_t>(set_it) << (n % mp_word_bits);
397 m_data.set_word_at(which, word_at(which) | mask);
398 }
399
400 /** ! operator, returns `true` iff this is zero, otherwise false. */
401 bool operator !() const noexcept { return is_zero(); }
402
403 bool operator==(const BigInt& b) const noexcept { return is_equal(b); }
404 bool operator!=(const BigInt& b) const noexcept { return !is_equal(b); }
405 bool operator<=(const BigInt& b) const noexcept { return cmp(b) <= 0; }
406 bool operator>=(const BigInt& b) const noexcept { return cmp(b) >= 0; }
407 bool operator<(const BigInt& b) const noexcept { return is_less_than(b); }
408 bool operator>(const BigInt& b) const noexcept { return b.is_less_than(*this); }
409
410 /**
411 std::strong_ordering operator<=>(const BigInt& b) const noexcept {
412 const int r = cmp(b);
413 return 0 == r ? std::strong_ordering::equal : ( 0 > r ? std::strong_ordering::less : std::strong_ordering::greater);
414 } */
415
416 BigInt& operator++() noexcept { return *this += 1; }
417
418 BigInt& operator--() noexcept { return *this -= 1; }
419
420 BigInt operator++(int) noexcept { BigInt x = (*this); ++(*this); return x; }
421
422 BigInt operator--(int) noexcept { BigInt x = (*this); --(*this); return x; }
423
424 BigInt& operator+=(const BigInt& y ) noexcept {
425 return add(y.data(), y.sig_words(), y.sign());
426 }
427
428 BigInt& operator-=(const BigInt& y ) noexcept {
429 return add(y.data(), y.sig_words(), y.sign() == positive ? negative : positive);
430 }
431
432 BigInt operator+(const BigInt& y ) const noexcept {
433 return add2(*this, y.data(), y.sig_words(), y.sign());
434 }
435
436 BigInt operator-(const BigInt& y ) const noexcept {
437 return add2(*this, y.data(), y.sig_words(), y.reverse_sign());
438 }
439
440 BigInt& operator<<=(size_t shift) noexcept {
441 const size_t shift_words = shift / mp_word_bits;
442 const size_t shift_bits = shift % mp_word_bits;
443 const size_t size = sig_words();
444
445 const size_t bits_free = top_bits_free();
446
447 const size_t new_size = size + shift_words + (bits_free < shift_bits);
448
449 m_data.grow_to(new_size);
450
451 ops::bigint_shl1(m_data.mutable_data(), new_size, size, shift_words, shift_bits);
452
453 return *this;
454 }
455
456 BigInt& operator>>=(size_t shift) noexcept {
457 const size_t shift_words = shift / mp_word_bits;
458 const size_t shift_bits = shift % mp_word_bits;
459
460 ops::bigint_shr1(m_data.mutable_data(), m_data.size(), shift_words, shift_bits);
461
462 if(is_negative() && is_zero()) {
464 }
465 return *this;
466 }
467
468 BigInt operator<<(size_t shift) const {
469 const size_t shift_words = shift / mp_word_bits;
470 const size_t shift_bits = shift % mp_word_bits;
471 const size_t x_sw = sig_words();
472
473 BigInt y = BigInt::with_capacity(x_sw + shift_words + (shift_bits ? 1 : 0));
474 ops::bigint_shl2(y.mutable_data(), data(), x_sw, shift_words, shift_bits);
475 y.set_sign(sign());
476 return y;
477 }
478
479 BigInt operator>>(size_t shift) const {
480 const size_t shift_words = shift / mp_word_bits;
481 const size_t shift_bits = shift % mp_word_bits;
482 const size_t x_sw = sig_words();
483
484 if(shift_words >= x_sw) {
485 return BigInt::zero();
486 }
487
488 BigInt y = BigInt::with_capacity(x_sw - shift_words);
489 ops::bigint_shr2(y.mutable_data(), data(), x_sw, shift_words, shift_bits);
490
491 if(is_negative() && y.is_zero()) {
493 } else {
494 y.set_sign(sign());
495 }
496 return y;
497 }
498
499
500 BigInt& operator*=(const BigInt& y) noexcept {
501 std::vector<mp_word_t> ws;
502 return this->mul(y, ws);
503 }
504
505 BigInt operator*(const BigInt& y) noexcept
506 {
507 const size_t x_sw = sig_words();
508 const size_t y_sw = y.sig_words();
509
510 BigInt z;
511 z.resize(size() + y.size());
512
513 if(x_sw == 1 && y_sw) {
514 ops::bigint_linmul3(z.mutable_data(), y.data(), y_sw, word_at(0));
515 } else if(y_sw == 1 && x_sw) {
516 ops::bigint_linmul3(z.mutable_data(), data(), x_sw, y.word_at(0));
517 } else if(x_sw && y_sw) {
518 ops::basecase_mul(z.mutable_data(), z.size(), data(), x_sw, y.data(), y_sw);
519 }
520 z.cond_flip_sign(x_sw > 0 && y_sw > 0 && sign() != y.sign());
521 return z;
522 }
523
525 if(y.sig_words() == 1 && jau::is_power_of_2(y.word_at(0))) {
526 (*this) >>= (y.bits() - 1);
527 } else {
528 (*this) = (*this) / y;
529 }
530 return (*this);
531 }
532
533 BigInt operator/(const BigInt& y) const {
534 if(y.sig_words() == 1) {
535 return *this / y.word_at(0);
536 }
537 BigInt q, r;
538 vartime_divide(*this, y, q, r);
539 return q;
540 }
541 /**
542 * Modulo operator
543 * @param y the modulus to reduce this by
544 */
545 BigInt& operator%=(const BigInt& mod) {
546 return (*this = (*this) % mod);
547 }
548
549 BigInt operator%(const BigInt& mod) {
550 if(mod.is_zero()) {
552 }
553 if(mod.is_negative()) {
554 throw jau::math::MathDomainError("mod < 0", E_FILE_LINE);
555 }
556 if(is_positive() && mod.is_positive() && *this < mod) {
557 return *this;
558 }
559 if(mod.sig_words() == 1) {
560 return from_word(*this % mod.word_at(0));
561 }
562 BigInt q, r;
563 vartime_divide(*this, mod, q, r);
564 return r;
565 }
566
567 /**
568 * Returns (*this)^e, or pow(*this, e)
569 *
570 * Implementation is not optimized and naive, i.e. O(n)
571 *
572 * @param e the exponent
573 */
575 const BigInt& b = *this;
576 if( b.is_zero() ) {
577 return BigInt::zero();
578 }
579 const BigInt one_v = BigInt::one();
580 BigInt r = one_v;
581 bool is_negative;
582 if( e.is_negative() ) {
583 is_negative = true;
584 e.flip_sign();
585 } else {
586 is_negative = false;
587 }
588
589 while( e.is_nonzero() ) {
590 r *= b;
591 --e;
592 }
593
594 if( is_negative ) {
595 return one_v / r;
596 } else {
597 return r;
598 }
599 }
600
601 /**
602 * Returns (*this)^e % m, or pow(*this, e) % m
603 *
604 * Implementation is not optimized and naive, i.e. O(n)
605 *
606 * @param e the exponent
607 */
609 const BigInt& b = *this;
610 if( b.is_zero() ) {
611 return BigInt::zero();
612 }
613 const BigInt one_v = BigInt::one();
614 BigInt r = one_v;
615 bool is_negative;
616 if( e.is_negative() ) {
617 is_negative = true;
618 e.flip_sign();
619 } else {
620 is_negative = false;
621 }
622
623 while( e.is_nonzero() ) {
624 r *= b;
625 r %= m;
626 --e;
627 }
628
629 if( is_negative ) {
630 return one_v / r;
631 } else {
632 return r;
633 }
634 }
635
636 /**
637 * Square value of *this
638 * @param ws a temp workspace
639 */
640 BigInt& square(std::vector<mp_word_t>& ws); // TODO
641
642 /**
643 * Set *this to y - *this
644 * @param y the big_int_t to subtract from
645 * @param ws a temp workspace
646 */
647 BigInt& rev_sub(const BigInt& y, std::vector<mp_word_t>& ws); // TODO
648
649 /**
650 * Set *this to (*this + y) % mod
651 * This function assumes *this is >= 0 && < mod
652 * @param y the big_int_t to add - assumed y >= 0 and y < mod
653 * @param mod the positive modulus
654 * @param ws a temp workspace
655 */
656 BigInt& mod_add(const BigInt& y, const BigInt& mod, std::vector<mp_word_t>& ws); // TODO
657
658 /**
659 * Set *this to (*this - y) % mod
660 * This function assumes *this is >= 0 && < mod
661 * @param y the big_int_t to subtract - assumed y >= 0 and y < mod
662 * @param mod the positive modulus
663 * @param ws a temp workspace
664 */
665 BigInt& mod_sub(const BigInt& y, const BigInt& mod, std::vector<mp_word_t>& ws); // TODO
666
667 /**
668 * Set *this to (*this * y) % mod
669 * This function assumes *this is >= 0 && < mod
670 * y should be small, less than 16
671 * @param y the small integer to multiply by
672 * @param mod the positive modulus
673 * @param ws a temp workspace
674 */
675 BigInt& mod_mul(uint8_t y, const BigInt& mod, std::vector<mp_word_t>& ws); // TODO
676
677 /**
678 * @param rng a random number generator
679 * @param min the minimum value (must be non-negative)
680 * @param max the maximum value (must be non-negative and > min)
681 * @return random integer in [min,max)
682 static big_int_t random_integer(RandomNumberGenerator& rng,
683 const big_int_t& min,
684 const big_int_t& max); // TODO
685 */
686
687 std::string to_dec_string(bool add_details=false) const {
688 // Use the largest power of 10 that fits in a mp_word_t
689 mp_word_t conversion_radix, radix_digits;
690 if constexpr ( 64 == mp_word_bits ) {
691 conversion_radix = 10000000000000000000U;
692 radix_digits = 19;
693 } else {
694 conversion_radix = 1000000000U;
695 radix_digits = 9;
696 }
697
698 // (over-)estimate of the number of digits needed; log2(10) ~ 3.3219
699 const size_t digit_estimate = static_cast<size_t>(1 + (static_cast<double>(this->bits()) / 3.32));
700
701 // (over-)estimate of db such that conversion_radix^db > *this
702 const size_t digit_blocks = (digit_estimate + radix_digits - 1) / radix_digits;
703
704 BigInt value = *this;
705 value.set_sign(positive);
706
707 // Extract groups of digits into words
708 std::vector<mp_word_t> digit_groups(digit_blocks);
709
710 for(size_t i = 0; i != digit_blocks; ++i) {
711 mp_word_t remainder = 0;
712 ct_divide_word(value, conversion_radix, value, remainder);
713 digit_groups[i] = remainder;
714 }
715 assert(value.is_zero());
716
717 // Extract digits from the groups
718 std::vector<uint8_t> digits(digit_blocks * radix_digits);
719
720 for(size_t i = 0; i != digit_blocks; ++i) {
721 mp_word_t remainder = digit_groups[i];
722 for(size_t j = 0; j != radix_digits; ++j)
723 {
724 // Compiler should convert div/mod by 10 into mul by magic constant
725 const mp_word_t digit = remainder % 10;
726 remainder /= 10;
727 digits[radix_digits*i + j] = static_cast<uint8_t>(digit);
728 }
729 }
730
731 // remove leading zeros
732 while(!digits.empty() && digits.back() == 0) {
733 digits.pop_back();
734 }
735
736 assert(digit_estimate >= digits.size());
737
738 // Reverse the digits to big-endian and format to text
739 std::string s;
740 s.reserve(1 + digits.size());
741
742 if(is_negative()) {
743 s += "-";
744 }
745
746 // Reverse and convert to textual digits
747 for(uint8_t d : digits) {
748 s.push_back(static_cast<char>(d + '0')); // assumes ASCII
749 }
750
751 if(s.empty()) {
752 s += "0";
753 }
754 if( add_details ) {
755 append_detail(s);
756 }
757 return s;
758 }
759
760 std::string to_hex_string(bool add_details=false) const noexcept {
761 std::vector<uint8_t> bits;
762 const uint8_t* data;
763 size_t data_len;
764
765 if( is_zero() ) {
766 bits.push_back(0);
767 data = bits.data();
768 data_len = bits.size();
769 } else {
770 data = reinterpret_cast<const uint8_t*>(m_data.const_data());
771 data_len = bytes();
772 }
773
774 std::string s;
775 if(is_negative()) {
776 s += "-";
777 }
778 s.append( jau::bytesHexString(data, data_len,
779 false /* lsbFirst */, true /* lowerCase */) );
780 if( add_details ) {
781 append_detail(s);
782 }
783 return s;
784 }
785
786 private:
787 class data_t
788 {
789 public:
790 data_t() noexcept
791 : m_reg(), m_sig_words(sig_words_npos) {}
792
793 data_t(const data_t& o) noexcept = default;
794
795 data_t(data_t&& o) noexcept {
796 swap(o);
797 }
798 data_t(std::vector<mp_word_t>&& reg) noexcept {
799 swap(reg);
800 }
801
802 ~data_t() noexcept = default;
803
804 data_t& operator=(const data_t& r) noexcept = default;
805
806 data_t& operator=(data_t&& other) noexcept {
807 if(this != &other) {
808 this->swap(other);
809 }
810 return *this;
811 }
812 data_t& operator=(std::vector<mp_word_t>&& other_reg) noexcept {
813 if(&m_reg != &other_reg) {
814 this->swap(other_reg);
815 }
816 return *this;
817 }
818
819 mp_word_t* mutable_data() noexcept {
820 invalidate_sig_words();
821 return m_reg.data();
822 }
823
824 const mp_word_t* const_data() const noexcept {
825 return m_reg.data();
826 }
827
828 std::vector<mp_word_t>& mutable_vector() noexcept {
829 invalidate_sig_words();
830 return m_reg;
831 }
832
833 const std::vector<mp_word_t>& const_vector() const noexcept {
834 return m_reg;
835 }
836
837 mp_word_t get_word_at(size_t n) const noexcept {
838 if(n < m_reg.size()) {
839 return m_reg[n];
840 }
841 return 0;
842 }
843
844 void set_word_at(size_t i, mp_word_t w) {
845 invalidate_sig_words();
846 if(i >= m_reg.size()) {
847 if(w == 0) {
848 return;
849 }
850 grow_to(i + 1);
851 }
852 m_reg[i] = w;
853 }
854
855 void set_words(const mp_word_t w[], size_t len) {
856 invalidate_sig_words();
857 m_reg.assign(w, w + len);
858 }
859
860 void set_to_zero() {
861 m_reg.resize(m_reg.capacity());
862 clear_mem(m_reg.data(), m_reg.size());
863 m_sig_words = 0;
864 }
865
866 void set_size(size_t s) {
867 invalidate_sig_words();
868 clear_mem(m_reg.data(), m_reg.size());
869 m_reg.resize(s + (8 - (s % 8)));
870 }
871
872 void mask_bits(size_t n) noexcept {
873 if(n == 0) { return set_to_zero(); }
874
875 const size_t top_word = n / mp_word_bits;
876
877 // if(top_word < sig_words()) ?
878 if(top_word < size())
879 {
880 const mp_word_t mask = (static_cast<mp_word_t>(1) << (n % mp_word_bits)) - 1;
881 const size_t len = size() - (top_word + 1);
882 if(len > 0)
883 {
884 clear_mem(&m_reg[top_word+1], len);
885 }
886 m_reg[top_word] &= mask;
887 invalidate_sig_words();
888 }
889 }
890
891 void grow_to(size_t n) const {
892 if(n > size()) {
893 if(n <= m_reg.capacity()) {
894 m_reg.resize(n);
895 } else {
896 m_reg.resize(n + (8 - (n % 8)));
897 }
898 }
899 }
900
901 size_t size() const noexcept { return m_reg.size(); }
902
903 void shrink_to_fit(size_t min_size = 0) {
904 const size_t words = std::max(min_size, sig_words());
905 m_reg.resize(words);
906 }
907
908 void resize(size_t s) {
909 m_reg.resize(s);
910 }
911
912 void swap(data_t& other) noexcept {
913 m_reg.swap(other.m_reg);
914 std::swap(m_sig_words, other.m_sig_words);
915 }
916
917 void swap(std::vector<mp_word_t>& reg) noexcept {
918 m_reg.swap(reg);
919 invalidate_sig_words();
920 }
921
922 void invalidate_sig_words() const noexcept {
923 m_sig_words = sig_words_npos;
924 }
925
926 size_t sig_words() const noexcept {
927 if(m_sig_words == sig_words_npos) {
928 m_sig_words = calc_sig_words();
929 } else {
930 assert(m_sig_words == calc_sig_words());
931 }
932 return m_sig_words;
933 }
934 private:
935 static const size_t sig_words_npos = static_cast<size_t>(-1);
936
937 size_t calc_sig_words() const noexcept {
938 const size_t sz = m_reg.size();
939 size_t sig = sz;
940
941 mp_word_t sub = 1;
942
943 for(size_t i = 0; i != sz; ++i)
944 {
945 const mp_word_t w = m_reg[sz - i - 1];
946 sub &= ct_is_zero(w);
947 sig -= sub;
948 }
949
950 /*
951 * This depends on the data so is poisoned, but unpoison it here as
952 * later conditionals are made on the size.
953 */
955
956 return sig;
957 }
958
959 mutable std::vector<mp_word_t> m_reg;
960 mutable size_t m_sig_words = sig_words_npos;
961 };
962 data_t m_data;
963 sign_t m_signedness = positive;
964
965 /**
966 * Byte extraction of big-endian value
967 * @param byte_num which byte to extract, 0 == highest byte
968 * @param input the value to extract from
969 * @return byte byte_num of input
970 */
971 template<typename T> static constexpr uint8_t get_byte_var_be(size_t byte_num, T input) noexcept {
972 return static_cast<uint8_t>( input >> (((~byte_num)&(sizeof(T)-1)) << 3) );
973 }
974 /**
975 * Byte extraction of little-endian value
976 * @param byte_num which byte to extract, 0 == lowest byte
977 * @param input the value to extract from
978 * @return byte byte_num of input
979 */
980 template<typename T> static constexpr uint8_t get_byte_var_le(size_t byte_num, T input) noexcept {
981 return static_cast<uint8_t>( input >> ( byte_num << 3 ) );
982 }
983
984 /**
985 * Zero out some bytes. Warning: use secure_scrub_memory instead if the
986 * memory is about to be freed or otherwise the compiler thinks it can
987 * elide the writes.
988 *
989 * @param ptr a pointer to memory to zero
990 * @param bytes the number of bytes to zero in ptr
991 */
992 static constexpr void clear_bytes(void* ptr, size_t bytes) noexcept {
993 if(bytes > 0) {
994 std::memset(ptr, 0, bytes);
995 }
996 }
997
998 /**
999 * Zero memory before use. This simply calls memset and should not be
1000 * used in cases where the compiler cannot see the call as a
1001 * side-effecting operation (for example, if calling clear_mem before
1002 * deallocating memory, the compiler would be allowed to omit the call
1003 * to memset entirely under the as-if rule.)
1004 *
1005 * @param ptr a pointer to an array of Ts to zero
1006 * @param n the number of Ts pointed to by ptr
1007 */
1008 template<typename T> static constexpr void clear_mem(T* ptr, size_t n) noexcept {
1009 clear_bytes(ptr, sizeof(T)*n);
1010 }
1011
1012 /**
1013 * Increase internal register buffer to at least n words
1014 * @param n new size of register
1015 */
1016 void grow_to(size_t n) const { m_data.grow_to(n); }
1017
1018 void resize(size_t s) { m_data.resize(s); }
1019
1020 void set_word_at(size_t i, mp_word_t w) {
1021 m_data.set_word_at(i, w);
1022 }
1023
1024 void set_words(const mp_word_t w[], size_t len) {
1025 m_data.set_words(w, len);
1026 }
1027 /**
1028 * Return a mutable pointer to the register
1029 * @result a pointer to the start of the internal register
1030 */
1031 mp_word_t* mutable_data() { return m_data.mutable_data(); }
1032
1033 /**
1034 * Set this number to the value in buf with given byte_len,
1035 * considering the given byte_order.
1036 *
1037 * The value is stored in the local storage format, see \ref bigint_storage_format
1038 */
1039 void binary_decode(const uint8_t buf[], size_t byte_len, const lb_endian_t byte_order) {
1040 const size_t full_words = byte_len / sizeof(mp_word_t);
1041 const size_t extra_bytes = byte_len % sizeof(mp_word_t);
1042
1043 // clear() + setting size
1044 m_signedness = positive;
1045 m_data.set_size( jau::round_up(full_words + (extra_bytes > 0 ? 1U : 0U), 8U) );
1046
1047 mp_word_t* sink = m_data.mutable_data();
1048 if( is_little_endian(byte_order) ) {
1049 // little-endian to local (words arranged as little-endian w/ word itself in native-endian)
1050 for(size_t i = 0; i < full_words; ++i) {
1051 sink[i] = jau::get_value<mp_word_t>(buf + sizeof(mp_word_t)*i, byte_order);
1052 }
1053 } else {
1054 // big-endian to local (words arranged as little-endian w/ word itself in native-endian)
1055 for(size_t i = 0; i < full_words; ++i) {
1056 sink[i] = jau::get_value<mp_word_t>(buf + byte_len - sizeof(mp_word_t)*(i+1), byte_order);
1057 }
1058 }
1059 mp_word_t le_w = 0;
1060 if( is_little_endian(byte_order) ) {
1061 for(size_t i = 0; i < extra_bytes; ++i) {
1062 le_w |= mp_word_t( buf[full_words*sizeof(mp_word_t) + i] ) << ( i * 8 ); // next lowest byte
1063 }
1064 } else {
1065 for(size_t i = 0; i < extra_bytes; ++i) {
1066 le_w = (le_w << 8) | mp_word_t( buf[i] ); // buf[0] highest byte
1067 }
1068 }
1069 sink[full_words] = jau::le_to_cpu( le_w );
1070 }
1071
1072 /**
1073 * Set this number to the decoded hex-string value of buf with given str_len,
1074 * considering the given byte_order.
1075 *
1076 * The value is stored in the local storage format, see \ref bigint_storage_format
1077 */
1078 static BigInt hex_decode(const uint8_t buf[], size_t str_len, const lb_endian_t byte_order) {
1079 BigInt r;
1080
1081 std::vector<uint8_t> bin_out;
1082 const size_t exp_blen = str_len / 2 + str_len % 2;
1083 const size_t blen = jau::hexStringBytes(bin_out, buf, str_len, is_little_endian(byte_order), false /* checkLeading0x */);
1084 if( exp_blen != blen ) {
1085 throw jau::math::MathDomainError("invalid hexadecimal char @ "+std::to_string(blen)+"/"+std::to_string(exp_blen)+" of '"+
1086 std::string(cast_uint8_ptr_to_char(buf), str_len)+"'", E_FILE_LINE);
1087 }
1088 r.binary_decode(bin_out.data(), bin_out.size(), lb_endian_t::little);
1089 return r;
1090 }
1091
1092 static BigInt dec_decode(const uint8_t buf[], size_t str_len) {
1093 BigInt r;
1094
1095 // This could be made faster using the same trick as to_dec_string
1096 for(size_t i = 0; i < str_len; ++i) {
1097 const char c = static_cast<char>(buf[i]);
1098
1099 if(c < '0' || c > '9') {
1100 throw jau::math::MathDomainError("invalid decimal char", E_FILE_LINE);
1101 }
1102 const uint8_t x = c - '0';
1103 assert(x < 10);
1104
1105 r *= 10;
1106 r += x;
1107 }
1108 return r;
1109 }
1110
1111 public:
1112 /**
1113 * Stores this number to the value in buf with given byte_len,
1114 * considering the given byte_order.
1115 *
1116 * The value is read from the local storage in its format, see \ref bigint_storage_format
1117 *
1118 * If byte_len is less than the byt-esize of this integer, i.e. bytes(),
1119 * then it will be truncated.
1120 *
1121 * If byte_len is greater than the byte-size of this integer, i.e. bytes(), it will be zero-padded.
1122 *
1123 * @return actual number of bytes copied, i.e. min(byte_len, bytes());
1124 */
1125 size_t binary_encode(uint8_t output[], size_t byte_len, const lb_endian_t byte_order) const noexcept {
1126 const size_t full_words = byte_len / sizeof(mp_word_t);
1127 const size_t extra_bytes = byte_len % sizeof(mp_word_t);
1128
1129 if( is_little_endian( byte_order ) ) {
1130 // to little-endian from local (words arranged as little-endian w/ word itself in native-endian)
1131 for(size_t i = 0; i < full_words; ++i) {
1132 jau::put_value<mp_word_t>(output + i*sizeof(mp_word_t), word_at(i), byte_order);
1133 }
1134 } else {
1135 // to big-endian from local (words arranged as little-endian w/ word itself in native-endian)
1136 for(size_t i = 0; i < full_words; ++i) {
1137 jau::put_value<mp_word_t>(output + byte_len - (i+1)*sizeof(mp_word_t), word_at(i), byte_order);
1138 }
1139 }
1140 if(extra_bytes > 0) {
1141 const mp_word_t le_w = jau::cpu_to_le( word_at(full_words) );
1142 if( is_little_endian( byte_order ) ) {
1143 for(size_t i = 0; i < extra_bytes; ++i) {
1144 output[full_words*sizeof(mp_word_t) + i] = get_byte_var_le(i, le_w); // next lowest byte
1145 }
1146 } else {
1147 for(size_t i = 0; i < extra_bytes; ++i) {
1148 output[extra_bytes-1-i] = get_byte_var_le(i, le_w); // output[0] highest byte
1149 }
1150 }
1151 }
1152 return extra_bytes + full_words * sizeof(mp_word_t);
1153 }
1154
1155 private:
1156 size_t top_bits_free() const noexcept {
1157 const size_t words = sig_words();
1158
1159 const mp_word_t top_word = word_at(words - 1);
1160 const size_t bits_used = jau::high_bit(top_word);
1161 CT::unpoison(bits_used);
1162 return mp_word_bits - bits_used;
1163 }
1164
1165 /**
1166 * Compares this instance against other, considering sign depending on check_signs
1167 * Returns
1168 * - -1 if this < other
1169 * - 0 if this == other
1170 * - 1 if this > other
1171 */
1172 int cmp(const BigInt& other, bool check_signs = true) const noexcept {
1173 if(check_signs) {
1174 if(other.is_positive() && this->is_negative()) {
1175 return -1;
1176 }
1177 if(other.is_negative() && this->is_positive()) {
1178 return 1;
1179 }
1180 if(other.is_negative() && this->is_negative()) {
1181 return (-ops::bigint_cmp(this->data(), this->size(),
1182 other.data(), other.size()));
1183 }
1184 }
1185 return ops::bigint_cmp(this->data(), this->size(),
1186 other.data(), other.size());
1187 }
1188
1189 bool is_equal(const BigInt& other) const noexcept {
1190 if(this->sign() != other.sign()) {
1191 return false;
1192 }
1193 return ops::bigint_ct_is_eq(this->data(), this->sig_words(),
1194 other.data(), other.sig_words()).is_set();
1195 }
1196
1197 bool is_less_than(const BigInt& other) const noexcept {
1198 if(this->is_negative() && other.is_positive()) {
1199 return true;
1200 }
1201 if(this->is_positive() && other.is_negative()) {
1202 return false;
1203 }
1204 if(other.is_negative() && this->is_negative()) {
1205 return ops::bigint_ct_is_lt(other.data(), other.sig_words(),
1206 this->data(), this->sig_words()).is_set();
1207 }
1208 return ops::bigint_ct_is_lt(this->data(), this->sig_words(),
1209 other.data(), other.sig_words()).is_set();
1210 }
1211
1212 BigInt& add(const mp_word_t y[], size_t y_words, sign_t y_sign) {
1213 const size_t x_sw = sig_words();
1214 grow_to(std::max(x_sw, y_words) + 1);
1215
1216 if(sign() == y_sign)
1217 {
1218 ops::bigint_add2(mutable_data(), size() - 1, y, y_words);
1219 } else {
1220 const int32_t relative_size = ops::bigint_cmp(data(), x_sw, y, y_words);
1221
1222 if(relative_size >= 0)
1223 {
1224 // *this >= y
1225 ops::bigint_sub2(mutable_data(), x_sw, y, y_words);
1226 } else {
1227 // *this < y
1228 ops::bigint_sub2_rev(mutable_data(), y, y_words);
1229 }
1230 //this->sign_fixup(relative_size, y_sign);
1231 if(relative_size < 0) {
1232 set_sign(y_sign);
1233 } else if(relative_size == 0) {
1235 }
1236 }
1237 return (*this);
1238 }
1239
1240 BigInt& operator+=(mp_word_t y) noexcept {
1241 return add(&y, 1, sign_t::positive);
1242 }
1243 BigInt& operator-=(mp_word_t y) noexcept {
1244 return add(&y, 1, sign_t::negative);
1245 }
1246
1247 static BigInt add2(const BigInt& x, const mp_word_t y[], size_t y_words, sign_t y_sign) {
1248 const size_t x_sw = x.sig_words();
1249
1250 BigInt z = BigInt::with_capacity(std::max(x_sw, y_words) + 1);
1251
1252 if(x.sign() == y_sign) {
1253 ops::bigint_add3(z.mutable_data(), x.data(), x_sw, y, y_words);
1254 z.set_sign(x.sign());
1255 } else {
1256 const int32_t relative_size = ops::bigint_sub_abs(z.mutable_data(), x.data(), x_sw, y, y_words);
1257
1258 //z.sign_fixup(relative_size, y_sign);
1259 if(relative_size < 0) {
1260 z.set_sign(y_sign);
1261 } else if(relative_size == 0) {
1262 z.set_sign(positive);
1263 } else {
1264 z.set_sign(x.sign());
1265 }
1266 }
1267 return z;
1268 }
1269
1270 BigInt& mul(const BigInt& y, std::vector<mp_word_t>& ws) noexcept {
1271 const size_t x_sw = sig_words();
1272 const size_t y_sw = y.sig_words();
1273 set_sign((sign() == y.sign()) ? positive : negative);
1274
1275 if(x_sw == 0 || y_sw == 0)
1276 {
1277 clear();
1279 }
1280 else if(x_sw == 1 && y_sw)
1281 {
1282 grow_to(y_sw + 1);
1283 ops::bigint_linmul3(mutable_data(), y.data(), y_sw, word_at(0));
1284 }
1285 else if(y_sw == 1 && x_sw)
1286 {
1287 mp_word_t carry = ops::bigint_linmul2(mutable_data(), x_sw, y.word_at(0));
1288 set_word_at(x_sw, carry);
1289 }
1290 else
1291 {
1292 const size_t new_size = x_sw + y_sw + 1;
1293 // ws.resize(new_size);
1294 (void)ws;
1295 std::vector<mp_word_t> z_reg(new_size);
1296
1297 ops::basecase_mul(z_reg.data(), z_reg.size(), data(), x_sw, y.data(), y_sw);
1298
1299 this->swap_reg(z_reg);
1300 }
1301
1302 return (*this);
1303 }
1304
1306 const size_t x_sw = sig_words();
1307 BigInt z = BigInt::with_capacity(x_sw + 1);
1308
1309 if(x_sw && y) {
1310 ops::bigint_linmul3(z.mutable_data(), data(), x_sw, y);
1311 z.set_sign(sign());
1312 }
1313 return z;
1314 }
1315
1316 void cond_flip_sign(bool predicate) noexcept {
1317 // This code is assuming Negative == 0, Positive == 1
1318 const auto mask = CT::Mask<uint8_t>::expand(predicate);
1319 const uint8_t current_sign = static_cast<uint8_t>(sign());
1320 const uint8_t new_sign = mask.select(current_sign ^ 1, current_sign);
1321 set_sign(static_cast<sign_t>(new_sign));
1322 }
1323
1324 /**
1325 * Return *this % mod
1326 *
1327 * Assumes that *this is (if anything) only slightly larger than
1328 * mod and performs repeated subtractions. It should not be used if
1329 * *this is much larger than mod, instead use modulo operator.
1330 */
1331 inline size_t reduce_below(const BigInt& p, std::vector<mp_word_t>& ws) {
1332 if(p.is_negative() || this->is_negative()) {
1333 std::string msg;
1334 if( p.is_negative() ) {
1335 msg.append("p < 0");
1336 }
1337 if( this->is_negative() ) {
1338 if( msg.length() > 0 ) {
1339 msg.append(" and ");
1340 }
1341 msg.append("*this < 0");
1342 }
1343 throw jau::math::MathDomainError(msg, E_FILE_LINE);
1344 }
1345 const size_t p_words = p.sig_words();
1346
1347 if(size() < p_words + 1) {
1348 grow_to(p_words + 1);
1349 }
1350 if(ws.size() < p_words + 1) {
1351 ws.resize(p_words + 1);
1352 }
1353 clear_mem(ws.data(), ws.size());
1354
1355 size_t reductions = 0;
1356
1357 for(;;)
1358 {
1359 mp_word_t borrow = ops::bigint_sub3(ws.data(), data(), p_words + 1, p.data(), p_words);
1360 if(borrow) {
1361 break;
1362 }
1363 ++reductions;
1364 swap_reg(ws);
1365 }
1366 return reductions;
1367 }
1368
1369 static void sign_fixup(const BigInt& x, const BigInt& y, BigInt& q, BigInt& r) {
1370 q.cond_flip_sign(x.sign() != y.sign());
1371
1372 if(x.is_negative() && r.is_nonzero())
1373 {
1374 q -= 1;
1375 r = y.abs() - r;
1376 }
1377 }
1378
1379 static bool division_check(mp_word_t q, mp_word_t y2, mp_word_t y1,
1380 mp_word_t x3, mp_word_t x2, mp_word_t x1) noexcept {
1381 /*
1382 Compute (y3,y2,y1) = (y2,y1) * q
1383 and return true if (y3,y2,y1) > (x3,x2,x1)
1384 */
1385
1386 mp_word_t y3 = 0;
1387 y1 = ops::word_madd2(q, y1, y3);
1388 y2 = ops::word_madd2(q, y2, y3);
1389
1390 const mp_word_t x[3] = { x1, x2, x3 };
1391 const mp_word_t y[3] = { y1, y2, y3 };
1392
1393 return ops::bigint_ct_is_lt(x, 3, y, 3).is_set();
1394 }
1395
1396 /*
1397 * Solve x = q * y + r
1398 *
1399 * See Handbook of Applied Cryptography section 14.2.5
1400 */
1401 static void vartime_divide(const BigInt& x, const BigInt& y_arg, BigInt& q_out, BigInt& r_out) {
1402 if( y_arg.is_zero() ) {
1403 throw jau::math::MathDivByZeroError("y_arg == 0", E_FILE_LINE);
1404 }
1405 const size_t y_words = y_arg.sig_words();
1406
1407 assert(y_words > 0);
1408
1409 BigInt y = y_arg;
1410
1411 BigInt r = x;
1412 BigInt q = BigInt::zero();
1413 std::vector<mp_word_t> ws;
1414
1415 r.set_sign(BigInt::positive);
1416 y.set_sign(BigInt::positive);
1417
1418 // Calculate shifts needed to normalize y with high bit set
1419 const size_t shifts = y.top_bits_free();
1420
1421 y <<= shifts;
1422 r <<= shifts;
1423
1424 // we know y has not changed size, since we only shifted up to set high bit
1425 const size_t t = y_words - 1;
1426 const size_t n = std::max(y_words, r.sig_words()) - 1; // r may have changed size however
1427 assert(n >= t);
1428
1429 q.grow_to(n - t + 1);
1430
1431 mp_word_t* q_words = q.mutable_data();
1432
1433 BigInt shifted_y = y << (mp_word_bits * (n-t));
1434
1435 // Set q_{n-t} to number of times r > shifted_y
1436 q_words[n-t] = r.reduce_below(shifted_y, ws);
1437
1438 const mp_word_t y_t0 = y.word_at(t);
1439 const mp_word_t y_t1 = y.word_at(t-1);
1440 assert((y_t0 >> (mp_word_bits-1)) == 1);
1441
1442 for(size_t j = n; j != t; --j) {
1443 const mp_word_t x_j0 = r.word_at(j);
1444 const mp_word_t x_j1 = r.word_at(j-1);
1445 const mp_word_t x_j2 = r.word_at(j-2);
1446
1447 mp_word_t qjt = ops::bigint_divop(x_j0, x_j1, y_t0);
1448
1449 qjt = CT::Mask<mp_word_t>::is_equal(x_j0, y_t0).select(mp_word_max, qjt);
1450
1451 // Per HAC 14.23, this operation is required at most twice
1452 qjt -= division_check(qjt, y_t0, y_t1, x_j0, x_j1, x_j2);
1453 qjt -= division_check(qjt, y_t0, y_t1, x_j0, x_j1, x_j2);
1454 assert(division_check(qjt, y_t0, y_t1, x_j0, x_j1, x_j2) == false);
1455
1456 shifted_y >>= mp_word_bits;
1457 // Now shifted_y == y << (BOTAN_MP_WORD_BITS * (j-t-1))
1458
1459 // TODO this sequence could be better
1460 r -= shifted_y * qjt;
1461 qjt -= r.is_negative();
1462 r += shifted_y * static_cast<mp_word_t>(r.is_negative());
1463
1464 q_words[j-t-1] = qjt;
1465 }
1466
1467 r >>= shifts;
1468
1469 sign_fixup(x, y_arg, q, r);
1470
1471 r_out = r;
1472 q_out = q;
1473 }
1474
1475 BigInt operator/(const mp_word_t& y) const {
1476 if(y == 0) {
1477 throw jau::math::MathDivByZeroError("y == 0", E_FILE_LINE);
1478 }
1479 BigInt q;
1480 mp_word_t r;
1481 ct_divide_word(*this, y, q, r);
1482 return q;
1483 }
1484
1485 static void ct_divide_word(const BigInt& x, mp_word_t y, BigInt& q_out, mp_word_t& r_out) {
1486 if(y == 0) {
1487 throw jau::math::MathDivByZeroError("y == 0", E_FILE_LINE);
1488 }
1489 const size_t x_words = x.sig_words();
1490 const size_t x_bits = x.bits();
1491
1492 BigInt q = BigInt::with_capacity(x_words);
1493 mp_word_t r = 0;
1494
1495 for(size_t i = 0; i != x_bits; ++i)
1496 {
1497 const size_t b = x_bits - 1 - i;
1498 const bool x_b = x.get_bit(b);
1499
1500 const auto r_carry = CT::Mask<mp_word_t>::expand(r >> (mp_word_bits - 1));
1501
1502 r *= 2;
1503 r += x_b;
1504
1505 const auto r_gte_y = CT::Mask<mp_word_t>::is_gte(r, y) | r_carry;
1506 q.conditionally_set_bit(b, r_gte_y.is_set());
1507 r = r_gte_y.select(r - y, r);
1508 }
1509
1510 if(x.is_negative()) {
1511 q.flip_sign();
1512 if(r != 0) {
1513 --q;
1514 r = y - r;
1515 }
1516 }
1517
1518 r_out = r;
1519 q_out = q;
1520 }
1521
1523 if(mod == 0) {
1524 throw jau::math::MathDivByZeroError("mod == 0", E_FILE_LINE);
1525 }
1526 if(mod == 1) {
1527 return 0;
1528 }
1529 mp_word_t remainder = 0;
1530
1531 if( jau::is_power_of_2(mod) ) {
1532 remainder = (word_at(0) & (mod - 1));
1533 } else {
1534 const size_t sw = sig_words();
1535 for(size_t i = sw; i > 0; --i) {
1536 remainder = ops::bigint_modop(remainder, word_at(i-1), mod);
1537 }
1538 }
1539 if(remainder && sign() == BigInt::negative) {
1540 return mod - remainder;
1541 }
1542 return remainder;
1543 }
1544
1545 void append_detail(std::string& s) const noexcept {
1546 s.append(", bits ").append(std::to_string(bits())).append(", ").append(std::to_string(sig_words())).append(" word(s): ");
1547 for(size_t i=0; i<sig_words(); ++i) {
1548 const mp_word_t w = word_at(i);
1549 s.append( jau::bytesHexString(&w, mp_word_bits/CHAR_BIT, false /* lsbFirst */, true /* lowerCase */) )
1550 .append(", ");
1551 }
1552 }
1553 };
1554
1555 /**@}*/
1556}
1557
1558namespace jau {
1559
1560 /** \addtogroup Integer
1561 *
1562 * @{
1563 */
1564
1565 inline mp::BigInt abs(const mp::BigInt& x) noexcept { return x.abs(); }
1566 inline mp::BigInt pow(mp::BigInt b, const mp::BigInt& e) { return b.pow(e); }
1567
1568 inline const mp::BigInt& min(const mp::BigInt& x, const mp::BigInt& y) noexcept {
1569 return x < y ? x : y;
1570 }
1571 inline const mp::BigInt& max(const mp::BigInt& x, const mp::BigInt& y) noexcept {
1572 return x > y ? x : y;
1573 }
1574 inline const mp::BigInt& clamp(const mp::BigInt& x, const mp::BigInt& min_val, const mp::BigInt& max_val) noexcept {
1575 return min(max(x, min_val), max_val);
1576 }
1577
1578 inline mp::BigInt& min(mp::BigInt& x, mp::BigInt& y) noexcept {
1579 return x < y ? x : y;
1580 }
1581 inline mp::BigInt& max(mp::BigInt& x, mp::BigInt& y) noexcept {
1582 return x > y ? x : y;
1583 }
1584 inline mp::BigInt& clamp(mp::BigInt& x, mp::BigInt& min_val, mp::BigInt& max_val) noexcept {
1585 return min(max(x, min_val), max_val);
1586 }
1587
1588 inline mp::BigInt gcd(const mp::BigInt& a, const mp::BigInt& b) noexcept {
1589 mp::BigInt a_ = abs(a);
1590 mp::BigInt b_ = abs(b);
1591 while( b_.is_nonzero() ) {
1592 const mp::BigInt t = b_;
1593 b_ = a_ % b_;
1594 a_ = t;
1595 }
1596 return a_;
1597 }
1598
1599 /**@}*/
1600}
1601
1602namespace std {
1603 inline std::ostream& operator<<(std::ostream& out, const jau::mp::BigInt& v) {
1604 return out << v.to_dec_string();
1605 }
1606}
1607
1608#endif /** JAU_BIG_INT_HPP_ */
#define E_FILE_LINE
static Mask< T > is_gte(T x, T y) noexcept
Return a Mask<T> which is set if x >= y.
Definition ct_utils.hpp:177
static Mask< T > is_equal(T x, T y) noexcept
Return a Mask<T> which is set if x == y.
Definition ct_utils.hpp:145
static Mask< T > expand(T v) noexcept
Return a Mask<T> which is set if v is != 0.
Definition ct_utils.hpp:119
math_error_t::div_by_zero, i.e.
math_error_t::invalid
Arbitrary precision integer type.
Definition big_int.hpp:49
BigInt & set_sign(sign_t sign) noexcept
Set sign of the integer.
Definition big_int.hpp:287
BigInt(uint64_t n)
Create big_int_t from an unsigned 64 bit integer.
Definition big_int.hpp:135
bool is_even() const noexcept
Test if the integer has an even value.
Definition big_int.hpp:347
size_t bits() const noexcept
Returns bit length of this integer.
Definition big_int.hpp:316
std::string to_dec_string(bool add_details=false) const
Definition big_int.hpp:687
bool operator==(const BigInt &b) const noexcept
Definition big_int.hpp:403
bool is_odd() const noexcept
Test if the integer has an odd value.
Definition big_int.hpp:353
bool get_bit(size_t n) const noexcept
Return bit value at specified position.
Definition big_int.hpp:374
BigInt & operator-=(const BigInt &y) noexcept
Definition big_int.hpp:428
BigInt abs() const noexcept
Returns absolute (positive) value of this instance.
Definition big_int.hpp:296
BigInt & operator++() noexcept
std::strong_ordering operator<=>(const BigInt& b) const noexcept { const int r = cmp(b); return 0 == ...
Definition big_int.hpp:416
BigInt operator-() const noexcept
Unary negation operator, returns new negative instance of this.
Definition big_int.hpp:222
BigInt & mod_add(const BigInt &y, const BigInt &mod, std::vector< mp_word_t > &ws)
Set *this to (*this + y) % mod This function assumes *this is >= 0 && < mod.
BigInt pow(BigInt e)
Returns (*this)^e, or pow(*this, e)
Definition big_int.hpp:574
BigInt & operator--() noexcept
Definition big_int.hpp:418
size_t bytes() const
Returns byte length of this integer.
Definition big_int.hpp:313
BigInt & operator<<=(size_t shift) noexcept
Definition big_int.hpp:440
BigInt(const uint8_t buf[], size_t byte_len, const lb_endian_t byte_order)
Create a big_int_t from an integer in a byte array with given byte_len, considering the given byte_or...
Definition big_int.hpp:184
BigInt & operator=(const BigInt &r)=default
BigInt operator>>(size_t shift) const
Definition big_int.hpp:479
std::string to_hex_string(bool add_details=false) const noexcept
Definition big_int.hpp:760
BigInt operator--(int) noexcept
Definition big_int.hpp:422
sign_t reverse_sign() const noexcept
Definition big_int.hpp:269
bool operator<=(const BigInt &b) const noexcept
Definition big_int.hpp:405
BigInt operator+(const BigInt &y) const noexcept
Definition big_int.hpp:432
sign_t
Sign symbol definitions for positive and negative numbers.
Definition big_int.hpp:54
BigInt operator%(const BigInt &mod)
Definition big_int.hpp:549
bool operator>(const BigInt &b) const noexcept
Definition big_int.hpp:408
const mp_word_t * data() const
Return a const pointer to the register.
Definition big_int.hpp:246
bool is_nonzero() const noexcept
Test if the integer is not zero.
Definition big_int.hpp:359
size_t size() const noexcept
Give size of internal register.
Definition big_int.hpp:304
static BigInt one()
Create a 1-value big_int.
Definition big_int.hpp:71
sign_t sign() const noexcept
Return the sign of the integer.
Definition big_int.hpp:264
bool is_negative() const noexcept
Tests if the sign of the integer is negative.
Definition big_int.hpp:252
BigInt & operator>>=(size_t shift) noexcept
Definition big_int.hpp:456
BigInt & rev_sub(const BigInt &y, std::vector< mp_word_t > &ws)
Set *this to y - *this.
BigInt & operator=(BigInt &&other) noexcept
Definition big_int.hpp:198
static BigInt zero()
Create a 0-value big_int.
Definition big_int.hpp:66
void conditionally_set_bit(size_t n, bool set_it) noexcept
Conditionally set bit at specified position.
Definition big_int.hpp:394
size_t sig_words() const noexcept
Return how many words we need to hold this value.
Definition big_int.hpp:310
void clear()
Zeroize the big_int.
Definition big_int.hpp:330
BigInt(std::vector< mp_word_t > &&other_reg) noexcept
Definition big_int.hpp:188
BigInt() noexcept=default
BigInt & operator%=(const BigInt &mod)
Modulo operator.
Definition big_int.hpp:545
bool operator!() const noexcept
!
Definition big_int.hpp:401
BigInt & mod_sub(const BigInt &y, const BigInt &mod, std::vector< mp_word_t > &ws)
Set *this to (*this - y) % mod This function assumes *this is >= 0 && < mod.
bool is_positive() const noexcept
Tests if the sign of the integer is positive.
Definition big_int.hpp:258
bool operator<(const BigInt &b) const noexcept
Definition big_int.hpp:407
bool operator!=(const BigInt &b) const noexcept
Definition big_int.hpp:404
bool operator>=(const BigInt &b) const noexcept
Definition big_int.hpp:406
size_t binary_encode(uint8_t output[], size_t byte_len, const lb_endian_t byte_order) const noexcept
Stores this number to the value in buf with given byte_len, considering the given byte_order.
Definition big_int.hpp:1125
void swap(BigInt &other) noexcept
Swap this value with another.
Definition big_int.hpp:209
int compare(const BigInt &b) const noexcept
Compares this instance against other considering both sign() value Returns.
Definition big_int.hpp:339
static BigInt from_s32(int32_t n)
Create big_int from a signed 32 bit integer.
Definition big_int.hpp:102
BigInt & square(std::vector< mp_word_t > &ws)
Square value of *this.
BigInt & operator/=(const BigInt &y)
Definition big_int.hpp:524
static BigInt from_word(mp_word_t n)
Create big_int from a mp_word_t (limb)
Definition big_int.hpp:92
static BigInt with_capacity(size_t n)
Create big_int of specified size, all zeros.
Definition big_int.hpp:114
bool is_zero() const noexcept
Test if the integer is zero.
Definition big_int.hpp:365
mp_word_t word_at(size_t n) const noexcept
Return the mp_word_t at a specified position of the internal register.
Definition big_int.hpp:238
BigInt operator/(const BigInt &y) const
Definition big_int.hpp:533
BigInt operator-(const BigInt &y) const noexcept
Definition big_int.hpp:436
BigInt & operator*=(const BigInt &y) noexcept
Definition big_int.hpp:500
BigInt & mod_mul(uint8_t y, const BigInt &mod, std::vector< mp_word_t > &ws)
Set *this to (*this * y) % mod This function assumes *this is >= 0 && < mod y should be small,...
uint8_t byte_at(size_t n) const noexcept
Definition big_int.hpp:228
BigInt mod_pow(BigInt e, const BigInt &m)
Returns (*this)^e % m, or pow(*this, e) % m.
Definition big_int.hpp:608
BigInt & flip_sign() noexcept
Flip the sign of this big_int.
Definition big_int.hpp:279
BigInt & operator+=(const BigInt &y) noexcept
Definition big_int.hpp:424
void set_bit(size_t n) noexcept
Set bit at specified position.
Definition big_int.hpp:382
BigInt operator*(const BigInt &y) noexcept
Definition big_int.hpp:505
BigInt operator<<(size_t shift) const
Definition big_int.hpp:468
BigInt(const std::string &str)
Construct a big_int_t from a string encoded as hexadecimal or decimal.
Definition big_int.hpp:151
static BigInt from_u64(uint64_t n)
Create big_int from an unsigned 64 bit integer.
Definition big_int.hpp:77
BigInt operator++(int) noexcept
Definition big_int.hpp:420
BigInt(BigInt &&other) noexcept
Definition big_int.hpp:192
static BigInt power_of_2(size_t n)
Create a power of two.
Definition big_int.hpp:125
constexpr bool is_little_endian() noexcept
Evaluates true if platform is running in little endian mode, i.e.
constexpr std::enable_if_t< std::is_standard_layout_v< T >, T > get_value(uint8_t const *buffer) noexcept
Returns a T value from the given byte address using packed_t to resolve a potential memory alignment ...
char * cast_uint8_ptr_to_char(uint8_t *b) noexcept
constexpr uint16_t le_to_cpu(uint16_t const l) noexcept
lb_endian_t
Simplified reduced endian type only covering little- and big-endian.
constexpr std::enable_if_t< std::is_standard_layout_v< T >, void > put_value(uint8_t *buffer, const T &v) noexcept
Put the given T value into the given byte address using packed_t to resolve a potential memory alignm...
constexpr uint16_t cpu_to_le(uint16_t const h) noexcept
const uint8_t * cast_char_ptr_to_uint8(const char *s) noexcept
@ little
Identifier for little endian, equivalent to endian::little.
@ big
Identifier for big endian, equivalent to endian::big.
constexpr T ct_is_zero(T x)
Returns ~0 (2-complement) if arg is zero, otherwise 0 (w/o branching) in O(1) and constant time (CT).
constexpr bool value(const Bool rhs) noexcept
void swap(cow_darray< Value_type, Size_type, Alloc_type > &rhs, cow_darray< Value_type, Size_type, Alloc_type > &lhs) noexcept
constexpr T clamp(const T x, const T min_val, const T max_val) noexcept
Returns constrained integral value to lie between given min- and maximum value (w/ branching) in O(1)...
constexpr T round_up(const T n, const U align_to)
Round up w/ branching in O(1)
Definition int_math.hpp:99
constexpr nsize_t high_bit(T x)
Return the index of the highest set bit w/ branching (loop) in O(n), actually O(n/2).
Definition int_math.hpp:162
constexpr T min(const T x, const T y) noexcept
Returns the minimum of two integrals (w/ branching) in O(1)
constexpr T gcd(T a, T b) noexcept
Returns the greatest common divisor (GCD) of the two given integer values following Euclid's algorith...
Definition int_math.hpp:330
constexpr T max(const T x, const T y) noexcept
Returns the maximum of two integrals (w/ branching) in O(1)
mp::BigInt pow(mp::BigInt b, const mp::BigInt &e)
Definition big_int.hpp:1566
constexpr T abs(const T x) noexcept
Returns the absolute value of an arithmetic number (w/ branching) in O(1)
constexpr bool is_power_of_2(const T x) noexcept
Power of 2 test (w/o branching ?) in O(1)
Definition int_math.hpp:138
std::string bytesHexString(const void *data, const nsize_t length, const bool lsbFirst, const bool lowerCase=true, const bool skipLeading0x=false) noexcept
Produce a hexadecimal string representation of the given byte values.
size_t hexStringBytes(std::vector< uint8_t > &out, const std::string &hexstr, const bool lsbFirst, const bool checkLeading0x) noexcept
Converts a given hexadecimal string representation into a byte vector.
void unpoison(const T *p, size_t n)
Definition ct_utils.hpp:54
CT::Mask< mp_word_t > bigint_ct_is_lt(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size, bool lt_or_equal=false) noexcept
Compare x and y Return ~0 if x[0:x_size] < y[0:y_size] or 0 otherwise If lt_or_equal is true,...
CT::Mask< mp_word_t > bigint_ct_is_eq(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
int bigint_cmp(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Compare unsigned x and y mp_word_t.
mp_word_t bigint_divop(mp_word_t n1, mp_word_t n0, mp_word_t d)
Computes ((n1<<bits) + n0) / d.
int32_t bigint_sub_abs(mp_word_t z[], const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Set z to abs(x-y), ie if x >= y, then compute z = x - y Otherwise compute z = y - x No borrow is poss...
mp_word_t bigint_add2(mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Two operand addition with carry out.
void bigint_shr1(mp_word_t x[], size_t x_size, size_t word_shift, size_t bit_shift) noexcept
mp_word_t bigint_linmul2(mp_word_t x[], size_t x_size, mp_word_t y) noexcept
void bigint_linmul3(mp_word_t z[], const mp_word_t x[], size_t x_size, mp_word_t y) noexcept
mp_word_t bigint_sub3(mp_word_t z[], const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Three operand subtraction.
void bigint_shl2(mp_word_t y[], const mp_word_t x[], size_t x_size, size_t word_shift, size_t bit_shift) noexcept
void bigint_add3(mp_word_t z[], const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Three operand addition.
mp_word_t bigint_sub2(mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
Two operand subtraction.
mp_word_t bigint_modop(mp_word_t n1, mp_word_t n0, mp_word_t d)
Compute ((n1<<bits) + n0) % d.
mp_word_t word_madd2(mp_word_t a, mp_word_t b, mp_word_t &c) noexcept
Word Multiply/Add.
void bigint_shl1(mp_word_t x[], size_t x_size, size_t x_words, size_t word_shift, size_t bit_shift) noexcept
void basecase_mul(mp_word_t z[], size_t z_size, const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
void bigint_shr2(mp_word_t y[], const mp_word_t x[], size_t x_size, size_t word_shift, size_t bit_shift) noexcept
void bigint_sub2_rev(mp_word_t x[], const mp_word_t y[], size_t y_size) noexcept
Two operand subtraction, x = y - x; assumes y >= x.
big_int_t (this file) (c) 2024 Gothel Software e.K.
Definition big_int.hpp:32
jau::uint_bytes< impl::best_word_byte_size()>::type mp_word_t
constexpr const size_t mp_word_bits
constexpr const mp_word_t mp_word_max
__pack(...): Produces MSVC, clang and gcc compatible lead-in and -out macros.
Definition backtrace.hpp:32
STL namespace.
std::ostream & operator<<(std::ostream &os, const jau::fraction_timespec &v) noexcept
Output stream operator for jau::fraction_timespec.