math_big

We use u128 as an intermediary.

Errors are a strict superset of runtime.Allocation_Error.

======================= END OF TUNABLES =======================

Base 10 extraction constants

internal_int_is_prime switchables.

Use Frobenius-Underwood for primality testing, or use Lucas-Selfridge (default).

These defaults were tuned on an AMD A8-6600K (64-bit) using libTomMath's make tune.

TODO(Jeroen): Port this tuning algorithm and tune them for more modern processors.

	It would also be cool if we collected some data across various processor families.
	This would let uss set reasonable defaults at runtime as this library initializes
	itself by using `cpuid` or the ARM equivalent.

floor of positive solution of (2^16) - 1 = (a + 4) (2 a + 5)

• TODO: Both values are smaller than N^(1/4), would have to use a bigint
• for a instead, but any a bigger than about 120 are already so rare that
• it is possible to ignore them and still get enough pseudoprimes.
• But it is still a restriction of the set of available pseudoprimes
• which makes this implementation less secure if used stand-alone.

Trade a smaller memory footprint for more processing overhead?

Maximum number of digits.

• Must be small enough such that _bit_count does not overflow.
• Must be small enough such that _radix_size for base 2 does not overflow.

`_radix_size` needs two additional bytes for zero termination and sign.

How do we store the Ints?

Minimum number of available digits in `Int`, `_DEFAULT_DIGIT_COUNT` >= `_MIN_DIGIT_COUNT`
	- Must be at least 3 for `_div_school`.
	- Must be large enough such that `init_integer` can store `u128` in the `Int` without growing.

Overestimate the size needed for the bigint to string conversion by a very small amount.

The error is about 10^-8; it will overestimate the result by at most 11 elements for
	a number of the size 2^(2^31)-1 which is currently the largest possible in this library.
	Some short tests gave no results larger than 5 (plus 2 for sign and EOS).
Table of {0, INT(log_2([1..64])*2^p)+1 } where p is the scale
	factor defined in MP_RADIX_SIZE_SCALE and INT() extracts the integer part (truncating).
	Good for 32 bit "int". Set MP_RADIX_SIZE_SCALE = 61 and recompute values
	for 64 bit "int".

Cutoff to switch to int_factorial_binary_split, and its max recursion level.

Largest N for which we'll compute N!

Initialize constants.

Initialize constants.

Initialize constants.

Initialize constants.

Initialize constants.

Initialize constants.

How many times we'll call internal_int_random during random prime generation before we bail out.

Set to 0 or less to try indefinitely.

TODO: Make the tunables runtime adjustable where practical.

This allows to benchmark and/or setting optimized values for a certain CPU without recompiling.

Characters used in radix conversions.

How many iterations we used for the last random prime.

Runtime tunable to use Miller-Rabin primality testing only and skip the above.

======================== End of private procedures =======================

===============================  Private tables  ===============================

	Tables used by `internal_*` and `_*`.

Stores a bignum as a ASCII string in a given radix (2..64)

The buffer must be appropriately sized. This routine doesn't check.

Old internal digit extraction procedure. We're keeping this around as ground truth for the tests.

Copies DIGITs from src to dest.

Assumes `src` and `dest` to not be `nil` and have been initialized.

Direct implementation of algorithms 1.8 "RecursiveDivRem" and 1.9 "UnbalancedDivision" from:

	Brent, Richard P., and Paul Zimmermann. "Modern computer arithmetic"
		Vol. 18. Cambridge University Press, 2010
		Available online at https://arxiv.org/pdf/1004.4710

	pages 19ff. in the above online document.

Determines if a number is a valid DR modulus.

Assumes `a` to not be `nil` and to have been initialized.

Divide by three (based on routine from MPI and the GMP manual).

Signed Integer Division

c*b + d == a [i.e. a/b, c=quotient, d=remainder], HAC pp.598 Algorithm 14.20

	Note that the description in HAC is horribly incomplete.
	For example, it doesn't consider the case where digits are removed from 'x' in
	the inner loop.

	It also doesn't consider the case that y has fewer than three digits, etc.
	The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.

Slower bit-bang division... also smaller.

Prefer `_int_div_school` for speed.

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.

Determines the setup value.

Assumes `a` to not be `nil` and to have been initialized.

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.

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.

Binary split factorial algo due to: http://www.luschny.de/math/factorial/binarysplitfact.html

Internal function computing both GCD using the binary method,

	and, if target isn't `nil`, also LCM.

	Expects the `a` and `b` to have been initialized
		and one or both of `res_gcd` or `res_lcm` not to be `nil`.

	If both `a` and `b` are zero, return zero.
	If either `a` or `b`, return the other one.

	The `gcd` and `lcm` wrappers have already done this test,
	but `gcd_lcm` wouldn't have, so we still need to perform it.

	If neither result is wanted, we have nothing to do.

Internal implementation of log.

Assumes `a` not to be `nil` and to have been initialized.

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

Computes xR**-1 == x (mod N) via Montgomery Reduction.

Assumes `x` and `n` not to be nil.

Sets up the Montgomery reduction stuff.

Multiplies |a| * |b| and only computes upto digs digits of result.

HAC pp. 595, Algorithm 14.12  Modified so you can control how
	many digits of output are created.

Single-digit multiplication with the smaller number as the single-digit.

Fast (comba) multiplier

This is the fast column-array [comba] multiplier.  It is
	designed to compute the columns of the product first
	then handle the carries afterwards.  This has the effect
	of making the nested loops that compute the columns very
	simple and schedulable on super-scalar processors.

	This has been modified to produce a variable number of
	digits of output so if say only a half-product is required
	you don't have to compute the upper half (a feature
	required for fast Barrett reduction).

	Based on Algorithm 14.12 on pp.595 of HAC.

Multiplies |a| * |b| and does not compute the lower digs digits

[meant to get the higher part of the product]

This is a modified version of _private_int_mul_comba that only produces output digits above digits.

See the comments for `_private_int_mul_comba` to see how it works.

	This is used in the Barrett reduction since for one of the multiplications
	only the higher digits were needed.  This essentially halves the work.

	Based on Algorithm 14.12 on pp.595 of HAC.

product = |a| * |b| using Karatsuba Multiplication using three half size multiplications.

Let `B` represent the radix [e.g. 2**_DIGIT_BITS] and let `n` represent
	half of the number of digits in the min(a,b)

	`a` = `a1` * `B`**`n` + `a0`
	`b` = `b`1 * `B`**`n` + `b0`

	Then, a * b => 1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0

	Note that a1b1 and a0b0 are used twice and only need to be computed once.
	So in total three half size (half # of digit) multiplications are performed,
		a0b0, a1b1 and (a1+b1)(a0+b0)

	Note that a multiplication of half the digits requires 1/4th the number of
	single precision multiplications, so in total after one call 25% of the
	single precision multiplications are saved.

	Note also that the call to `internal_mul` can end up back in this function
	if the a0, a1, b0, or b1 are above the threshold.

	This is known as divide-and-conquer and leads to the famous O(N**lg(3)) or O(N**1.584)
	work which is asymptopically lower than the standard O(N**2) that the
	baseline/comba methods use. Generally though, the overhead of this method doesn't pay off
	until a certain size is reached, of around 80 used DIGITs.

Multiplication using the Toom-Cook 3-way algorithm.

Much more complicated than Karatsuba but has a lower asymptotic running time of O(N**1.464).
	This algorithm is only particularly useful on VERY large inputs.
	(We're talking 1000s of digits here...).

	This file contains code from J. Arndt's book  "Matters Computational"
	and the accompanying FXT-library with permission of the author.

	Setup from:
		Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
		18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.

	The interpolation from above needed one temporary variable more than the interpolation here:

		Bodrato, Marco, and Alberto Zanoni. "What about Toom-Cook matrices optimality."
		Centro Vito Volterra Universita di Roma Tor Vergata (2006)

Recursive product used by binary split factorial algorithm.

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.

Reduces a modulo n, where n is of the form 2**p - d.

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.

Determines the setup value.

Assumes `a` is not `nil`.

Determines the setup value.

Assumes `mu` and `P` are not `nil`.

	d := (1 << a.bits) - a;

Determines if internal_int_reduce_2k can be used.

Asssumes `a` not to be `nil` and to have been initialized.

Determines if internal_int_reduce_2k_l can be used.

Asssumes `a` not to be `nil` and to have been initialized.

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

Shift left by digits * _DIGIT_BITS bits.

Shift right by digits * _DIGIT_BITS bits.

Low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16

Assumes `dest` and `src` to not be `nil`, and `src` to have been initialized.

The jist of squaring...

You do like mult except the offset of the tmpx [one that starts closer to zero] can't equal the offset of tmpy.
	So basically you set up iy like before then you min it with (ty-tx) so that it never happens.
	You double all those you add in the inner loop. After that loop you do the squares and add them in.

	Assumes `dest` and `src` not to be `nil` and `src` to have been initialized.

Karatsuba squaring, computes dest = src * src using three half-size squarings.

See comments of `_private_int_mul_karatsuba` for details.
	It is essentially the same algorithm but merely tuned to perform recursive squarings.

Squaring using Toom-Cook 3-way algorithm.

Setup and interpolation from algorithm SQR_3 in Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
	  18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.

hac 14.61, pp608

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.

Returns the log2 of an Int.

Assumes `a` not to be `nil` and to have been initialized.
	Also assumes `base` is a power of two.

Computes xR**-1 == x (mod N) via Montgomery Reduction.

This is an optimized implementation of `internal_montgomery_reduce`
	which uses the comba method to quickly calculate the columns of the reduction.
	Based on Algorithm 14.32 on pp.601 of HAC.

Internal helpers.

Trim unused digits.

This is used to ensure that leading zero digits are trimmed and the leading "used" digit will be non-zero.
	Typically very fast.  Also fixes the sign if there are no more leading digits.

With n items, calculate how many ways that r of them can be chosen.

Also known as the multiset coefficient or (n multichoose k).

With n items, calculate how many ways that r of them can be chosen without any repeats.

Also known as the binomial coefficient or (n choose k).

Count bits in an Int.

Set dest to |src|.

===========================

	User-level routines    
	===========================
High-level addition. Handles sign.

Adds the unsigned DIGIT immediate to an Int,

such that the `DIGIT` doesn't have to be turned into an `Int` first.

	dest = a + digit;

remainder = (number + addend) % modulus.

Read a string [ASCII] in a given radix.

Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.

Made available under Odin's license.

	An arbitrary precision mathematics implementation in Odin.
	For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
	The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.

	This file contains logical operations like `and`, `or` and `xor`.
The `and`, `or` and `xor` binops differ in two lines only.
	We could handle those with a switch, but that adds overhead.

	TODO: Implement versions that take a DIGIT immediate.
2's complement `and`, returns `dest = a & b;`

dest = ~src

2's complement or, returns dest = a | b;

2's complement xor, returns dest = a ^ b;

Helpers to extract values from the Int.

Number of ways to choose k items from n items.

Also known as the binomial coefficient.

	TODO: Speed up.

	Could be done faster by reusing code from factorial and reusing the common "prefix" results for n!, k! and n-k!
	We know that n >= k, otherwise we early out with res = 0.

	So:
		n-k, keep result
		n, start from previous result
		k, start from previous result

Clear Int and resize it to the default size.

Compare two Ints, signed.

Compare an Int to an unsigned number upto the size of the backing type.

Compare the magnitude of two Ints, unsigned.

Copy one Int to another.

Returns the number of trailing zeroes before the first one.

Differs from regular `ctz` in that 0 returns 0.

TODO: Int.flags and Constants like ONE, NAN, etc, are not yet properly handled everywhere. Deallocates the backing memory of one or more Ints.

divmod.

Both the quotient and remainder are optional and may be passed a nil.

dest = src * 2

dest = src << 1

bool := a == b

bool := |a| == |b| Compares the magnitudes only, ignores the sign.

bool := a == b

Read Int from a Big Endian binary representation.

Sign is detected from the first byte if `signed` is true.

Read Int from a Big Endian Python binary representation.

Sign is detected from the first byte if `signed` is true.

Read Int from a Little Endian binary representation.

Sign is detected from the last byte if `signed` is true.

Read Int from a Little Endian Python binary representation.

Sign is detected from the first byte if `signed` is true.