jaulib v1.3.8
Jau Support Library (C++, Java, ..)
Loading...
Searching...
No Matches
big_int_ops.hpp
Go to the documentation of this file.
1/*
2 * big_int_ops.hpp (this file)
3 * (c) 2024 Gothel Software e.K.
4 *
5 * Includes code from: Botan Lowest Level MPI Algorithms (mp_asmi.h)
6 * (C) 1999-2010 Jack Lloyd
7 * 2006 Luca Piccarreta
8 *
9 * Includes code from: Botan MPI Algorithms (mp_core.h)
10 * (C) 1999-2010,2018 Jack Lloyd
11 * 2006 Luca Piccarreta
12 * 2016 Matthias Gierlings
13 *
14 * jaulib including this code is released under the MIT License (see COPYING)
15 * Botan itself is released under the Simplified BSD License (see COPYING)
16 *
17 * SPDX-License-Identifier: MIT
18 *
19 * This Source Code Form is subject to the terms of the MIT License
20 * If a copy of the MIT was not distributed with this file,
21 * you can obtain one at https://opensource.org/license/mit/.
22 */
23
24#ifndef JAU_BIG_INT_OPS_HPP_
25#define JAU_BIG_INT_OPS_HPP_
26
27#include <cstdint>
28#include <cassert>
29
30#include <jau/cpp_lang_util.hpp>
31#include <jau/cpuid.hpp>
32#include <jau/ct_utils.hpp>
34
35namespace jau::mp {
36 namespace impl {
37 constexpr size_t best_word_byte_size() {
38 if constexpr ( 64 == jau::cpu::pointer_bit_size() && is_builtin_int128_available() ) {
39 return 8;
40 } else {
41 return 4;
42 }
43 }
44 }
45 #undef JAU_FORCE_MP_WORD_32_BITS
46 #if !defined( JAU_FORCE_MP_WORD_32_BITS )
47 constexpr const size_t mp_word_bits = impl::best_word_byte_size() * CHAR_BIT;
50 constexpr const bool has_mp_dword = is_builtin_int128_available();
51 #elif 1
52 constexpr const size_t mp_word_bits = 32;
53 typedef uint32_t mp_word_t;
54 constexpr const bool has_mp_dword = true;
55 typedef uint64_t mp_dword_t;
56 #elif 0
57 constexpr const size_t mp_word_bits = 64;
58 typedef uint64_t mp_word_t;
59 #if defined(__SIZEOF_INT128__)
60 constexpr const bool has_mp_dword = true;
61 typedef uint128_t mp_dword_t;
62 #else
63 constexpr const bool has_mp_dword = false;
64 #endif
65 #endif
66 constexpr const mp_word_t mp_word_max = ~static_cast<mp_word_t>(0);
67}
68
69namespace jau::mp::ops {
70
71 int bigint_cmp(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept;
72
73 inline mp_word_t word_add(const mp_word_t x, const mp_word_t y, mp_word_t& carry) noexcept {
74 mp_word_t z = x + y;
75 mp_word_t c1 = (z < x);
76 z += carry;
77 carry = c1 | (z < carry);
78 return z;
79 }
80
81 inline mp_word_t word8_add2(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept {
82 x[0] = word_add(x[0], y[0], carry);
83 x[1] = word_add(x[1], y[1], carry);
84 x[2] = word_add(x[2], y[2], carry);
85 x[3] = word_add(x[3], y[3], carry);
86 x[4] = word_add(x[4], y[4], carry);
87 x[5] = word_add(x[5], y[5], carry);
88 x[6] = word_add(x[6], y[6], carry);
89 x[7] = word_add(x[7], y[7], carry);
90 return carry;
91 }
92
93 inline mp_word_t word8_add3(mp_word_t z[8], const mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept
94 {
95 z[0] = word_add(x[0], y[0], carry);
96 z[1] = word_add(x[1], y[1], carry);
97 z[2] = word_add(x[2], y[2], carry);
98 z[3] = word_add(x[3], y[3], carry);
99 z[4] = word_add(x[4], y[4], carry);
100 z[5] = word_add(x[5], y[5], carry);
101 z[6] = word_add(x[6], y[6], carry);
102 z[7] = word_add(x[7], y[7], carry);
103 return carry;
104 }
105
106 inline mp_word_t word_sub(mp_word_t x, mp_word_t y, mp_word_t& carry) noexcept {
107 mp_word_t t0 = x - y;
108 mp_word_t c1 = (t0 > x);
109 mp_word_t z = t0 - carry;
110 carry = c1 | (z > t0);
111 return z;
112 }
113
114 inline mp_word_t word8_sub2(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept {
115 x[0] = word_sub(x[0], y[0], carry);
116 x[1] = word_sub(x[1], y[1], carry);
117 x[2] = word_sub(x[2], y[2], carry);
118 x[3] = word_sub(x[3], y[3], carry);
119 x[4] = word_sub(x[4], y[4], carry);
120 x[5] = word_sub(x[5], y[5], carry);
121 x[6] = word_sub(x[6], y[6], carry);
122 x[7] = word_sub(x[7], y[7], carry);
123 return carry;
124 }
125
126 inline mp_word_t word8_sub2_rev(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept {
127 x[0] = word_sub(y[0], x[0], carry);
128 x[1] = word_sub(y[1], x[1], carry);
129 x[2] = word_sub(y[2], x[2], carry);
130 x[3] = word_sub(y[3], x[3], carry);
131 x[4] = word_sub(y[4], x[4], carry);
132 x[5] = word_sub(y[5], x[5], carry);
133 x[6] = word_sub(y[6], x[6], carry);
134 x[7] = word_sub(y[7], x[7], carry);
135 return carry;
136 }
137
138 inline mp_word_t word8_sub3(mp_word_t z[8], const mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept {
139 z[0] = word_sub(x[0], y[0], carry);
140 z[1] = word_sub(x[1], y[1], carry);
141 z[2] = word_sub(x[2], y[2], carry);
142 z[3] = word_sub(x[3], y[3], carry);
143 z[4] = word_sub(x[4], y[4], carry);
144 z[5] = word_sub(x[5], y[5], carry);
145 z[6] = word_sub(x[6], y[6], carry);
146 z[7] = word_sub(x[7], y[7], carry);
147 return carry;
148 }
149
150 /** 64x64->128 bit multiplication */
151 inline void mul64x64_128(const uint64_t a, const uint64_t b, uint64_t& lo, uint64_t& hi) noexcept {
152 if constexpr ( is_builtin_int128_available() ) {
153 const uint128_t r = static_cast<uint128_t>(a) * b;
154 hi = (r >> 64) & 0xFFFFFFFFFFFFFFFFUL;
155 lo = (r ) & 0xFFFFFFFFFFFFFFFFUL;
156 } else {
157 /*
158 * Do a 64x64->128 multiply using four 32x32->64 multiplies plus
159 * some adds and shifts. Last resort for CPUs like UltraSPARC (with
160 * 64-bit registers/ALU, but no 64x64->128 multiply) or 32-bit CPUs.
161 */
162 constexpr const size_t HWORD_BITS = 32;
163 constexpr const uint32_t HWORD_MASK = 0xFFFFFFFFU;
164
165 const uint32_t a_hi = (a >> HWORD_BITS);
166 const uint32_t a_lo = (a & HWORD_MASK);
167 const uint32_t b_hi = (b >> HWORD_BITS);
168 const uint32_t b_lo = (b & HWORD_MASK);
169
170 uint64_t x0 = static_cast<uint64_t>(a_hi) * b_hi;
171 uint64_t x1 = static_cast<uint64_t>(a_lo) * b_hi;
172 uint64_t x2 = static_cast<uint64_t>(a_hi) * b_lo;
173 uint64_t x3 = static_cast<uint64_t>(a_lo) * b_lo;
174
175 // this cannot overflow as (2^32-1)^2 + 2^32-1 < 2^64-1
176 x2 += x3 >> HWORD_BITS;
177
178 // this one can overflow
179 x2 += x1;
180
181 // propagate the carry if any
182 x0 += static_cast<uint64_t>(static_cast<bool>(x2 < x1)) << HWORD_BITS;
183
184 hi = x0 + (x2 >> HWORD_BITS);
185 lo = ((x2 & HWORD_MASK) << HWORD_BITS) + (x3 & HWORD_MASK);
186 }
187 }
188
189 /** Word Multiply/Add */
191 if constexpr ( has_mp_dword ) {
192 const mp_dword_t s = static_cast<mp_dword_t>(a) * b + c;
193 c = static_cast<mp_word_t>(s >> mp_word_bits);
194 return static_cast<mp_word_t>(s);
195 } else if constexpr ( 64 == mp_word_bits ) {
196 uint64_t hi, lo;
197
198 mul64x64_128(a, b, lo, hi);
199
200 lo += c;
201 hi += (lo < c); // carry?
202
203 c = hi;
204 return lo;
205 } else {
206 static_assert( has_mp_dword || 64 == mp_word_bits );
207 return 0;
208 }
209 }
210
211 /** Word Multiply/Add */
213 if constexpr ( has_mp_dword ) {
214 const mp_dword_t s = static_cast<mp_dword_t>(a) * b + c + d;
215 d = static_cast<mp_word_t>(s >> mp_word_bits);
216 return static_cast<mp_word_t>(s);
217 } else if constexpr ( 64 == mp_word_bits ) {
218 uint64_t hi, lo;
219
220 mul64x64_128(a, b, lo, hi);
221
222 lo += c;
223 hi += (lo < c); // carry?
224
225 lo += d;
226 hi += (lo < d); // carry?
227
228 d = hi;
229 return lo;
230 } else {
231 static_assert( has_mp_dword || 64 == mp_word_bits );
232 return 0;
233 }
234 }
235
236 /*
237 * Eight Word Block Multiply/Add
238 */
239 inline mp_word_t word8_madd3(mp_word_t z[8], const mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept {
240 z[0] = word_madd3(x[0], y, z[0], carry);
241 z[1] = word_madd3(x[1], y, z[1], carry);
242 z[2] = word_madd3(x[2], y, z[2], carry);
243 z[3] = word_madd3(x[3], y, z[3], carry);
244 z[4] = word_madd3(x[4], y, z[4], carry);
245 z[5] = word_madd3(x[5], y, z[5], carry);
246 z[6] = word_madd3(x[6], y, z[6], carry);
247 z[7] = word_madd3(x[7], y, z[7], carry);
248 return carry;
249 }
250
251 /*
252 * Eight Word Block Linear Multiplication
253 */
254 inline mp_word_t word8_linmul2(mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept {
255 x[0] = word_madd2(x[0], y, carry);
256 x[1] = word_madd2(x[1], y, carry);
257 x[2] = word_madd2(x[2], y, carry);
258 x[3] = word_madd2(x[3], y, carry);
259 x[4] = word_madd2(x[4], y, carry);
260 x[5] = word_madd2(x[5], y, carry);
261 x[6] = word_madd2(x[6], y, carry);
262 x[7] = word_madd2(x[7], y, carry);
263 return carry;
264 }
265
266 /*
267 * Eight Word Block Linear Multiplication
268 */
269 inline mp_word_t word8_linmul3(mp_word_t z[8], const mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept {
270 z[0] = word_madd2(x[0], y, carry);
271 z[1] = word_madd2(x[1], y, carry);
272 z[2] = word_madd2(x[2], y, carry);
273 z[3] = word_madd2(x[3], y, carry);
274 z[4] = word_madd2(x[4], y, carry);
275 z[5] = word_madd2(x[5], y, carry);
276 z[6] = word_madd2(x[6], y, carry);
277 z[7] = word_madd2(x[7], y, carry);
278 return carry;
279 }
280
281 /** Two operand addition with carry out */
282 inline mp_word_t bigint_add2(mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept {
283 mp_word_t carry = 0;
284
285 assert(x_size >= y_size);
286
287 const size_t blocks = y_size - (y_size % 8);
288
289 for(size_t i = 0; i != blocks; i += 8) {
290 carry = word8_add2(x + i, y + i, carry);
291 }
292 for(size_t i = blocks; i != y_size; ++i) {
293 x[i] = word_add(x[i], y[i], carry);
294 }
295 for(size_t i = y_size; i != x_size; ++i) {
296 x[i] = word_add(x[i], 0, carry);
297 }
298 return carry;
299 }
300
301 /** Three operand addition with carry out */
302 inline mp_word_t bigint_add3_nc(mp_word_t z[], const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept {
303 if(x_size < y_size)
304 { return bigint_add3_nc(z, y, y_size, x, x_size); }
305
306 mp_word_t carry = 0;
307
308 const size_t blocks = y_size - (y_size % 8);
309
310 for(size_t i = 0; i != blocks; i += 8) {
311 carry = word8_add3(z + i, x + i, y + i, carry);
312 }
313 for(size_t i = blocks; i != y_size; ++i) {
314 z[i] = word_add(x[i], y[i], carry);
315 }
316 for(size_t i = y_size; i != x_size; ++i) {
317 z[i] = word_add(x[i], 0, carry);
318 }
319 return carry;
320 }
321 /** Three operand addition */
322 inline 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 {
323 z[x_size > y_size ? x_size : y_size] +=
324 bigint_add3_nc(z, x, x_size, y, y_size);
325 }
326
327 /** Two operand subtraction */
328 inline mp_word_t bigint_sub2(mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept {
329 mp_word_t borrow = 0;
330
331 assert(x_size >= y_size);
332
333 const size_t blocks = y_size - (y_size % 8);
334
335 for(size_t i = 0; i != blocks; i += 8) {
336 borrow = word8_sub2(x + i, y + i, borrow);
337 }
338 for(size_t i = blocks; i != y_size; ++i) {
339 x[i] = word_sub(x[i], y[i], borrow);
340 }
341 for(size_t i = y_size; i != x_size; ++i) {
342 x[i] = word_sub(x[i], 0, borrow);
343 }
344 return borrow;
345 }
346
347 /** Two operand subtraction, x = y - x; assumes y >= x */
348 inline void bigint_sub2_rev(mp_word_t x[], const mp_word_t y[], size_t y_size) noexcept {
349 mp_word_t borrow = 0;
350
351 const size_t blocks = y_size - (y_size % 8);
352
353 for(size_t i = 0; i != blocks; i += 8) {
354 borrow = word8_sub2_rev(x + i, y + i, borrow);
355 }
356 for(size_t i = blocks; i != y_size; ++i) {
357 x[i] = word_sub(y[i], x[i], borrow);
358 }
359 assert(borrow == 0); // y must be greater than x
360 }
361
362 /** Three operand subtraction */
363 inline 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 {
364 mp_word_t borrow = 0;
365
366 assert(x_size >= y_size); // Expected sizes
367
368 const size_t blocks = y_size - (y_size % 8);
369
370 for(size_t i = 0; i != blocks; i += 8) {
371 borrow = word8_sub3(z + i, x + i, y + i, borrow);
372 }
373 for(size_t i = blocks; i != y_size; ++i) {
374 z[i] = word_sub(x[i], y[i], borrow);
375 }
376 for(size_t i = y_size; i != x_size; ++i) {
377 z[i] = word_sub(x[i], 0, borrow);
378 }
379 return borrow;
380 }
381
382 /**
383 * Set z to abs(x-y), ie if x >= y, then compute z = x - y
384 * Otherwise compute z = y - x
385 * No borrow is possible since the result is always >= 0
386 *
387 * Return the relative size of x vs y (-1, 0, 1)
388 *
389 * @param z output array of max(x_size,y_size) words
390 * @param x input param
391 * @param x_size length of x
392 * @param y input param
393 * @param y_size length of y
394 */
395 inline int32_t
396 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 {
397 const int32_t relative_size = bigint_cmp(x, x_size, y, y_size);
398
399 // Swap if relative_size == -1
400 const bool need_swap = relative_size < 0;
401 CT::conditional_swap_ptr(need_swap, x, y);
402 CT::conditional_swap(need_swap, x_size, y_size);
403
404 /*
405 * We know at this point that x >= y so if y_size is larger than
406 * x_size, we are guaranteed they are just leading zeros which can
407 * be ignored
408 */
409 y_size = std::min(x_size, y_size);
410
411 bigint_sub3(z, x, x_size, y, y_size);
412
413 return relative_size;
414 }
415
416 /*
417 * Linear Multiply - returns the carry
418 */
419 [[nodiscard]] inline mp_word_t bigint_linmul2(mp_word_t x[], size_t x_size, mp_word_t y) noexcept {
420 const size_t blocks = x_size - (x_size % 8);
421
422 mp_word_t carry = 0;
423
424 for(size_t i = 0; i != blocks; i += 8) {
425 carry = word8_linmul2(x + i, y, carry);
426 }
427 for(size_t i = blocks; i != x_size; ++i) {
428 x[i] = word_madd2(x[i], y, carry);
429 }
430 return carry;
431 }
432
433 inline void bigint_linmul3(mp_word_t z[], const mp_word_t x[], size_t x_size, mp_word_t y) noexcept {
434 const size_t blocks = x_size - (x_size % 8);
435
436 mp_word_t carry = 0;
437
438 for(size_t i = 0; i != blocks; i += 8) {
439 carry = word8_linmul3(z + i, x + i, y, carry);
440 }
441 for(size_t i = blocks; i != x_size; ++i) {
442 z[i] = word_madd2(x[i], y, carry);
443 }
444 z[x_size] = carry;
445 }
446
447 inline void bigint_shl1(mp_word_t x[], size_t x_size, size_t x_words, size_t word_shift, size_t bit_shift) noexcept {
448 std::memmove(x + word_shift, x, sizeof(mp_word_t)*x_words);
449 std::memset(x, 0, sizeof(mp_word_t)*word_shift);
450
451 const auto carry_mask = CT::Mask<mp_word_t>::expand(bit_shift);
452 const size_t carry_shift = carry_mask.if_set_return(mp_word_bits - bit_shift);
453
454 mp_word_t carry = 0;
455 for(size_t i = word_shift; i != x_size; ++i) {
456 const mp_word_t w = x[i];
457 x[i] = (w << bit_shift) | carry;
458 carry = carry_mask.if_set_return(w >> carry_shift);
459 }
460 }
461
462 inline void bigint_shr1(mp_word_t x[], size_t x_size, size_t word_shift, size_t bit_shift) noexcept {
463 const size_t top = x_size >= word_shift ? (x_size - word_shift) : 0;
464
465 if(top > 0) {
466 std::memmove(x, x + word_shift, sizeof(mp_word_t)*top);
467 }
468 std::memset(x + top, 0, sizeof(mp_word_t)*jau::min(word_shift, x_size));
469
470 const auto carry_mask = CT::Mask<mp_word_t>::expand(bit_shift);
471 const size_t carry_shift = carry_mask.if_set_return(mp_word_bits - bit_shift);
472
473 mp_word_t carry = 0;
474 for(size_t i = 0; i != top; ++i) {
475 const mp_word_t w = x[top - i - 1];
476 x[top-i-1] = (w >> bit_shift) | carry;
477 carry = carry_mask.if_set_return(w << carry_shift);
478 }
479 }
480
481 inline 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 {
482 std::memmove(y + word_shift, x, sizeof(mp_word_t)*x_size);
483
484 const auto carry_mask = CT::Mask<mp_word_t>::expand(bit_shift);
485 const size_t carry_shift = carry_mask.if_set_return(mp_word_bits - bit_shift);
486
487 mp_word_t carry = 0;
488 for(size_t i = word_shift; i != x_size + word_shift + 1; ++i) {
489 const mp_word_t w = y[i];
490 y[i] = (w << bit_shift) | carry;
491 carry = carry_mask.if_set_return(w >> carry_shift);
492 }
493 }
494 inline 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 {
495 const size_t new_size = x_size < word_shift ? 0 : (x_size - word_shift);
496
497 if(new_size > 0) {
498 std::memmove(y, x + word_shift, sizeof(mp_word_t)*new_size);
499 }
500 const auto carry_mask = CT::Mask<mp_word_t>::expand(bit_shift);
501 const size_t carry_shift = carry_mask.if_set_return(mp_word_bits - bit_shift);
502
503 mp_word_t carry = 0;
504 for(size_t i = new_size; i > 0; --i) {
505 mp_word_t w = y[i-1];
506 y[i-1] = (w >> bit_shift) | carry;
507 carry = carry_mask.if_set_return(w << carry_shift);
508 }
509 }
510
511 /*
512 * O(n*n) multiplication
513 */
514 inline void basecase_mul(mp_word_t z[], size_t z_size,
515 const mp_word_t x[], size_t x_size,
516 const mp_word_t y[], size_t y_size) noexcept
517 {
518 assert(z_size >= x_size + y_size); // z_size too small
519
520 const size_t x_size_8 = x_size - (x_size % 8);
521
522 for(size_t i=0; i<z_size; ++i) { z[i]=0; }
523
524 for(size_t i = 0; i != y_size; ++i) {
525 const mp_word_t y_i = y[i];
526
527 mp_word_t carry = 0;
528
529 for(size_t j = 0; j != x_size_8; j += 8) {
530 carry = word8_madd3(z + i + j, x + j, y_i, carry);
531 }
532 for(size_t j = x_size_8; j != x_size; ++j) {
533 z[i+j] = word_madd3(x[j], y_i, z[i+j], carry);
534 }
535 z[x_size+i] = carry;
536 }
537 }
538
539 /** Computes ((n1<<bits) + n0) / d */
541 if(d == 0) {
543 }
544 if constexpr ( has_mp_dword ) {
545 return static_cast<mp_word_t>( ( ( static_cast<mp_dword_t>(n1) << mp_word_bits ) | n0 ) / d );
546 } else {
547 mp_word_t high = n1 % d;
548 mp_word_t quotient = 0;
549
550 for(size_t i = 0; i != mp_word_bits; ++i) {
551 const mp_word_t high_top_bit = high >> ( mp_word_bits - 1 );
552
553 high <<= 1;
554 high |= (n0 >> (mp_word_bits-1-i)) & 1;
555 quotient <<= 1;
556
557 if(high_top_bit || high >= d) {
558 high -= d;
559 quotient |= 1;
560 }
561 }
562 return quotient;
563 }
564 }
565
566 /** Compute ((n1<<bits) + n0) % d */
568 if(d == 0) {
570 }
571 if constexpr ( has_mp_dword ) {
572 return ( ( static_cast<mp_dword_t>(n1) << mp_word_bits) | n0 ) % d;
573 } else {
574 mp_word_t z = bigint_divop(n1, n0, d);
575 mp_word_t dummy = 0;
576 z = word_madd2(z, d, dummy);
577 return (n0-z);
578 }
579 }
580
582 bigint_ct_is_eq(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept {
583 const size_t common_elems = std::min(x_size, y_size);
584 mp_word_t diff = 0;
585
586 for(size_t i = 0; i != common_elems; i++) {
587 diff |= (x[i] ^ y[i]);
588 }
589
590 // If any bits were set in high part of x/y, then they are not equal
591 if(x_size < y_size)
592 {
593 for(size_t i = x_size; i != y_size; i++) {
594 diff |= y[i];
595 }
596 } else if(y_size < x_size) {
597 for(size_t i = y_size; i != x_size; i++) {
598 diff |= x[i];
599 }
600 }
601 return CT::Mask<mp_word_t>::is_zero(diff);
602 }
603
604 /**
605 * Compare x and y
606 * Return ~0 if x[0:x_size] < y[0:y_size] or 0 otherwise
607 * If lt_or_equal is true, returns ~0 also for x == y
608 */
610 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
611 {
612 const size_t common_elems = jau::min(x_size, y_size);
613
614 auto is_lt = CT::Mask<mp_word_t>::expand(lt_or_equal);
615
616 for(size_t i = 0; i != common_elems; i++)
617 {
618 const auto eq = CT::Mask<mp_word_t>::is_equal(x[i], y[i]);
619 const auto lt = CT::Mask<mp_word_t>::is_lt(x[i], y[i]);
620 is_lt = eq.select_mask(is_lt, lt);
621 }
622
623 if(x_size < y_size) {
624 mp_word_t mask = 0;
625 for(size_t i = x_size; i != y_size; i++) {
626 mask |= y[i];
627 }
628 // If any bits were set in high part of y, then is_lt should be forced true
629 is_lt |= CT::Mask<mp_word_t>::expand(mask);
630 } else if(y_size < x_size) {
631 mp_word_t mask = 0;
632 for(size_t i = y_size; i != x_size; i++) {
633 mask |= x[i];
634 }
635
636 // If any bits were set in high part of x, then is_lt should be false
637 is_lt &= CT::Mask<mp_word_t>::is_zero(mask);
638 }
639
640 return is_lt;
641 }
642
643 /**
644 * Compare unsigned x and y mp_word_t
645 *
646 * Returns
647 * - -1 if x < y
648 * - 0 if x == y
649 * - 1 if x > y
650 */
651 inline int bigint_cmp(const mp_word_t x[], size_t x_size, const mp_word_t y[], size_t y_size) noexcept
652 {
653 constexpr const mp_word_t LT = static_cast<mp_word_t>(-1);
654 constexpr const mp_word_t EQ = 0;
655 constexpr const mp_word_t GT = 1;
656 const size_t common_elems = jau::min(x_size, y_size);
657
658 mp_word_t result = EQ; // until found otherwise
659
660 for(size_t i = 0; i != common_elems; i++) {
661 const auto is_eq = CT::Mask<mp_word_t>::is_equal(x[i], y[i]);
662 const auto is_lt = CT::Mask<mp_word_t>::is_lt(x[i], y[i]);
663 result = is_eq.select(result, is_lt.select(LT, GT));
664 }
665 if(x_size < y_size) {
666 mp_word_t mask = 0;
667 for(size_t i = x_size; i != y_size; i++) {
668 mask |= y[i];
669 }
670 // If any bits were set in high part of y, then x < y
671 result = CT::Mask<mp_word_t>::is_zero(mask).select(result, LT);
672 } else if(y_size < x_size) {
673 mp_word_t mask = 0;
674 for(size_t i = y_size; i != x_size; i++) {
675 mask |= x[i];
676 }
677 // If any bits were set in high part of x, then x > y
678 result = CT::Mask<mp_word_t>::is_zero(mask).select(result, GT);
679 }
680 CT::unpoison(result);
681 return static_cast<int>(result);
682 }
683
684} // namespace jau::mp::ops
685
686#endif /* JAU_BIG_INT_OPS_HPP_ */
#define E_FILE_LINE
A Mask type used for constant-time operations.
Definition ct_utils.hpp:81
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
static Mask< T > is_zero(T x) noexcept
Return a Mask<T> which is set if v is == 0 or cleared otherwise.
Definition ct_utils.hpp:137
static Mask< T > is_lt(T x, T y) noexcept
Return a Mask<T> which is set if x < y.
Definition ct_utils.hpp:153
math_error_t::div_by_zero, i.e.
consteval_cxx20 bool is_builtin_int128_available() noexcept
constexpr T min(const T x, const T y) noexcept
Returns the minimum of two integrals (w/ branching) in O(1)
constexpr size_t pointer_bit_size() noexcept
Returns the compile time pointer architecture size in bits.
Definition cpuid.hpp:48
void unpoison(const T *p, size_t n)
Definition ct_utils.hpp:54
void conditional_swap(bool cnd, T &x, T &y) noexcept
Definition ct_utils.hpp:369
void conditional_swap_ptr(bool cnd, T &x, T &y) noexcept
Definition ct_utils.hpp:380
constexpr size_t best_word_byte_size()
mp_word_t word8_madd3(mp_word_t z[8], const mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept
mp_word_t word8_sub2(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept
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
mp_word_t word8_sub2_rev(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) 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.
mp_word_t word8_linmul3(mp_word_t z[8], const mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept
mp_word_t word_sub(mp_word_t x, mp_word_t y, mp_word_t &carry) noexcept
mp_word_t word8_sub3(mp_word_t z[8], const mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept
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.
mp_word_t bigint_add3_nc(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 with carry out.
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.
void mul64x64_128(const uint64_t a, const uint64_t b, uint64_t &lo, uint64_t &hi) noexcept
64x64->128 bit multiplication
mp_word_t word8_linmul2(mp_word_t x[8], mp_word_t y, mp_word_t carry) noexcept
mp_word_t word8_add3(mp_word_t z[8], const mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) noexcept
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_madd3(mp_word_t a, mp_word_t b, mp_word_t c, mp_word_t &d) noexcept
Word Multiply/Add.
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
mp_word_t word8_add2(mp_word_t x[8], const mp_word_t y[8], mp_word_t carry) 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
mp_word_t word_add(const mp_word_t x, const mp_word_t y, mp_word_t &carry) 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
constexpr const bool has_mp_dword
jau::uint_bytes< impl::best_word_byte_size()>::type mp_word_t
constexpr const size_t mp_word_bits
jau::uint_bytes< impl::best_word_byte_size() *2 >::type mp_dword_t
constexpr const mp_word_t mp_word_max