From f72a0de0740d1f6146251ab3e3b7fa55b180826f Mon Sep 17 00:00:00 2001 From: Jeroen van Rijn Date: Fri, 13 Aug 2021 01:41:33 +0200 Subject: [PATCH] big: Add inverse mod. --- core/math/big/api.odin | 2 - core/math/big/example.odin | 20 +-- core/math/big/prime.odin | 45 ++++++ core/math/big/private.odin | 294 +++++++++++++++++++++++++++++++++++++ core/math/big/tune.odin | 1 + 5 files changed, 347 insertions(+), 15 deletions(-) diff --git a/core/math/big/api.odin b/core/math/big/api.odin index e1b687c14..1f2eab8d7 100644 --- a/core/math/big/api.odin +++ b/core/math/big/api.odin @@ -9,9 +9,7 @@ package math_big The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks. This file collects public proc maps and their aliases. -/* -*/ === === === === === === === === === === === === === === === === === === === === === === === === Basic arithmetic. See `public.odin`. diff --git a/core/math/big/example.odin b/core/math/big/example.odin index ae900d0c9..ff4d71059 100644 --- a/core/math/big/example.odin +++ b/core/math/big/example.odin @@ -206,19 +206,13 @@ demo :: proc() { a, b, c, d, e, f := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}; defer destroy(a, b, c, d, e, f); - set(a, 64336); - fmt.println("--- --- --- ---"); - int_to_byte(a); - fmt.println("--- --- --- ---"); - int_to_byte_little(a); - fmt.println("--- --- --- ---"); - - set(b, -64336); - fmt.println("--- --- --- ---"); - int_to_byte(b); - fmt.println("--- --- --- ---"); - int_to_byte_little(b); - fmt.println("--- --- --- ---"); + { + SCOPED_TIMING(.rm_trials); + for bits in 0..10242 { + _ = number_of_rabin_miller_trials(bits); + } + } + SCOPED_COUNT_ADD(.rm_trials, 10242); } main :: proc() { diff --git a/core/math/big/prime.odin b/core/math/big/prime.odin index 388ac3f98..81fa2f69b 100644 --- a/core/math/big/prime.odin +++ b/core/math/big/prime.odin @@ -31,3 +31,48 @@ int_prime_is_divisible :: proc(a: ^Int, allocator := context.allocator) -> (res: */ return false, nil; } + +number_of_rabin_miller_trials :: proc(bit_size: int) -> (number_of_trials: int) { + switch { + case bit_size <= 80: + return - 1; /* Use deterministic algorithm for size <= 80 bits */ + case bit_size >= 81 && bit_size < 96: + return 37; /* max. error = 2^(-96) */ + case bit_size >= 96 && bit_size < 128: + return 32; /* max. error = 2^(-96) */ + case bit_size >= 128 && bit_size < 160: + return 40; /* max. error = 2^(-112) */ + case bit_size >= 160 && bit_size < 256: + return 35; /* max. error = 2^(-112) */ + case bit_size >= 256 && bit_size < 384: + return 27; /* max. error = 2^(-128) */ + case bit_size >= 384 && bit_size < 512: + return 16; /* max. error = 2^(-128) */ + case bit_size >= 512 && bit_size < 768: + return 18; /* max. error = 2^(-160) */ + case bit_size >= 768 && bit_size < 896: + return 11; /* max. error = 2^(-160) */ + case bit_size >= 896 && bit_size < 1_024: + return 10; /* max. error = 2^(-160) */ + case bit_size >= 1_024 && bit_size < 1_536: + return 12; /* max. error = 2^(-192) */ + case bit_size >= 1_536 && bit_size < 2_048: + return 8; /* max. error = 2^(-192) */ + case bit_size >= 2_048 && bit_size < 3_072: + return 6; /* max. error = 2^(-192) */ + case bit_size >= 3_072 && bit_size < 4_096: + return 4; /* max. error = 2^(-192) */ + case bit_size >= 4_096 && bit_size < 5_120: + return 5; /* max. error = 2^(-256) */ + case bit_size >= 5_120 && bit_size < 6_144: + return 4; /* max. error = 2^(-256) */ + case bit_size >= 6_144 && bit_size < 8_192: + return 4; /* max. error = 2^(-256) */ + case bit_size >= 8_192 && bit_size < 9_216: + return 3; /* max. error = 2^(-256) */ + case bit_size >= 9_216 && bit_size < 10_240: + return 3; /* max. error = 2^(-256) */ + case: + return 2; /* For keysizes bigger than 10_240 use always at least 2 Rounds */ + } +} \ No newline at end of file diff --git a/core/math/big/private.odin b/core/math/big/private.odin index 0a5cd6163..20e9c2819 100644 --- a/core/math/big/private.odin +++ b/core/math/big/private.odin @@ -1100,6 +1100,300 @@ _private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) - return; } + +/* + hac 14.61, pp608 +*/ +_private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) { + context.allocator = allocator; + x, y, u, v, A, B, C, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}; + defer destroy(x, y, u, v, A, B, C, D); + + /* + `b` cannot be negative. + */ + if b.sign == .Negative || internal_is_zero(b) { return .Invalid_Argument; } + + /* + init temps. + */ + if err = internal_init_multi(x, y, u, v, A, B, C, D); err != nil { return err; } + + /* + `x` = `a` % `b`, `y` = `b` + */ + if err = internal_mod(x, a, b); err != nil { return err; } + if err = internal_copy(y, b); err != nil { return err; } + + /* + 2. [modified] if x,y are both even then return an error! + */ + if internal_is_even(x) && internal_is_even(y) { return .Invalid_Argument; } + + /* + 3. u=x, v=y, A=1, B=0, C=0, D=1 + */ + if err = internal_copy(u, x); err != nil { return err; } + if err = internal_copy(v, y); err != nil { return err; } + if err = internal_one(A); err != nil { return err; } + if err = internal_one(D); err != nil { return err; } + + for { + /* + 4. while `u` is even do: + */ + for internal_is_even(u) { + /* + 4.1 `u` = `u` / 2 + */ + if err = internal_int_shr1(u, u); err != nil { return err; } + + /* + 4.2 if `A` or `B` is odd then: + */ + if internal_is_odd(A) || internal_is_odd(B) { + /* + `A` = (`A`+`y`) / 2, `B` = (`B`-`x`) / 2 + */ + if err = internal_add(A, A, y); err != nil { return err; } + if err = internal_add(B, B, x); err != nil { return err; } + } + /* + `A` = `A` / 2, `B` = `B` / 2 + */ + if err = internal_int_shr1(A, A); err != nil { return err; } + if err = internal_int_shr1(B, B); err != nil { return err; } + } + + /* + 5. while `v` is even do: + */ + for internal_is_even(v) { + /* + 5.1 `v` = `v` / 2 + */ + if err = internal_int_shr1(v, v); err != nil { return err; } + + /* + 5.2 if `C` or `D` is odd then: + */ + if internal_is_odd(C) || internal_is_odd(D) { + /* + `C` = (`C`+`y`) / 2, `D` = (`D`-`x`) / 2 + */ + if err = internal_add(C, C, y); err != nil { return err; } + if err = internal_add(D, D, x); err != nil { return err; } + } + /* + `C` = `C` / 2, `D` = `D` / 2 + */ + if err = internal_int_shr1(C, C); err != nil { return err; } + if err = internal_int_shr1(D, D); err != nil { return err; } + } + + /* + 6. if `u` >= `v` then: + */ + if internal_cmp(u, v) != -1 { + /* + `u` = `u` - `v`, `A` = `A` - `C`, `B` = `B` - `D` + */ + if err = internal_sub(u, u, v); err != nil { return err; } + if err = internal_sub(A, A, C); err != nil { return err; } + if err = internal_sub(B, B, D); err != nil { return err; } + } else { + /* v - v - u, C = C - A, D = D - B */ + if err = internal_sub(v, v, u); err != nil { return err; } + if err = internal_sub(C, C, A); err != nil { return err; } + if err = internal_sub(D, D, B); err != nil { return err; } + } + + /* + If not zero goto step 4 + */ + if internal_is_zero(u) { break; } + } + + /* + Now `a` = `C`, `b` = `D`, `gcd` == `g`*`v` + */ + + /* + If `v` != `1` then there is no inverse. + */ + if internal_cmp(v, 1) != 0 { return .Invalid_Argument; } + + /* + If its too low. + */ + if internal_cmp(C, 0) == -1 { + if err = internal_add(C, C, b); err != nil { return err; } + } + + /* + Too big. + */ + if internal_cmp(C, 0) != -1 { + if err = internal_sub(C, C, b); err != nil { return err; } + } + + /* + `C` is now the inverse. + */ + swap(dest, C); + + return; +} + +/* + Computes the modular inverse via binary extended Euclidean algorithm, that is `dest` = 1 / `a` mod `b`. + + Based on slow invmod except this is optimized for the case where `b` is odd, + as per HAC Note 14.64 on pp. 610. +*/ +_private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) { + context.allocator = allocator; + x, y, u, v, B, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}; + defer destroy(x, y, u, v, B, D); + + sign: Sign; + + /* + 2. [modified] `b` must be odd. + */ + if internal_is_even(b) { return .Invalid_Argument; } + + /* + Init all our temps. + */ + if err = internal_init_multi(x, y, u, v, B, D); err != nil { return err; } + + /* + `x` == modulus, `y` == value to invert. + */ + if err = internal_copy(x, b); err != nil { return err; } + + /* + We need `y` = `|a|`. + */ + if err = internal_mod(y, a, b); err != nil { return err; } + + /* + If one of `x`, `y` is zero return an error! + */ + if internal_is_zero(x) || internal_is_zero(y) { return .Invalid_Argument; } + + /* + 3. `u` = `x`, `v` = `y`, `A` = 1, `B` = 0, `C` = 0, `D` = 1 + */ + if err = internal_copy(u, x); err != nil { return err; } + if err = internal_copy(v, y); err != nil { return err; } + + if err = internal_one(D); err != nil { return err; } + + for { + /* + 4. while `u` is even do. + */ + for internal_is_even(u) { + /* + 4.1 `u` = `u` / 2 + */ + if err = internal_int_shr1(u, u); err != nil { return err; } + + /* + 4.2 if `B` is odd then: + */ + if internal_is_odd(B) { + /* + `B` = (`B` - `x`) / 2 + */ + if err = internal_sub(B, B, x); err != nil { return err; } + } + + /* + `B` = `B` / 2 + */ + if err = internal_int_shr1(B, B); err != nil { return err; } + } + + /* + 5. while `v` is even do: + */ + for internal_is_even(v) { + /* + 5.1 `v` = `v` / 2 + */ + if err = internal_int_shr1(v, v); err != nil { return err; } + + /* + 5.2 if `D` is odd then: + */ + if internal_is_odd(D) { + /* + `D` = (`D` - `x`) / 2 + */ + if err = internal_sub(D, D, x); err != nil { return err; } + } + /* + `D` = `D` / 2 + */ + if err = internal_int_shr1(D, D); err != nil { return err; } + } + + /* + 6. if `u` >= `v` then: + */ + if internal_cmp(u, v) != -1 { + /* + `u` = `u` - `v`, `B` = `B` - `D` + */ + if err = internal_sub(u, u, v); err != nil { return err; } + if err = internal_sub(B, B, D); err != nil { return err; } + } else { + /* + `v` - `v` - `u`, `D` = `D` - `B` + */ + if err = internal_sub(v, v, u); err != nil { return err; } + if err = internal_sub(D, D, B); err != nil { return err; } + } + + /* + If not zero goto step 4. + */ + if internal_is_zero(u) { break; } + } + + /* + Now `a` = C, `b` = D, gcd == g*v + */ + + /* + if `v` != 1 then there is no inverse + */ + if internal_cmp(v, 1) != 0 { return .Invalid_Argument; } + + /* + `b` is now the inverse. + */ + sign = a.sign; + for internal_int_is_negative(D) { + if err = internal_add(D, D, b); err != nil { return err; } + } + + /* + Too big. + */ + for internal_cmp_mag(D, b) != -1 { + if err = internal_sub(D, D, b); err != nil { return err; } + } + + swap(dest, D); + dest.sign = sign; + return nil; +} + + /* Returns the log2 of an `Int`. Assumes `a` not to be `nil` and to have been initialized. diff --git a/core/math/big/tune.odin b/core/math/big/tune.odin index a51c04a78..14e10e483 100644 --- a/core/math/big/tune.odin +++ b/core/math/big/tune.odin @@ -23,6 +23,7 @@ Category :: enum { ctz, sqr, bitfield_extract, + rm_trials, }; Event :: struct {