big: Add inverse mod.

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`.

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() {

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 */
+	}
+}

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.

core/math/big/tune.odin

@@ -23,6 +23,7 @@ Category :: enum {
 	ctz,
 	sqr,
 	bitfield_extract,
+	rm_trials,
 };
 
 Event :: struct {