diff options
| author | Jeroen van Rijn <Kelimion@users.noreply.github.com> | 2021-09-01 00:04:55 +0200 |
|---|---|---|
| committer | Jeroen van Rijn <Kelimion@users.noreply.github.com> | 2021-09-01 19:13:47 +0200 |
| commit | a056e1943435ed37330bc6d5c0eac241d323170b (patch) | |
| tree | 12b9347795c180e4f8f5e299c8754a192d70e54d /core/math/big/private.odin | |
| parent | 7d0dedf951b615835cd84819460f3f29a89a5c9b (diff) | |
big: Cue up `internal_int_exponent_mod` wrapper function.
Diffstat (limited to 'core/math/big/private.odin')
| -rw-r--r-- | core/math/big/private.odin | 1009 |
1 files changed, 1009 insertions, 0 deletions
diff --git a/core/math/big/private.odin b/core/math/big/private.odin index 552b100cf..d2878fcc1 100644 --- a/core/math/big/private.odin +++ b/core/math/big/private.odin @@ -1851,6 +1851,1015 @@ _private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := co }
/*
+ Computes xR**-1 == x (mod N) via Montgomery Reduction.
+ Assumes `x` and `n` not to be nil.
+*/
+_private_int_montgomery_reduce :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+ /*
+ Can the fast reduction [comba] method be used?
+ Note that unlike in mul, you're safely allowed *less* than the available columns [255 per default],
+ since carries are fixed up in the inner loop.
+ */
+ internal_clear_if_uninitialized(x, n) or_return;
+
+ digs := (n.used * 2) + 1;
+ if digs < _WARRAY && x.used <= _WARRAY && n.used < _MAX_COMBA {
+ return _private_montgomery_reduce_comba(x, n, rho);
+ }
+
+ /*
+ Grow the input as required
+ */
+ internal_grow(x, digs) or_return;
+ x.used = digs;
+
+ for ix := 0; ix < n.used; ix += 1 {
+ /*
+ `mu = ai * rho mod b`
+ The value of rho must be precalculated via `int_montgomery_setup()`,
+ such that it equals -1/n0 mod b this allows the following inner loop
+ to reduce the input one digit at a time.
+ */
+
+ mu := DIGIT((_WORD(x.digit[ix]) * _WORD(rho)) & _WORD(_MASK));
+
+ /*
+ a = a + mu * m * b**i
+ Multiply and add in place.
+ */
+ u := DIGIT(0);
+ iy := int(0);
+ for ; iy < n.used; iy += 1 {
+ /*
+ Compute product and sum.
+ */
+ r := (_WORD(mu) * _WORD(n.digit[iy]) + _WORD(u) + _WORD(x.digit[ix + iy]));
+
+ /*
+ Get carry.
+ */
+ u = DIGIT(r >> _DIGIT_BITS);
+
+ /*
+ Fix digit.
+ */
+ x.digit[ix + iy] = DIGIT(r & _WORD(_MASK));
+ }
+
+ /*
+ At this point the ix'th digit of x should be zero.
+ Propagate carries upwards as required.
+ */
+ for u != 0 {
+ x.digit[ix + iy] += u;
+ u = x.digit[ix + iy] >> _DIGIT_BITS;
+ x.digit[ix + iy] &= _MASK;
+ iy += 1;
+ }
+ }
+
+ /*
+ At this point the n.used'th least significant digits of x are all zero,
+ which means we can shift x to the right by n.used digits and the
+ residue is unchanged.
+
+ x = x/b**n.used.
+ */
+ internal_clamp(x);
+ internal_shr_digit(x, n.used);
+
+ /*
+ if x >= n then x = x - n
+ */
+ if internal_cmp_mag(x, n) != -1 {
+ return internal_sub(x, x, n);
+ }
+
+ return nil;
+}
+
+/*
+ Shifts with subtractions when the result is greater than b.
+
+ The method is slightly modified to shift B unconditionally upto just under
+ the leading bit of b. This saves alot of multiple precision shifting.
+
+ Assumes `a` and `b` not to be `nil`.
+*/
+_private_int_montgomery_calc_normalization :: proc(a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+ /*
+ How many bits of last digit does b use.
+ */
+ internal_clear_if_uninitialized(a, b) or_return;
+
+ bits := internal_count_bits(b) % _DIGIT_BITS;
+
+ if b.used > 1 {
+ power := ((b.used - 1) * _DIGIT_BITS) + bits - 1;
+ internal_int_power_of_two(a, power) or_return;
+ } else {
+ internal_one(a) or_return;
+ bits = 1;
+ }
+
+ /*
+ Now compute C = A * B mod b.
+ */
+ for x := bits - 1; x < _DIGIT_BITS; x += 1 {
+ internal_int_shl1(a, a) or_return;
+ if internal_cmp_mag(a, b) != -1 {
+ internal_sub(a, a, b) or_return;
+ }
+ }
+ return nil;
+}
+
+/*
+ Sets up the Montgomery reduction stuff.
+*/
+_private_int_montgomery_setup :: proc(n: ^Int, allocator := context.allocator) -> (rho: DIGIT, err: Error) {
+ /*
+ Fast inversion mod 2**k
+ Based on the fact that:
+
+ XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
+ => 2*X*A - X*X*A*A = 1
+ => 2*(1) - (1) = 1
+ */
+ internal_clear_if_uninitialized(n, allocator) or_return;
+
+ b := n.digit[0];
+ if b & 1 == 0 { return 0, .Invalid_Argument; }
+
+ x := (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
+ x *= 2 - (b * x); /* here x*a==1 mod 2**8 */
+ x *= 2 - (b * x); /* here x*a==1 mod 2**16 */
+
+ when _DIGIT_TYPE_BITS == 64 {
+ x *= 2 - (b * x); /* here x*a==1 mod 2**32 */
+ x *= 2 - (b * x); /* here x*a==1 mod 2**64 */
+ }
+
+ /*
+ rho = -1/m mod b
+ */
+ rho = DIGIT(((_WORD(1) << _WORD(_DIGIT_BITS)) - _WORD(x)) & _WORD(_MASK));
+ return rho, nil;
+}
+
+/*
+ Reduces `x` mod `m`, assumes 0 < x < m**2, mu is precomputed via reduce_setup.
+ From HAC pp.604 Algorithm 14.42
+
+ Assumes `x`, `m` and `mu` all not to be `nil` and have been initialized.
+*/
+_private_int_reduce :: proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ q := &Int{};
+ defer internal_destroy(q);
+ um := m.used;
+
+ /*
+ q = x
+ */
+ internal_copy(q, x) or_return;
+
+ /*
+ q1 = x / b**(k-1)
+ */
+ internal_shr_digit(q, um - 1);
+
+ /*
+ According to HAC this optimization is ok.
+ */
+ if DIGIT(um) > DIGIT(1) << (_DIGIT_BITS - 1) {
+ internal_mul(q, q, mu) or_return;
+ } else {
+ _private_int_mul_high(q, q, mu, um) or_return;
+ }
+
+ /*
+ q3 = q2 / b**(k+1)
+ */
+ internal_shr_digit(q, um + 1);
+
+ /*
+ x = x mod b**(k+1), quick (no division)
+ */
+ internal_int_mod_bits(x, x, _DIGIT_BITS * (um + 1)) or_return;
+
+ /*
+ q = q * m mod b**(k+1), quick (no division)
+ */
+ _private_int_mul(q, q, m, um + 1) or_return;
+
+ /*
+ x = x - q
+ */
+ internal_sub(x, x, q) or_return;
+
+ /*
+ If x < 0, add b**(k+1) to it.
+ */
+ if internal_cmp(x, 0) == -1 {
+ internal_set(q, 1) or_return;
+ internal_shl_digit(q, um + 1) or_return;
+ internal_add(x, x, q) or_return;
+ }
+
+ /*
+ Back off if it's too big.
+ */
+ for internal_cmp(x, m) != -1 {
+ internal_sub(x, x, m) or_return;
+ }
+
+ return nil;
+}
+
+/*
+ Reduces `a` modulo `n`, where `n` is of the form 2**p - d.
+*/
+_private_int_reduce_2k :: proc(a, n: ^Int, d: DIGIT, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ q := &Int{};
+ defer internal_destroy(q);
+
+ internal_zero(q) or_return;
+
+ p := internal_count_bits(n);
+
+ for {
+ /*
+ q = a/2**p, a = a mod 2**p
+ */
+ internal_shrmod(q, a, a, p) or_return;
+
+ if d != 1 {
+ /*
+ q = q * d
+ */
+ internal_mul(q, q, d) or_return;
+ }
+
+ /*
+ a = a + q
+ */
+ internal_add(a, a, q) or_return;
+ if internal_cmp_mag(a, n) == -1 { break; }
+ internal_sub(a, a, n) or_return;
+ }
+
+ return nil;
+}
+
+/*
+ Reduces `a` modulo `n` where `n` is of the form 2**p - d
+ This differs from reduce_2k since "d" can be larger than a single digit.
+*/
+_private_int_reduce_2k_l :: proc(a, n, d: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ q := &Int{};
+ defer internal_destroy(q);
+
+ internal_zero(q) or_return;
+
+ p := internal_count_bits(n);
+
+ for {
+ /*
+ q = a/2**p, a = a mod 2**p
+ */
+ internal_shrmod(q, a, a, p) or_return;
+
+ /*
+ q = q * d
+ */
+ internal_mul(q, q, d) or_return;
+
+ /*
+ a = a + q
+ */
+ internal_add(a, a, q) or_return;
+ if internal_cmp_mag(a, n) == -1 { break; }
+ internal_sub(a, a, n) or_return;
+ }
+
+ return nil;
+}
+
+/*
+ Determines if `internal_int_reduce_2k` can be used.
+ Asssumes `a` not to be `nil` and to have been initialized.
+*/
+_private_int_reduce_is_2k :: proc(a: ^Int) -> (reducible: bool, err: Error) {
+ assert_if_nil(a);
+
+ if internal_is_zero(a) {
+ return false, nil;
+ } else if a.used == 1 {
+ return true, nil;
+ } else if a.used > 1 {
+ iy := internal_count_bits(a);
+ iw := 1;
+ iz := DIGIT(1);
+
+ /*
+ Test every bit from the second digit up, must be 1.
+ */
+ for ix := _DIGIT_BITS; ix < iy; ix += 1 {
+ if a.digit[iw] & iz == 0 {
+ return false, nil;
+ }
+
+ iz <<= 1;
+ if iz > _DIGIT_MAX {
+ iw += 1;
+ iz = 1;
+ }
+ }
+ return true, nil;
+ } else {
+ return true, nil;
+ }
+}
+
+/*
+ Determines if `internal_int_reduce_2k_l` can be used.
+ Asssumes `a` not to be `nil` and to have been initialized.
+*/
+_private_int_reduce_is_2k_l :: proc(a: ^Int) -> (reducible: bool, err: Error) {
+ assert_if_nil(a);
+
+ if internal_int_is_zero(a) {
+ return false, nil;
+ } else if a.used == 1 {
+ return true, nil;
+ } else if a.used > 1 {
+ /*
+ If more than half of the digits are -1 we're sold.
+ */
+ ix := 0;
+ iy := 0;
+
+ for ; ix < a.used; ix += 1 {
+ if a.digit[ix] == _DIGIT_MAX {
+ iy += 1;
+ }
+ }
+ return iy >= (a.used / 2), nil;
+ } else {
+ return false, nil;
+ }
+}
+
+/*
+ Determines the setup value.
+ Assumes `a` is not `nil`.
+*/
+_private_int_reduce_2k_setup :: proc(a: ^Int, allocator := context.allocator) -> (d: DIGIT, err: Error) {
+ context.allocator = allocator;
+
+ tmp := &Int{};
+ defer internal_destroy(tmp);
+ internal_zero(tmp) or_return;
+
+ internal_int_power_of_two(tmp, internal_count_bits(a)) or_return;
+ internal_sub(tmp, tmp, a) or_return;
+
+ return tmp.digit[0], nil;
+}
+
+/*
+ Determines the setup value.
+ Assumes `mu` and `P` are not `nil`.
+
+ d := (1 << a.bits) - a;
+*/
+_private_int_reduce_2k_setup_l :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ tmp := &Int{};
+ defer internal_destroy(tmp);
+ internal_zero(tmp) or_return;
+
+ internal_int_power_of_two(tmp, internal_count_bits(P)) or_return;
+ internal_sub(mu, tmp, P) or_return;
+
+ return nil;
+}
+
+/*
+ Pre-calculate the value required for Barrett reduction.
+ For a given modulus "P" it calulates the value required in "mu"
+ Assumes `mu` and `P` are not `nil`.
+*/
+_private_int_reduce_setup :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ internal_int_power_of_two(mu, P.used * 2 * _DIGIT_BITS) or_return;
+ return internal_int_div(mu, mu, P);
+}
+
+/*
+ Determines the setup value.
+ Assumes `a` to not be `nil` and to have been initialized.
+*/
+_private_int_dr_setup :: proc(a: ^Int) -> (d: DIGIT) {
+ /*
+ The casts are required if _DIGIT_BITS is one less than
+ the number of bits in a DIGIT [e.g. _DIGIT_BITS==31].
+ */
+ return DIGIT((1 << _DIGIT_BITS) - a.digit[0]);
+}
+
+/*
+ Determines if a number is a valid DR modulus.
+ Assumes `a` to not be `nil` and to have been initialized.
+*/
+_private_dr_is_modulus :: proc(a: ^Int) -> (res: bool) {
+ /*
+ Must be at least two digits.
+ */
+ if a.used < 2 { return false; }
+
+ /*
+ Must be of the form b**k - a [a <= b] so all but the first digit must be equal to -1 (mod b).
+ */
+ for ix := 1; ix < a.used; ix += 1 {
+ if a.digit[ix] != _MASK {
+ return false;
+ }
+ }
+ return true;
+}
+
+/*
+ Reduce "x" in place modulo "n" using the Diminished Radix algorithm.
+ Based on algorithm from the paper
+
+ "Generating Efficient Primes for Discrete Log Cryptosystems"
+ Chae Hoon Lim, Pil Joong Lee,
+ POSTECH Information Research Laboratories
+
+ The modulus must be of a special format [see manual].
+ Has been modified to use algorithm 7.10 from the LTM book instead
+
+ Input x must be in the range 0 <= x <= (n-1)**2
+ Assumes `x` and `n` to not be `nil` and to have been initialized.
+*/
+_private_int_dr_reduce :: proc(x, n: ^Int, k: DIGIT, allocator := context.allocator) -> (err: Error) {
+ /*
+ m = digits in modulus.
+ */
+ m := n.used;
+
+ /*
+ Ensure that "x" has at least 2m digits.
+ */
+ internal_grow(x, m + m) or_return;
+
+ /*
+ Top of loop, this is where the code resumes if another reduction pass is required.
+ */
+ for {
+ i: int;
+ mu := DIGIT(0);
+
+ /*
+ Compute (x mod B**m) + k * [x/B**m] inline and inplace.
+ */
+ for i = 0; i < m; i += 1 {
+ r := _WORD(x.digit[i + m]) * _WORD(k) + _WORD(x.digit[i] + mu);
+ x.digit[i] = DIGIT(r & _WORD(_MASK));
+ mu = DIGIT(r >> _WORD(_DIGIT_BITS));
+ }
+
+ /*
+ Set final carry.
+ */
+ x.digit[i] = mu;
+
+ /*
+ Zero words above m.
+ */
+ mem.zero_slice(x.digit[m + 1:][:x.used - m]);
+
+ /*
+ Clamp, sub and return.
+ */
+ internal_clamp(x) or_return;
+
+ /*
+ If x >= n then subtract and reduce again.
+ Each successive "recursion" makes the input smaller and smaller.
+ */
+ if internal_cmp_mag(x, n) == -1 { break; }
+
+ internal_sub(x, x, n) or_return;
+ }
+ return nil;
+}
+
+/*
+ Computes res == G**X mod P.
+ Assumes `res`, `G`, `X` and `P` to not be `nil` and for `G`, `X` and `P` to have been initialized.
+*/
+_private_int_exponent_mod :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ M := [_TAB_SIZE]Int{};
+ winsize: uint;
+
+ /*
+ Use a pointer to the reduction algorithm.
+ This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.
+ */
+ redux: #type proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error);
+
+ defer {
+ internal_destroy(&M[1]);
+ for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+ internal_destroy(&M[x]);
+ }
+ }
+
+ /*
+ Find window size.
+ */
+ x := internal_count_bits(X);
+ switch {
+ case x <= 7:
+ winsize = 2;
+ case x <= 36:
+ winsize = 3;
+ case x <= 140:
+ winsize = 4;
+ case x <= 450:
+ winsize = 5;
+ case x <= 1303:
+ winsize = 6;
+ case x <= 3529:
+ winsize = 7;
+ case:
+ winsize = 8;
+ }
+
+ winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize;
+
+ /*
+ Init M array.
+ Init first cell.
+ */
+ internal_zero(&M[1]) or_return;
+
+ /*
+ Now init the second half of the array.
+ */
+ for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+ internal_zero(&M[x]) or_return;
+ }
+
+ /*
+ Create `mu`, used for Barrett reduction.
+ */
+ mu := &Int{};
+ defer internal_destroy(mu);
+ internal_zero(mu) or_return;
+
+ if redmode == 0 {
+ _private_int_reduce_setup(mu, P) or_return;
+ redux = _private_int_reduce;
+ } else {
+ _private_int_reduce_2k_setup_l(mu, P) or_return;
+ redux = _private_int_reduce_2k_l;
+ }
+
+ /*
+ Create M table.
+
+ The M table contains powers of the base, e.g. M[x] = G**x mod P.
+ The first half of the table is not computed, though, except for M[0] and M[1].
+ */
+ internal_int_mod(&M[1], G, P) or_return;
+
+ /*
+ Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.
+
+ TODO: This can probably be replaced by computing the power and using `pow` to raise to it
+ instead of repeated squaring.
+ */
+ slot := 1 << (winsize - 1);
+ internal_copy(&M[slot], &M[1]) or_return;
+
+ for x = 0; x < int(winsize - 1); x += 1 {
+ /*
+ Square it.
+ */
+ internal_sqr(&M[slot], &M[slot]) or_return;
+
+ /*
+ Reduce modulo P
+ */
+ redux(&M[slot], P, mu) or_return;
+ }
+
+ /*
+ Create upper table, that is M[x] = M[x-1] * M[1] (mod P)
+ for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
+ */
+ for x = slot + 1; x < (1 << winsize); x += 1 {
+ internal_mul(&M[x], &M[x - 1], &M[1]) or_return;
+ redux(&M[x], P, mu) or_return;
+ }
+
+ /*
+ Setup result.
+ */
+ internal_one(res) or_return;
+
+ /*
+ Set initial mode and bit cnt.
+ */
+ mode := 0;
+ bitcnt := 1;
+ buf := DIGIT(0);
+ digidx := X.used - 1;
+ bitcpy := uint(0);
+ bitbuf := DIGIT(0);
+
+ for {
+ /*
+ Grab next digit as required.
+ */
+ bitcnt -= 1;
+ if bitcnt == 0 {
+ /*
+ If digidx == -1 we are out of digits.
+ */
+ if digidx == -1 { break; }
+
+ /*
+ Read next digit and reset the bitcnt.
+ */
+ buf = X.digit[digidx];
+ digidx -= 1;
+ bitcnt = _DIGIT_BITS;
+ }
+
+ /*
+ Grab the next msb from the exponent.
+ */
+ y := buf >> (_DIGIT_BITS - 1) & 1;
+ buf <<= 1;
+
+ /*
+ If the bit is zero and mode == 0 then we ignore it.
+ These represent the leading zero bits before the first 1 bit
+ in the exponent. Technically this opt is not required but it
+ does lower the # of trivial squaring/reductions used.
+ */
+ if mode == 0 && y == 0 {
+ continue;
+ }
+
+ /*
+ If the bit is zero and mode == 1 then we square.
+ */
+ if mode == 1 && y == 0 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, mu) or_return;
+ continue;
+ }
+
+ /*
+ Else we add it to the window.
+ */
+ bitcpy += 1;
+ bitbuf |= (y << (winsize - bitcpy));
+ mode = 2;
+
+ if (bitcpy == winsize) {
+ /*
+ Window is filled so square as required and multiply.
+ Square first.
+ */
+ for x = 0; x < int(winsize); x += 1 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, mu) or_return;
+ }
+
+ /*
+ Then multiply.
+ */
+ internal_mul(res, res, &M[bitbuf]) or_return;
+ redux(res, P, mu) or_return;
+
+ /*
+ Empty window and reset.
+ */
+ bitcpy = 0;
+ bitbuf = 0;
+ mode = 1;
+ }
+ }
+
+ /*
+ If bits remain then square/multiply.
+ */
+ if mode == 2 && bitcpy > 0 {
+ /*
+ Square then multiply if the bit is set.
+ */
+ for x = 0; x < int(bitcpy); x += 1 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, mu) or_return;
+
+ bitbuf <<= 1;
+ if ((bitbuf & (1 << winsize)) != 0) {
+ /*
+ Then multiply.
+ */
+ internal_mul(res, res, &M[1]) or_return;
+ redux(res, P, mu) or_return;
+ }
+ }
+ }
+ return err;
+}
+
+/*
+ Computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
+
+ Uses a left-to-right `k`-ary sliding window to compute the modular exponentiation.
+ The value of `k` changes based on the size of the exponent.
+
+ Uses Montgomery or Diminished Radix reduction [whichever appropriate]
+
+ Assumes `res`, `G`, `X` and `P` to not be `nil` and for `G`, `X` and `P` to have been initialized.
+*/
+_private_int_exponent_mod_fast :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ M := [_TAB_SIZE]Int{};
+ winsize: uint;
+
+ /*
+ Use a pointer to the reduction algorithm.
+ This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.
+ */
+ redux: #type proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error);
+
+ defer {
+ internal_destroy(&M[1]);
+ for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+ internal_destroy(&M[x]);
+ }
+ }
+
+ /*
+ Find window size.
+ */
+ x := internal_count_bits(X);
+ switch {
+ case x <= 7:
+ winsize = 2;
+ case x <= 36:
+ winsize = 3;
+ case x <= 140:
+ winsize = 4;
+ case x <= 450:
+ winsize = 5;
+ case x <= 1303:
+ winsize = 6;
+ case x <= 3529:
+ winsize = 7;
+ case:
+ winsize = 8;
+ }
+
+ winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize;
+
+ /*
+ Init M array
+ Init first cell.
+ */
+ cap := internal_int_allocated_cap(P);
+ internal_grow(&M[1], cap) or_return;
+
+ /*
+ Now init the second half of the array.
+ */
+ for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+ internal_grow(&M[x], cap) or_return;
+ }
+
+ /*
+ Determine and setup reduction code.
+ */
+ rho: DIGIT;
+
+ if redmode == 0 {
+ /*
+ Now setup Montgomery.
+ */
+ rho = _private_int_montgomery_setup(P) or_return;
+
+ /*
+ Automatically pick the comba one if available (saves quite a few calls/ifs).
+ */
+ if ((P.used * 2) + 1) < _WARRAY && P.used < _MAX_COMBA {
+ redux = _private_montgomery_reduce_comba;
+ } else {
+ /*
+ Use slower baseline Montgomery method.
+ */
+ redux = _private_int_montgomery_reduce;
+ }
+ } else if redmode == 1 {
+ /*
+ Setup DR reduction for moduli of the form B**k - b.
+ */
+ rho = _private_int_dr_setup(P);
+ redux = _private_int_dr_reduce;
+ } else {
+ /*
+ Setup DR reduction for moduli of the form 2**k - b.
+ */
+ rho = _private_int_reduce_2k_setup(P) or_return;
+ redux = _private_int_reduce_2k;
+ }
+
+ /*
+ Setup result.
+ */
+ internal_grow(res, cap) or_return;
+
+ /*
+ Create M table
+ The first half of the table is not computed, though, except for M[0] and M[1]
+ */
+
+ if redmode == 0 {
+ /*
+ Now we need R mod m.
+ */
+ _private_int_montgomery_calc_normalization(res, P) or_return;
+
+ /*
+ Now set M[1] to G * R mod m.
+ */
+ internal_mulmod(&M[1], G, res, P) or_return;
+ } else {
+ internal_one(res) or_return;
+ internal_mod(&M[1], G, P) or_return;
+ }
+
+ /*
+ Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.
+ */
+ slot := 1 << (winsize - 1);
+ internal_copy(&M[slot], &M[1]) or_return;
+
+ for x = 0; x < int(winsize - 1); x += 1 {
+ internal_sqr(&M[slot], &M[slot]) or_return;
+ redux(&M[slot], P, rho) or_return;
+ }
+
+ /*
+ Create upper table.
+ */
+ for x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x += 1 {
+ internal_mul(&M[x], &M[x - 1], &M[1]) or_return;
+ redux(&M[x], P, rho) or_return;
+ }
+
+ /*
+ Set initial mode and bit cnt.
+ */
+ mode := 0;
+ bitcnt := 1;
+ buf := DIGIT(0);
+ digidx := X.used - 1;
+ bitcpy := 0;
+ bitbuf := DIGIT(0);
+
+ for {
+ /*
+ Grab next digit as required.
+ */
+ bitcnt -= 1;
+ if bitcnt == 0 {
+ /*
+ If digidx == -1 we are out of digits so break.
+ */
+ if digidx == -1 { break; }
+
+ /*
+ Read next digit and reset the bitcnt.
+ */
+ buf = X.digit[digidx];
+ digidx -= 1;
+ bitcnt = _DIGIT_BITS;
+ }
+
+ /*
+ Grab the next msb from the exponent.
+ */
+ y := (buf >> (_DIGIT_BITS - 1)) & 1;
+ buf <<= 1;
+
+ /*
+ If the bit is zero and mode == 0 then we ignore it.
+ These represent the leading zero bits before the first 1 bit in the exponent.
+ Technically this opt is not required but it does lower the # of trivial squaring/reductions used.
+ */
+ if mode == 0 && y == 0 { continue; }
+
+ /*
+ If the bit is zero and mode == 1 then we square.
+ */
+ if mode == 1 && y == 0 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, rho) or_return;
+ continue;
+ }
+
+ /*
+ Else we add it to the window.
+ */
+ bitcpy += 1;
+ bitbuf |= (y << (winsize - uint(bitcpy)));
+ mode = 2;
+
+ if bitcpy == int(winsize) {
+ /*
+ Window is filled so square as required and multiply
+ Square first.
+ */
+ for x = 0; x < int(winsize); x += 1 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, rho) or_return;
+ }
+
+ /*
+ Then multiply.
+ */
+ internal_mul(res, res, &M[bitbuf]) or_return;
+ redux(res, P, rho) or_return;
+
+ /*
+ Empty window and reset.
+ */
+ bitcpy = 0;
+ bitbuf = 0;
+ mode = 1;
+ }
+ }
+
+ /*
+ If bits remain then square/multiply.
+ */
+ if mode == 2 && bitcpy > 0 {
+ /*
+ Square then multiply if the bit is set.
+ */
+ for x = 0; x < bitcpy; x += 1 {
+ internal_sqr(res, res) or_return;
+ redux(res, P, rho) or_return;
+
+ /*
+ Get next bit of the window.
+ */
+ bitbuf <<= 1;
+ if bitbuf & (1 << winsize) != 0 {
+ /*
+ Then multiply.
+ */
+ internal_mul(res, res, &M[1]) or_return;
+ redux(res, P, rho) or_return;
+ }
+ }
+ }
+
+ if redmode == 0 {
+ /*
+ Fixup result if Montgomery reduction is used.
+ Recall that any value in a Montgomery system is actually multiplied by R mod n.
+ So we have to reduce one more time to cancel out the factor of R.
+ */
+ redux(res, P, rho) or_return;
+ }
+
+ return nil;
+}
+
+/*
hac 14.61, pp608
*/
_private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
|