diff options
| author | gingerBill <bill@gingerbill.org> | 2021-08-31 22:21:13 +0100 |
|---|---|---|
| committer | gingerBill <bill@gingerbill.org> | 2021-08-31 22:21:13 +0100 |
| commit | 251da264ed6e0f039931683c7b0d4b97e88c8d99 (patch) | |
| tree | c7a9a088477d2452c2cf850458c62d994a211df6 /core/math/big/private.odin | |
| parent | b176af27427a6c39448a71a8023e4a9877f0a51c (diff) | |
Remove unneeded semicolons from the core library
Diffstat (limited to 'core/math/big/private.odin')
| -rw-r--r-- | core/math/big/private.odin | 1294 |
1 files changed, 647 insertions, 647 deletions
diff --git a/core/math/big/private.odin b/core/math/big/private.odin index d71946ce2..81bae57f0 100644 --- a/core/math/big/private.odin +++ b/core/math/big/private.odin @@ -27,35 +27,35 @@ import "core:mem" many digits of output are created.
*/
_private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
/*
Can we use the fast multiplier?
*/
if digits < _WARRAY && min(a.used, b.used) < _MAX_COMBA {
- return #force_inline _private_int_mul_comba(dest, a, b, digits);
+ return #force_inline _private_int_mul_comba(dest, a, b, digits)
}
/*
Set up temporary output `Int`, which we'll swap for `dest` when done.
*/
- t := &Int{};
+ t := &Int{}
- internal_grow(t, max(digits, _DEFAULT_DIGIT_COUNT)) or_return;
- t.used = digits;
+ internal_grow(t, max(digits, _DEFAULT_DIGIT_COUNT)) or_return
+ t.used = digits
/*
Compute the digits of the product directly.
*/
- pa := a.used;
+ pa := a.used
for ix := 0; ix < pa; ix += 1 {
/*
Limit ourselves to `digits` DIGITs of output.
*/
- pb := min(b.used, digits - ix);
- carry := _WORD(0);
- iy := 0;
+ pb := min(b.used, digits - ix)
+ carry := _WORD(0)
+ iy := 0
/*
Compute the column of the output and propagate the carry.
@@ -64,29 +64,29 @@ _private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.all /*
Compute the column as a _WORD.
*/
- column := _WORD(t.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + carry;
+ column := _WORD(t.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + carry
/*
The new column is the lower part of the result.
*/
- t.digit[ix + iy] = DIGIT(column & _WORD(_MASK));
+ t.digit[ix + iy] = DIGIT(column & _WORD(_MASK))
/*
Get the carry word from the result.
*/
- carry = column >> _DIGIT_BITS;
+ carry = column >> _DIGIT_BITS
}
/*
Set carry if it is placed below digits
*/
if ix + iy < digits {
- t.digit[ix + pb] = DIGIT(carry);
+ t.digit[ix + pb] = DIGIT(carry)
}
}
- internal_swap(dest, t);
- internal_destroy(t);
- return internal_clamp(dest);
+ internal_swap(dest, t)
+ internal_destroy(t)
+ return internal_clamp(dest)
}
@@ -110,121 +110,121 @@ _private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.all Centro Vito Volterra Universita di Roma Tor Vergata (2006)
*/
_private_int_mul_toom :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- S1, S2, T1, a0, a1, a2, b0, b1, b2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(S1, S2, T1, a0, a1, a2, b0, b1, b2);
+ S1, S2, T1, a0, a1, a2, b0, b1, b2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(S1, S2, T1, a0, a1, a2, b0, b1, b2)
/*
Init temps.
*/
- internal_init_multi(S1, S2, T1) or_return;
+ internal_init_multi(S1, S2, T1) or_return
/*
B
*/
- B := min(a.used, b.used) / 3;
+ B := min(a.used, b.used) / 3
/*
a = a2 * x^2 + a1 * x + a0;
*/
- internal_grow(a0, B) or_return;
- internal_grow(a1, B) or_return;
- internal_grow(a2, a.used - 2 * B) or_return;
+ internal_grow(a0, B) or_return
+ internal_grow(a1, B) or_return
+ internal_grow(a2, a.used - 2 * B) or_return
- a0.used, a1.used = B, B;
- a2.used = a.used - 2 * B;
+ a0.used, a1.used = B, B
+ a2.used = a.used - 2 * B
- internal_copy_digits(a0, a, a0.used) or_return;
- internal_copy_digits(a1, a, a1.used, B) or_return;
- internal_copy_digits(a2, a, a2.used, 2 * B) or_return;
+ internal_copy_digits(a0, a, a0.used) or_return
+ internal_copy_digits(a1, a, a1.used, B) or_return
+ internal_copy_digits(a2, a, a2.used, 2 * B) or_return
- internal_clamp(a0);
- internal_clamp(a1);
- internal_clamp(a2);
+ internal_clamp(a0)
+ internal_clamp(a1)
+ internal_clamp(a2)
/*
b = b2 * x^2 + b1 * x + b0;
*/
- internal_grow(b0, B) or_return;
- internal_grow(b1, B) or_return;
- internal_grow(b2, b.used - 2 * B) or_return;
+ internal_grow(b0, B) or_return
+ internal_grow(b1, B) or_return
+ internal_grow(b2, b.used - 2 * B) or_return
- b0.used, b1.used = B, B;
- b2.used = b.used - 2 * B;
+ b0.used, b1.used = B, B
+ b2.used = b.used - 2 * B
- internal_copy_digits(b0, b, b0.used) or_return;
- internal_copy_digits(b1, b, b1.used, B) or_return;
- internal_copy_digits(b2, b, b2.used, 2 * B) or_return;
+ internal_copy_digits(b0, b, b0.used) or_return
+ internal_copy_digits(b1, b, b1.used, B) or_return
+ internal_copy_digits(b2, b, b2.used, 2 * B) or_return
- internal_clamp(b0);
- internal_clamp(b1);
- internal_clamp(b2);
+ internal_clamp(b0)
+ internal_clamp(b1)
+ internal_clamp(b2)
/*
\\ S1 = (a2+a1+a0) * (b2+b1+b0);
*/
- internal_add(T1, a2, a1) or_return; /* T1 = a2 + a1; */
- internal_add(S2, T1, a0) or_return; /* S2 = T1 + a0; */
- internal_add(dest, b2, b1) or_return; /* dest = b2 + b1; */
- internal_add(S1, dest, b0) or_return; /* S1 = c + b0; */
- internal_mul(S1, S1, S2) or_return; /* S1 = S1 * S2; */
+ internal_add(T1, a2, a1) or_return /* T1 = a2 + a1; */
+ internal_add(S2, T1, a0) or_return /* S2 = T1 + a0; */
+ internal_add(dest, b2, b1) or_return /* dest = b2 + b1; */
+ internal_add(S1, dest, b0) or_return /* S1 = c + b0; */
+ internal_mul(S1, S1, S2) or_return /* S1 = S1 * S2; */
/*
\\S2 = (4*a2+2*a1+a0) * (4*b2+2*b1+b0);
*/
- internal_add(T1, T1, a2) or_return; /* T1 = T1 + a2; */
- internal_int_shl1(T1, T1) or_return; /* T1 = T1 << 1; */
- internal_add(T1, T1, a0) or_return; /* T1 = T1 + a0; */
- internal_add(dest, dest, b2) or_return; /* c = c + b2; */
- internal_int_shl1(dest, dest) or_return; /* c = c << 1; */
- internal_add(dest, dest, b0) or_return; /* c = c + b0; */
- internal_mul(S2, T1, dest) or_return; /* S2 = T1 * c; */
+ internal_add(T1, T1, a2) or_return /* T1 = T1 + a2; */
+ internal_int_shl1(T1, T1) or_return /* T1 = T1 << 1; */
+ internal_add(T1, T1, a0) or_return /* T1 = T1 + a0; */
+ internal_add(dest, dest, b2) or_return /* c = c + b2; */
+ internal_int_shl1(dest, dest) or_return /* c = c << 1; */
+ internal_add(dest, dest, b0) or_return /* c = c + b0; */
+ internal_mul(S2, T1, dest) or_return /* S2 = T1 * c; */
/*
\\S3 = (a2-a1+a0) * (b2-b1+b0);
*/
- internal_sub(a1, a2, a1) or_return; /* a1 = a2 - a1; */
- internal_add(a1, a1, a0) or_return; /* a1 = a1 + a0; */
- internal_sub(b1, b2, b1) or_return; /* b1 = b2 - b1; */
- internal_add(b1, b1, b0) or_return; /* b1 = b1 + b0; */
- internal_mul(a1, a1, b1) or_return; /* a1 = a1 * b1; */
- internal_mul(b1, a2, b2) or_return; /* b1 = a2 * b2; */
+ internal_sub(a1, a2, a1) or_return /* a1 = a2 - a1; */
+ internal_add(a1, a1, a0) or_return /* a1 = a1 + a0; */
+ internal_sub(b1, b2, b1) or_return /* b1 = b2 - b1; */
+ internal_add(b1, b1, b0) or_return /* b1 = b1 + b0; */
+ internal_mul(a1, a1, b1) or_return /* a1 = a1 * b1; */
+ internal_mul(b1, a2, b2) or_return /* b1 = a2 * b2; */
/*
\\S2 = (S2 - S3) / 3;
*/
- internal_sub(S2, S2, a1) or_return; /* S2 = S2 - a1; */
- _private_int_div_3(S2, S2) or_return; /* S2 = S2 / 3; \\ this is an exact division */
- internal_sub(a1, S1, a1) or_return; /* a1 = S1 - a1; */
- internal_int_shr1(a1, a1) or_return; /* a1 = a1 >> 1; */
- internal_mul(a0, a0, b0) or_return; /* a0 = a0 * b0; */
- internal_sub(S1, S1, a0) or_return; /* S1 = S1 - a0; */
- internal_sub(S2, S2, S1) or_return; /* S2 = S2 - S1; */
- internal_int_shr1(S2, S2) or_return; /* S2 = S2 >> 1; */
- internal_sub(S1, S1, a1) or_return; /* S1 = S1 - a1; */
- internal_sub(S1, S1, b1) or_return; /* S1 = S1 - b1; */
- internal_int_shl1(T1, b1) or_return; /* T1 = b1 << 1; */
- internal_sub(S2, S2, T1) or_return; /* S2 = S2 - T1; */
- internal_sub(a1, a1, S2) or_return; /* a1 = a1 - S2; */
+ internal_sub(S2, S2, a1) or_return /* S2 = S2 - a1; */
+ _private_int_div_3(S2, S2) or_return /* S2 = S2 / 3; \\ this is an exact division */
+ internal_sub(a1, S1, a1) or_return /* a1 = S1 - a1; */
+ internal_int_shr1(a1, a1) or_return /* a1 = a1 >> 1; */
+ internal_mul(a0, a0, b0) or_return /* a0 = a0 * b0; */
+ internal_sub(S1, S1, a0) or_return /* S1 = S1 - a0; */
+ internal_sub(S2, S2, S1) or_return /* S2 = S2 - S1; */
+ internal_int_shr1(S2, S2) or_return /* S2 = S2 >> 1; */
+ internal_sub(S1, S1, a1) or_return /* S1 = S1 - a1; */
+ internal_sub(S1, S1, b1) or_return /* S1 = S1 - b1; */
+ internal_int_shl1(T1, b1) or_return /* T1 = b1 << 1; */
+ internal_sub(S2, S2, T1) or_return /* S2 = S2 - T1; */
+ internal_sub(a1, a1, S2) or_return /* a1 = a1 - S2; */
/*
P = b1*x^4+ S2*x^3+ S1*x^2+ a1*x + a0;
*/
- internal_shl_digit(b1, 4 * B) or_return;
- internal_shl_digit(S2, 3 * B) or_return;
- internal_add(b1, b1, S2) or_return;
- internal_shl_digit(S1, 2 * B) or_return;
- internal_add(b1, b1, S1) or_return;
- internal_shl_digit(a1, 1 * B) or_return;
- internal_add(b1, b1, a1) or_return;
- internal_add(dest, b1, a0) or_return;
+ internal_shl_digit(b1, 4 * B) or_return
+ internal_shl_digit(S2, 3 * B) or_return
+ internal_add(b1, b1, S2) or_return
+ internal_shl_digit(S1, 2 * B) or_return
+ internal_add(b1, b1, S1) or_return
+ internal_shl_digit(a1, 1 * B) or_return
+ internal_add(b1, b1, a1) or_return
+ internal_add(dest, b1, a0) or_return
/*
a * b - P
*/
- return nil;
+ return nil
}
/*
@@ -255,76 +255,76 @@ _private_int_mul_toom :: proc(dest, a, b: ^Int, allocator := context.allocator) until a certain size is reached, of around 80 used DIGITs.
*/
_private_int_mul_karatsuba :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- x0, x1, y0, y1, t1, x0y0, x1y1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(x0, x1, y0, y1, t1, x0y0, x1y1);
+ x0, x1, y0, y1, t1, x0y0, x1y1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(x0, x1, y0, y1, t1, x0y0, x1y1)
/*
min # of digits, divided by two.
*/
- B := min(a.used, b.used) >> 1;
+ B := min(a.used, b.used) >> 1
/*
Init all the temps.
*/
- internal_grow(x0, B) or_return;
- internal_grow(x1, a.used - B) or_return;
- internal_grow(y0, B) or_return;
- internal_grow(y1, b.used - B) or_return;
- internal_grow(t1, B * 2) or_return;
- internal_grow(x0y0, B * 2) or_return;
- internal_grow(x1y1, B * 2) or_return;
+ internal_grow(x0, B) or_return
+ internal_grow(x1, a.used - B) or_return
+ internal_grow(y0, B) or_return
+ internal_grow(y1, b.used - B) or_return
+ internal_grow(t1, B * 2) or_return
+ internal_grow(x0y0, B * 2) or_return
+ internal_grow(x1y1, B * 2) or_return
/*
Now shift the digits.
*/
- x0.used, y0.used = B, B;
- x1.used = a.used - B;
- y1.used = b.used - B;
+ x0.used, y0.used = B, B
+ x1.used = a.used - B
+ y1.used = b.used - B
/*
We copy the digits directly instead of using higher level functions
since we also need to shift the digits.
*/
- internal_copy_digits(x0, a, x0.used);
- internal_copy_digits(y0, b, y0.used);
- internal_copy_digits(x1, a, x1.used, B);
- internal_copy_digits(y1, b, y1.used, B);
+ internal_copy_digits(x0, a, x0.used)
+ internal_copy_digits(y0, b, y0.used)
+ internal_copy_digits(x1, a, x1.used, B)
+ internal_copy_digits(y1, b, y1.used, B)
/*
Only need to clamp the lower words since by definition the
upper words x1/y1 must have a known number of digits.
*/
- clamp(x0);
- clamp(y0);
+ clamp(x0)
+ clamp(y0)
/*
Now calc the products x0y0 and x1y1,
after this x0 is no longer required, free temp [x0==t2]!
*/
- internal_mul(x0y0, x0, y0) or_return; /* x0y0 = x0*y0 */
- internal_mul(x1y1, x1, y1) or_return; /* x1y1 = x1*y1 */
- internal_add(t1, x1, x0) or_return; /* now calc x1+x0 and */
- internal_add(x0, y1, y0) or_return; /* t2 = y1 + y0 */
- internal_mul(t1, t1, x0) or_return; /* t1 = (x1 + x0) * (y1 + y0) */
+ internal_mul(x0y0, x0, y0) or_return /* x0y0 = x0*y0 */
+ internal_mul(x1y1, x1, y1) or_return /* x1y1 = x1*y1 */
+ internal_add(t1, x1, x0) or_return /* now calc x1+x0 and */
+ internal_add(x0, y1, y0) or_return /* t2 = y1 + y0 */
+ internal_mul(t1, t1, x0) or_return /* t1 = (x1 + x0) * (y1 + y0) */
/*
Add x0y0.
*/
- internal_add(x0, x0y0, x1y1) or_return; /* t2 = x0y0 + x1y1 */
- internal_sub(t1, t1, x0) or_return; /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
+ internal_add(x0, x0y0, x1y1) or_return /* t2 = x0y0 + x1y1 */
+ internal_sub(t1, t1, x0) or_return /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
/*
shift by B.
*/
- internal_shl_digit(t1, B) or_return; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
- internal_shl_digit(x1y1, B * 2) or_return; /* x1y1 = x1y1 << 2*B */
+ internal_shl_digit(t1, B) or_return /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
+ internal_shl_digit(x1y1, B * 2) or_return /* x1y1 = x1y1 << 2*B */
- internal_add(t1, x0y0, t1) or_return; /* t1 = x0y0 + t1 */
- internal_add(dest, t1, x1y1) or_return; /* t1 = x0y0 + t1 + x1y1 */
+ internal_add(t1, x0y0, t1) or_return /* t1 = x0y0 + t1 */
+ internal_add(dest, t1, x1y1) or_return /* t1 = x0y0 + t1 + x1y1 */
- return nil;
+ return nil
}
@@ -346,84 +346,84 @@ _private_int_mul_karatsuba :: proc(dest, a, b: ^Int, allocator := context.alloca Based on Algorithm 14.12 on pp.595 of HAC.
*/
_private_int_mul_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
/*
Set up array.
*/
- W: [_WARRAY]DIGIT = ---;
+ W: [_WARRAY]DIGIT = ---
/*
Grow the destination as required.
*/
- internal_grow(dest, digits) or_return;
+ internal_grow(dest, digits) or_return
/*
Number of output digits to produce.
*/
- pa := min(digits, a.used + b.used);
+ pa := min(digits, a.used + b.used)
/*
Clear the carry
*/
- _W := _WORD(0);
+ _W := _WORD(0)
- ix: int;
+ ix: int
for ix = 0; ix < pa; ix += 1 {
- tx, ty, iy, iz: int;
+ tx, ty, iy, iz: int
/*
Get offsets into the two bignums.
*/
- ty = min(b.used - 1, ix);
- tx = ix - ty;
+ ty = min(b.used - 1, ix)
+ tx = ix - ty
/*
This is the number of times the loop will iterate, essentially.
while (tx++ < a->used && ty-- >= 0) { ... }
*/
- iy = min(a.used - tx, ty + 1);
+ iy = min(a.used - tx, ty + 1)
/*
Execute loop.
*/
#no_bounds_check for iz = 0; iz < iy; iz += 1 {
- _W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz]);
+ _W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])
}
/*
Store term.
*/
- W[ix] = DIGIT(_W) & _MASK;
+ W[ix] = DIGIT(_W) & _MASK
/*
Make next carry.
*/
- _W = _W >> _WORD(_DIGIT_BITS);
+ _W = _W >> _WORD(_DIGIT_BITS)
}
/*
Setup dest.
*/
- old_used := dest.used;
- dest.used = pa;
+ old_used := dest.used
+ dest.used = pa
/*
Now extract the previous digit [below the carry].
*/
- copy_slice(dest.digit[0:], W[:pa]);
+ copy_slice(dest.digit[0:], W[:pa])
/*
Clear unused digits [that existed in the old copy of dest].
*/
- internal_zero_unused(dest, old_used);
+ internal_zero_unused(dest, old_used)
/*
Adjust dest.used based on leading zeroes.
*/
- return internal_clamp(dest);
+ return internal_clamp(dest)
}
/*
@@ -431,42 +431,42 @@ _private_int_mul_comba :: proc(dest, a, b: ^Int, digits: int, allocator := conte [meant to get the higher part of the product]
*/
_private_int_mul_high :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
/*
Can we use the fast multiplier?
*/
if a.used + b.used + 1 < _WARRAY && min(a.used, b.used) < _MAX_COMBA {
- return _private_int_mul_high_comba(dest, a, b, digits);
+ return _private_int_mul_high_comba(dest, a, b, digits)
}
- internal_grow(dest, a.used + b.used + 1) or_return;
- dest.used = a.used + b.used + 1;
+ internal_grow(dest, a.used + b.used + 1) or_return
+ dest.used = a.used + b.used + 1
- pa := a.used;
- pb := b.used;
+ pa := a.used
+ pb := b.used
for ix := 0; ix < pa; ix += 1 {
- carry := DIGIT(0);
+ carry := DIGIT(0)
for iy := digits - ix; iy < pb; iy += 1 {
/*
Calculate the double precision result.
*/
- r := _WORD(dest.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + _WORD(carry);
+ r := _WORD(dest.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + _WORD(carry)
/*
Get the lower part.
*/
- dest.digit[ix + iy] = DIGIT(r & _WORD(_MASK));
+ dest.digit[ix + iy] = DIGIT(r & _WORD(_MASK))
/*
Carry the carry.
*/
- carry = DIGIT(r >> _WORD(_DIGIT_BITS));
+ carry = DIGIT(r >> _WORD(_DIGIT_BITS))
}
- dest.digit[ix + pb] = carry;
+ dest.digit[ix + pb] = carry
}
- return internal_clamp(dest);
+ return internal_clamp(dest)
}
/*
@@ -479,140 +479,140 @@ _private_int_mul_high :: proc(dest, a, b: ^Int, digits: int, allocator := contex Based on Algorithm 14.12 on pp.595 of HAC.
*/
_private_int_mul_high_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- W: [_WARRAY]DIGIT = ---;
- _W: _WORD = 0;
+ W: [_WARRAY]DIGIT = ---
+ _W: _WORD = 0
/*
Number of output digits to produce. Grow the destination as required.
*/
- pa := a.used + b.used;
- internal_grow(dest, pa) or_return;
+ pa := a.used + b.used
+ internal_grow(dest, pa) or_return
- ix: int;
+ ix: int
for ix = digits; ix < pa; ix += 1 {
/*
Get offsets into the two bignums.
*/
- ty := min(b.used - 1, ix);
- tx := ix - ty;
+ ty := min(b.used - 1, ix)
+ tx := ix - ty
/*
This is the number of times the loop will iterrate, essentially it's
while (tx++ < a->used && ty-- >= 0) { ... }
*/
- iy := min(a.used - tx, ty + 1);
+ iy := min(a.used - tx, ty + 1)
/*
Execute loop.
*/
for iz := 0; iz < iy; iz += 1 {
- _W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz]);
+ _W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])
}
/*
Store term.
*/
- W[ix] = DIGIT(_W) & DIGIT(_MASK);
+ W[ix] = DIGIT(_W) & DIGIT(_MASK)
/*
Make next carry.
*/
- _W = _W >> _WORD(_DIGIT_BITS);
+ _W = _W >> _WORD(_DIGIT_BITS)
}
/*
Setup dest
*/
- old_used := dest.used;
- dest.used = pa;
+ old_used := dest.used
+ dest.used = pa
for ix = digits; ix < pa; ix += 1 {
/*
Now extract the previous digit [below the carry].
*/
- dest.digit[ix] = W[ix];
+ dest.digit[ix] = W[ix]
}
/*
Zero remainder.
*/
- internal_zero_unused(dest, old_used);
+ internal_zero_unused(dest, old_used)
/*
Adjust dest.used based on leading zeroes.
*/
- return internal_clamp(dest);
+ return internal_clamp(dest)
}
/*
Single-digit multiplication with the smaller number as the single-digit.
*/
_private_int_mul_balance :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
- a, b := a, b;
+ context.allocator = allocator
+ a, b := a, b
- a0, tmp, r := &Int{}, &Int{}, &Int{};
- defer internal_destroy(a0, tmp, r);
+ a0, tmp, r := &Int{}, &Int{}, &Int{}
+ defer internal_destroy(a0, tmp, r)
- b_size := min(a.used, b.used);
- n_blocks := max(a.used, b.used) / b_size;
+ b_size := min(a.used, b.used)
+ n_blocks := max(a.used, b.used) / b_size
- internal_grow(a0, b_size + 2) or_return;
- internal_init_multi(tmp, r) or_return;
+ internal_grow(a0, b_size + 2) or_return
+ internal_init_multi(tmp, r) or_return
/*
Make sure that `a` is the larger one.
*/
if a.used < b.used {
- a, b = b, a;
+ a, b = b, a
}
- assert(a.used >= b.used);
+ assert(a.used >= b.used)
- i, j := 0, 0;
+ i, j := 0, 0
for ; i < n_blocks; i += 1 {
/*
Cut a slice off of `a`.
*/
- a0.used = b_size;
- internal_copy_digits(a0, a, a0.used, j);
- j += a0.used;
- internal_clamp(a0);
+ a0.used = b_size
+ internal_copy_digits(a0, a, a0.used, j)
+ j += a0.used
+ internal_clamp(a0)
/*
Multiply with `b`.
*/
- internal_mul(tmp, a0, b) or_return;
+ internal_mul(tmp, a0, b) or_return
/*
Shift `tmp` to the correct position.
*/
- internal_shl_digit(tmp, b_size * i) or_return;
+ internal_shl_digit(tmp, b_size * i) or_return
/*
Add to output. No carry needed.
*/
- internal_add(r, r, tmp) or_return;
+ internal_add(r, r, tmp) or_return
}
/*
The left-overs; there are always left-overs.
*/
if j < a.used {
- a0.used = a.used - j;
- internal_copy_digits(a0, a, a0.used, j);
- j += a0.used;
- internal_clamp(a0);
-
- internal_mul(tmp, a0, b) or_return;
- internal_shl_digit(tmp, b_size * i) or_return;
- internal_add(r, r, tmp) or_return;
+ a0.used = a.used - j
+ internal_copy_digits(a0, a, a0.used, j)
+ j += a0.used
+ internal_clamp(a0)
+
+ internal_mul(tmp, a0, b) or_return
+ internal_shl_digit(tmp, b_size * i) or_return
+ internal_add(r, r, tmp) or_return
}
- internal_swap(dest, r);
- return;
+ internal_swap(dest, r)
+ return
}
/*
@@ -620,68 +620,68 @@ _private_int_mul_balance :: proc(dest, a, b: ^Int, allocator := context.allocato Assumes `dest` and `src` to not be `nil`, and `src` to have been initialized.
*/
_private_int_sqr :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
- pa := src.used;
+ context.allocator = allocator
+ pa := src.used
- t := &Int{}; ix, iy: int;
+ t := &Int{}; ix, iy: int
/*
Grow `t` to maximum needed size, or `_DEFAULT_DIGIT_COUNT`, whichever is bigger.
*/
- internal_grow(t, max((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)) or_return;
- t.used = (2 * pa) + 1;
+ internal_grow(t, max((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)) or_return
+ t.used = (2 * pa) + 1
#no_bounds_check for ix = 0; ix < pa; ix += 1 {
- carry := DIGIT(0);
+ carry := DIGIT(0)
/*
First calculate the digit at 2*ix; calculate double precision result.
*/
- r := _WORD(t.digit[ix+ix]) + (_WORD(src.digit[ix]) * _WORD(src.digit[ix]));
+ r := _WORD(t.digit[ix+ix]) + (_WORD(src.digit[ix]) * _WORD(src.digit[ix]))
/*
Store lower part in result.
*/
- t.digit[ix+ix] = DIGIT(r & _WORD(_MASK));
+ t.digit[ix+ix] = DIGIT(r & _WORD(_MASK))
/*
Get the carry.
*/
- carry = DIGIT(r >> _DIGIT_BITS);
+ carry = DIGIT(r >> _DIGIT_BITS)
#no_bounds_check for iy = ix + 1; iy < pa; iy += 1 {
/*
First calculate the product.
*/
- r = _WORD(src.digit[ix]) * _WORD(src.digit[iy]);
+ r = _WORD(src.digit[ix]) * _WORD(src.digit[iy])
/* Now calculate the double precision result. Nóte we use
* addition instead of *2 since it's easier to optimize
*/
- r = _WORD(t.digit[ix+iy]) + r + r + _WORD(carry);
+ r = _WORD(t.digit[ix+iy]) + r + r + _WORD(carry)
/*
Store lower part.
*/
- t.digit[ix+iy] = DIGIT(r & _WORD(_MASK));
+ t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))
/*
Get carry.
*/
- carry = DIGIT(r >> _DIGIT_BITS);
+ carry = DIGIT(r >> _DIGIT_BITS)
}
/*
Propagate upwards.
*/
#no_bounds_check for carry != 0 {
- r = _WORD(t.digit[ix+iy]) + _WORD(carry);
- t.digit[ix+iy] = DIGIT(r & _WORD(_MASK));
- carry = DIGIT(r >> _WORD(_DIGIT_BITS));
- iy += 1;
+ r = _WORD(t.digit[ix+iy]) + _WORD(carry)
+ t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))
+ carry = DIGIT(r >> _WORD(_DIGIT_BITS))
+ iy += 1
}
}
- err = internal_clamp(t);
- internal_swap(dest, t);
- internal_destroy(t);
- return err;
+ err = internal_clamp(t)
+ internal_swap(dest, t)
+ internal_destroy(t)
+ return err
}
/*
@@ -693,94 +693,94 @@ _private_int_sqr :: proc(dest, src: ^Int, allocator := context.allocator) -> (er Assumes `dest` and `src` not to be `nil` and `src` to have been initialized.
*/
_private_int_sqr_comba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- W: [_WARRAY]DIGIT = ---;
+ W: [_WARRAY]DIGIT = ---
/*
Grow the destination as required.
*/
- pa := uint(src.used) + uint(src.used);
- internal_grow(dest, int(pa)) or_return;
+ pa := uint(src.used) + uint(src.used)
+ internal_grow(dest, int(pa)) or_return
/*
Number of output digits to produce.
*/
- W1 := _WORD(0);
- _W : _WORD = ---;
- ix := uint(0);
+ W1 := _WORD(0)
+ _W : _WORD = ---
+ ix := uint(0)
#no_bounds_check for ; ix < pa; ix += 1 {
/*
Clear counter.
*/
- _W = {};
+ _W = {}
/*
Get offsets into the two bignums.
*/
- ty := min(uint(src.used) - 1, ix);
- tx := ix - ty;
+ ty := min(uint(src.used) - 1, ix)
+ tx := ix - ty
/*
This is the number of times the loop will iterate,
essentially while (tx++ < a->used && ty-- >= 0) { ... }
*/
- iy := min(uint(src.used) - tx, ty + 1);
+ iy := min(uint(src.used) - tx, ty + 1)
/*
Now for squaring, tx can never equal ty.
We halve the distance since they approach at a rate of 2x,
and we have to round because odd cases need to be executed.
*/
- iy = min(iy, ((ty - tx) + 1) >> 1 );
+ iy = min(iy, ((ty - tx) + 1) >> 1 )
/*
Execute loop.
*/
#no_bounds_check for iz := uint(0); iz < iy; iz += 1 {
- _W += _WORD(src.digit[tx + iz]) * _WORD(src.digit[ty - iz]);
+ _W += _WORD(src.digit[tx + iz]) * _WORD(src.digit[ty - iz])
}
/*
Double the inner product and add carry.
*/
- _W = _W + _W + W1;
+ _W = _W + _W + W1
/*
Even columns have the square term in them.
*/
if ix & 1 == 0 {
- _W += _WORD(src.digit[ix >> 1]) * _WORD(src.digit[ix >> 1]);
+ _W += _WORD(src.digit[ix >> 1]) * _WORD(src.digit[ix >> 1])
}
/*
Store it.
*/
- W[ix] = DIGIT(_W & _WORD(_MASK));
+ W[ix] = DIGIT(_W & _WORD(_MASK))
/*
Make next carry.
*/
- W1 = _W >> _DIGIT_BITS;
+ W1 = _W >> _DIGIT_BITS
}
/*
Setup dest.
*/
- old_used := dest.used;
- dest.used = src.used + src.used;
+ old_used := dest.used
+ dest.used = src.used + src.used
#no_bounds_check for ix = 0; ix < pa; ix += 1 {
- dest.digit[ix] = W[ix] & _MASK;
+ dest.digit[ix] = W[ix] & _MASK
}
/*
Clear unused digits [that existed in the old copy of dest].
*/
- internal_zero_unused(dest, old_used);
+ internal_zero_unused(dest, old_used)
- return internal_clamp(dest);
+ return internal_clamp(dest)
}
/*
@@ -790,31 +790,31 @@ _private_int_sqr_comba :: proc(dest, src: ^Int, allocator := context.allocator) It is essentially the same algorithm but merely tuned to perform recursive squarings.
*/
_private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- x0, x1, t1, t2, x0x0, x1x1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(x0, x1, t1, t2, x0x0, x1x1);
+ x0, x1, t1, t2, x0x0, x1x1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(x0, x1, t1, t2, x0x0, x1x1)
/*
Min # of digits, divided by two.
*/
- B := src.used >> 1;
+ B := src.used >> 1
/*
Init temps.
*/
- internal_grow(x0, B) or_return;
- internal_grow(x1, src.used - B) or_return;
- internal_grow(t1, src.used * 2) or_return;
- internal_grow(t2, src.used * 2) or_return;
- internal_grow(x0x0, B * 2 ) or_return;
- internal_grow(x1x1, (src.used - B) * 2) or_return;
+ internal_grow(x0, B) or_return
+ internal_grow(x1, src.used - B) or_return
+ internal_grow(t1, src.used * 2) or_return
+ internal_grow(t2, src.used * 2) or_return
+ internal_grow(x0x0, B * 2 ) or_return
+ internal_grow(x1x1, (src.used - B) * 2) or_return
/*
Now shift the digits.
*/
- x0.used = B;
- x1.used = src.used - B;
+ x0.used = B
+ x1.used = src.used - B
#force_inline internal_copy_digits(x0, src, x0.used);
#force_inline mem.copy_non_overlapping(&x1.digit[0], &src.digit[B], size_of(DIGIT) * x1.used);
@@ -823,30 +823,30 @@ _private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocat /*
Now calc the products x0*x0 and x1*x1.
*/
- internal_sqr(x0x0, x0) or_return;
- internal_sqr(x1x1, x1) or_return;
+ internal_sqr(x0x0, x0) or_return
+ internal_sqr(x1x1, x1) or_return
/*
Now calc (x1+x0)^2
*/
- internal_add(t1, x0, x1) or_return;
- internal_sqr(t1, t1) or_return;
+ internal_add(t1, x0, x1) or_return
+ internal_sqr(t1, t1) or_return
/*
Add x0y0
*/
- internal_add(t2, x0x0, x1x1) or_return;
- internal_sub(t1, t1, t2) or_return;
+ internal_add(t2, x0x0, x1x1) or_return
+ internal_sub(t1, t1, t2) or_return
/*
Shift by B.
*/
- internal_shl_digit(t1, B) or_return;
- internal_shl_digit(x1x1, B * 2) or_return;
- internal_add(t1, t1, x0x0) or_return;
- internal_add(dest, t1, x1x1) or_return;
+ internal_shl_digit(t1, B) or_return
+ internal_shl_digit(x1x1, B * 2) or_return
+ internal_add(t1, t1, x0x0) or_return
+ internal_add(dest, t1, x1x1) or_return
- return #force_inline internal_clamp(dest);
+ return #force_inline internal_clamp(dest)
}
/*
@@ -856,160 +856,160 @@ _private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocat 18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.
*/
_private_int_sqr_toom :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- S0, a0, a1, a2 := &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(S0, a0, a1, a2);
+ S0, a0, a1, a2 := &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(S0, a0, a1, a2)
/*
Init temps.
*/
- internal_zero(S0) or_return;
+ internal_zero(S0) or_return
/*
B
*/
- B := src.used / 3;
+ B := src.used / 3
/*
a = a2 * x^2 + a1 * x + a0;
*/
- internal_grow(a0, B) or_return;
- internal_grow(a1, B) or_return;
- internal_grow(a2, src.used - (2 * B)) or_return;
+ internal_grow(a0, B) or_return
+ internal_grow(a1, B) or_return
+ internal_grow(a2, src.used - (2 * B)) or_return
- a0.used = B;
- a1.used = B;
- a2.used = src.used - 2 * B;
+ a0.used = B
+ a1.used = B
+ a2.used = src.used - 2 * B
#force_inline mem.copy_non_overlapping(&a0.digit[0], &src.digit[ 0], size_of(DIGIT) * a0.used);
#force_inline mem.copy_non_overlapping(&a1.digit[0], &src.digit[ B], size_of(DIGIT) * a1.used);
#force_inline mem.copy_non_overlapping(&a2.digit[0], &src.digit[2 * B], size_of(DIGIT) * a2.used);
- internal_clamp(a0);
- internal_clamp(a1);
- internal_clamp(a2);
+ internal_clamp(a0)
+ internal_clamp(a1)
+ internal_clamp(a2)
/** S0 = a0^2; */
- internal_sqr(S0, a0) or_return;
+ internal_sqr(S0, a0) or_return
/** \\S1 = (a2 + a1 + a0)^2 */
/** \\S2 = (a2 - a1 + a0)^2 */
/** \\S1 = a0 + a2; */
/** a0 = a0 + a2; */
- internal_add(a0, a0, a2) or_return;
+ internal_add(a0, a0, a2) or_return
/** \\S2 = S1 - a1; */
/** b = a0 - a1; */
- internal_sub(dest, a0, a1) or_return;
+ internal_sub(dest, a0, a1) or_return
/** \\S1 = S1 + a1; */
/** a0 = a0 + a1; */
- internal_add(a0, a0, a1) or_return;
+ internal_add(a0, a0, a1) or_return
/** \\S1 = S1^2; */
/** a0 = a0^2; */
- internal_sqr(a0, a0) or_return;
+ internal_sqr(a0, a0) or_return
/** \\S2 = S2^2; */
/** b = b^2; */
- internal_sqr(dest, dest) or_return;
+ internal_sqr(dest, dest) or_return
/** \\ S3 = 2 * a1 * a2 */
/** \\S3 = a1 * a2; */
/** a1 = a1 * a2; */
- internal_mul(a1, a1, a2) or_return;
+ internal_mul(a1, a1, a2) or_return
/** \\S3 = S3 << 1; */
/** a1 = a1 << 1; */
- internal_shl(a1, a1, 1) or_return;
+ internal_shl(a1, a1, 1) or_return
/** \\S4 = a2^2; */
/** a2 = a2^2; */
- internal_sqr(a2, a2) or_return;
+ internal_sqr(a2, a2) or_return
/** \\ tmp = (S1 + S2)/2 */
/** \\tmp = S1 + S2; */
/** b = a0 + b; */
- internal_add(dest, a0, dest) or_return;
+ internal_add(dest, a0, dest) or_return
/** \\tmp = tmp >> 1; */
/** b = b >> 1; */
- internal_shr(dest, dest, 1) or_return;
+ internal_shr(dest, dest, 1) or_return
/** \\ S1 = S1 - tmp - S3 */
/** \\S1 = S1 - tmp; */
/** a0 = a0 - b; */
- internal_sub(a0, a0, dest) or_return;
+ internal_sub(a0, a0, dest) or_return
/** \\S1 = S1 - S3; */
/** a0 = a0 - a1; */
- internal_sub(a0, a0, a1) or_return;
+ internal_sub(a0, a0, a1) or_return
/** \\S2 = tmp - S4 -S0 */
/** \\S2 = tmp - S4; */
/** b = b - a2; */
- internal_sub(dest, dest, a2) or_return;
+ internal_sub(dest, dest, a2) or_return
/** \\S2 = S2 - S0; */
/** b = b - S0; */
- internal_sub(dest, dest, S0) or_return;
+ internal_sub(dest, dest, S0) or_return
/** \\P = S4*x^4 + S3*x^3 + S2*x^2 + S1*x + S0; */
/** P = a2*x^4 + a1*x^3 + b*x^2 + a0*x + S0; */
- internal_shl_digit( a2, 4 * B) or_return;
- internal_shl_digit( a1, 3 * B) or_return;
- internal_shl_digit(dest, 2 * B) or_return;
- internal_shl_digit( a0, 1 * B) or_return;
-
- internal_add(a2, a2, a1) or_return;
- internal_add(dest, dest, a2) or_return;
- internal_add(dest, dest, a0) or_return;
- internal_add(dest, dest, S0) or_return;
+ internal_shl_digit( a2, 4 * B) or_return
+ internal_shl_digit( a1, 3 * B) or_return
+ internal_shl_digit(dest, 2 * B) or_return
+ internal_shl_digit( a0, 1 * B) or_return
+
+ internal_add(a2, a2, a1) or_return
+ internal_add(dest, dest, a2) or_return
+ internal_add(dest, dest, a0) or_return
+ internal_add(dest, dest, S0) or_return
/** a^2 - P */
- return #force_inline internal_clamp(dest);
+ return #force_inline internal_clamp(dest)
}
/*
Divide by three (based on routine from MPI and the GMP manual).
*/
_private_int_div_3 :: proc(quotient, numerator: ^Int, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
/*
b = 2^_DIGIT_BITS / 3
*/
- b := _WORD(1) << _WORD(_DIGIT_BITS) / _WORD(3);
+ b := _WORD(1) << _WORD(_DIGIT_BITS) / _WORD(3)
- q := &Int{};
- internal_grow(q, numerator.used) or_return;
- q.used = numerator.used;
- q.sign = numerator.sign;
+ q := &Int{}
+ internal_grow(q, numerator.used) or_return
+ q.used = numerator.used
+ q.sign = numerator.sign
- w, t: _WORD;
+ w, t: _WORD
#no_bounds_check for ix := numerator.used; ix >= 0; ix -= 1 {
- w = (w << _WORD(_DIGIT_BITS)) | _WORD(numerator.digit[ix]);
+ w = (w << _WORD(_DIGIT_BITS)) | _WORD(numerator.digit[ix])
if w >= 3 {
/*
Multiply w by [1/3].
*/
- t = (w * b) >> _WORD(_DIGIT_BITS);
+ t = (w * b) >> _WORD(_DIGIT_BITS)
/*
Now subtract 3 * [w/3] from w, to get the remainder.
*/
- w -= t+t+t;
+ w -= t+t+t
/*
Fixup the remainder as required since the optimization is not exact.
*/
for w >= 3 {
- t += 1;
- w -= 3;
+ t += 1
+ w -= 3
}
} else {
- t = 0;
+ t = 0
}
- q.digit[ix] = DIGIT(t);
+ q.digit[ix] = DIGIT(t)
}
- remainder = DIGIT(w);
+ remainder = DIGIT(w)
/*
[optional] store the quotient.
*/
if quotient != nil {
- err = clamp(q);
- internal_swap(q, quotient);
+ err = clamp(q)
+ internal_swap(q, quotient)
}
- internal_destroy(q);
- return remainder, nil;
+ internal_destroy(q)
+ return remainder, nil
}
/*
@@ -1025,64 +1025,64 @@ _private_int_div_3 :: proc(quotient, numerator: ^Int, allocator := context.alloc The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.
*/
_private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- error_if_immutable(quotient, remainder) or_return;
+ error_if_immutable(quotient, remainder) or_return
- q, x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(q, x, y, t1, t2);
+ q, x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(q, x, y, t1, t2)
- internal_grow(q, numerator.used + 2) or_return;
- q.used = numerator.used + 2;
+ internal_grow(q, numerator.used + 2) or_return
+ q.used = numerator.used + 2
- internal_init_multi(t1, t2) or_return;
- internal_copy(x, numerator) or_return;
- internal_copy(y, denominator) or_return;
+ internal_init_multi(t1, t2) or_return
+ internal_copy(x, numerator) or_return
+ internal_copy(y, denominator) or_return
/*
Fix the sign.
*/
- neg := numerator.sign != denominator.sign;
- x.sign = .Zero_or_Positive;
- y.sign = .Zero_or_Positive;
+ neg := numerator.sign != denominator.sign
+ x.sign = .Zero_or_Positive
+ y.sign = .Zero_or_Positive
/*
Normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT]
*/
- norm := internal_count_bits(y) % _DIGIT_BITS;
+ norm := internal_count_bits(y) % _DIGIT_BITS
if norm < _DIGIT_BITS - 1 {
- norm = (_DIGIT_BITS - 1) - norm;
- internal_shl(x, x, norm) or_return;
- internal_shl(y, y, norm) or_return;
+ norm = (_DIGIT_BITS - 1) - norm
+ internal_shl(x, x, norm) or_return
+ internal_shl(y, y, norm) or_return
} else {
- norm = 0;
+ norm = 0
}
/*
Note: HAC does 0 based, so if used==5 then it's 0,1,2,3,4, i.e. use 4
*/
- n := x.used - 1;
- t := y.used - 1;
+ n := x.used - 1
+ t := y.used - 1
/*
while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} }
y = y*b**{n-t}
*/
- internal_shl_digit(y, n - t) or_return;
+ internal_shl_digit(y, n - t) or_return
- c := internal_cmp(x, y);
+ c := internal_cmp(x, y)
for c != -1 {
- q.digit[n - t] += 1;
- internal_sub(x, x, y) or_return;
- c = internal_cmp(x, y);
+ q.digit[n - t] += 1
+ internal_sub(x, x, y) or_return
+ c = internal_cmp(x, y)
}
/*
Reset y by shifting it back down.
*/
- internal_shr_digit(y, n - t);
+ internal_shr_digit(y, n - t)
/*
Step 3. for i from n down to (t + 1).
@@ -1094,16 +1094,16 @@ _private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^In step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt
*/
if x.digit[i] == y.digit[t] {
- q.digit[(i - t) - 1] = 1 << (_DIGIT_BITS - 1);
+ q.digit[(i - t) - 1] = 1 << (_DIGIT_BITS - 1)
} else {
- tmp := _WORD(x.digit[i]) << _DIGIT_BITS;
- tmp |= _WORD(x.digit[i - 1]);
- tmp /= _WORD(y.digit[t]);
+ tmp := _WORD(x.digit[i]) << _DIGIT_BITS
+ tmp |= _WORD(x.digit[i - 1])
+ tmp /= _WORD(y.digit[t])
if tmp > _WORD(_MASK) {
- tmp = _WORD(_MASK);
+ tmp = _WORD(_MASK)
}
- q.digit[(i - t) - 1] = DIGIT(tmp & _WORD(_MASK));
+ q.digit[(i - t) - 1] = DIGIT(tmp & _WORD(_MASK))
}
/* while (q{i-t-1} * (yt * b + y{t-1})) >
@@ -1112,53 +1112,53 @@ _private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^In do q{i-t-1} -= 1;
*/
- iter := 0;
+ iter := 0
- q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] + 1) & _MASK;
+ q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] + 1) & _MASK
#no_bounds_check for {
- q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK;
+ q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK
/*
Find left hand.
*/
- internal_zero(t1);
- t1.digit[0] = ((t - 1) < 0) ? 0 : y.digit[t - 1];
- t1.digit[1] = y.digit[t];
- t1.used = 2;
- internal_mul(t1, t1, q.digit[(i - t) - 1]) or_return;
+ internal_zero(t1)
+ t1.digit[0] = ((t - 1) < 0) ? 0 : y.digit[t - 1]
+ t1.digit[1] = y.digit[t]
+ t1.used = 2
+ internal_mul(t1, t1, q.digit[(i - t) - 1]) or_return
/*
Find right hand.
*/
- t2.digit[0] = ((i - 2) < 0) ? 0 : x.digit[i - 2];
- t2.digit[1] = x.digit[i - 1]; /* i >= 1 always holds */
- t2.digit[2] = x.digit[i];
- t2.used = 3;
+ t2.digit[0] = ((i - 2) < 0) ? 0 : x.digit[i - 2]
+ t2.digit[1] = x.digit[i - 1] /* i >= 1 always holds */
+ t2.digit[2] = x.digit[i]
+ t2.used = 3
if t1_t2 := internal_cmp_mag(t1, t2); t1_t2 != 1 {
- break;
+ break
}
iter += 1; if iter > 100 {
- return .Max_Iterations_Reached;
+ return .Max_Iterations_Reached
}
}
/*
Step 3.3 x = x - q{i-t-1} * y * b**{i-t-1}
*/
- int_mul_digit(t1, y, q.digit[(i - t) - 1]) or_return;
- internal_shl_digit(t1, (i - t) - 1) or_return;
- internal_sub(x, x, t1) or_return;
+ int_mul_digit(t1, y, q.digit[(i - t) - 1]) or_return
+ internal_shl_digit(t1, (i - t) - 1) or_return
+ internal_sub(x, x, t1) or_return
/*
if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; }
*/
if x.sign == .Negative {
- internal_copy(t1, y) or_return;
- internal_shl_digit(t1, (i - t) - 1) or_return;
- internal_add(x, x, t1) or_return;
+ internal_copy(t1, y) or_return
+ internal_shl_digit(t1, (i - t) - 1) or_return
+ internal_add(x, x, t1) or_return
- q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK;
+ q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK
}
}
@@ -1166,21 +1166,21 @@ _private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^In Now q is the quotient and x is the remainder, [which we have to normalize]
Get sign before writing to c.
*/
- z, _ := is_zero(x);
- x.sign = .Zero_or_Positive if z else numerator.sign;
+ z, _ := is_zero(x)
+ x.sign = .Zero_or_Positive if z else numerator.sign
if quotient != nil {
- internal_clamp(q);
- internal_swap(q, quotient);
- quotient.sign = .Negative if neg else .Zero_or_Positive;
+ internal_clamp(q)
+ internal_swap(q, quotient)
+ quotient.sign = .Negative if neg else .Zero_or_Positive
}
if remainder != nil {
- internal_shr(x, x, norm) or_return;
- internal_swap(x, remainder);
+ internal_shr(x, x, norm) or_return
+ internal_swap(x, remainder)
}
- return nil;
+ return nil
}
/*
@@ -1193,49 +1193,49 @@ _private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^In pages 19ff. in the above online document.
*/
_private_div_recursion :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- A1, A2, B1, B0, Q1, Q0, R1, R0, t := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(A1, A2, B1, B0, Q1, Q0, R1, R0, t);
+ A1, A2, B1, B0, Q1, Q0, R1, R0, t := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(A1, A2, B1, B0, Q1, Q0, R1, R0, t)
- m := a.used - b.used;
- k := m / 2;
+ m := a.used - b.used
+ k := m / 2
if m < MUL_KARATSUBA_CUTOFF {
- return _private_int_div_school(quotient, remainder, a, b);
+ return _private_int_div_school(quotient, remainder, a, b)
}
- internal_init_multi(A1, A2, B1, B0, Q1, Q0, R1, R0, t) or_return;
+ internal_init_multi(A1, A2, B1, B0, Q1, Q0, R1, R0, t) or_return
/*
`B1` = `b` / `beta`^`k`, `B0` = `b` % `beta`^`k`
*/
- internal_shrmod(B1, B0, b, k * _DIGIT_BITS) or_return;
+ internal_shrmod(B1, B0, b, k * _DIGIT_BITS) or_return
/*
(Q1, R1) = RecursiveDivRem(A / beta^(2k), B1)
*/
- internal_shrmod(A1, t, a, 2 * k * _DIGIT_BITS) or_return;
- _private_div_recursion(Q1, R1, A1, B1) or_return;
+ internal_shrmod(A1, t, a, 2 * k * _DIGIT_BITS) or_return
+ _private_div_recursion(Q1, R1, A1, B1) or_return
/*
A1 = (R1 * beta^(2k)) + (A % beta^(2k)) - (Q1 * B0 * beta^k)
*/
- internal_shl_digit(R1, 2 * k) or_return;
- internal_add(A1, R1, t) or_return;
- internal_mul(t, Q1, B0) or_return;
+ internal_shl_digit(R1, 2 * k) or_return
+ internal_add(A1, R1, t) or_return
+ internal_mul(t, Q1, B0) or_return
/*
While A1 < 0 do Q1 = Q1 - 1, A1 = A1 + (beta^k * B)
*/
if internal_cmp(A1, 0) == -1 {
- internal_shl(t, b, k * _DIGIT_BITS) or_return;
+ internal_shl(t, b, k * _DIGIT_BITS) or_return
for {
- internal_decr(Q1) or_return;
- internal_add(A1, A1, t) or_return;
+ internal_decr(Q1) or_return
+ internal_add(A1, A1, t) or_return
if internal_cmp(A1, 0) != -1 {
- break;
+ break
}
}
}
@@ -1243,72 +1243,72 @@ _private_div_recursion :: proc(quotient, remainder, a, b: ^Int, allocator := con /*
(Q0, R0) = RecursiveDivRem(A1 / beta^(k), B1)
*/
- internal_shrmod(A1, t, A1, k * _DIGIT_BITS) or_return;
- _private_div_recursion(Q0, R0, A1, B1) or_return;
+ internal_shrmod(A1, t, A1, k * _DIGIT_BITS) or_return
+ _private_div_recursion(Q0, R0, A1, B1) or_return
/*
A2 = (R0*beta^k) + (A1 % beta^k) - (Q0*B0)
*/
- internal_shl_digit(R0, k) or_return;
- internal_add(A2, R0, t) or_return;
- internal_mul(t, Q0, B0) or_return;
- internal_sub(A2, A2, t) or_return;
+ internal_shl_digit(R0, k) or_return
+ internal_add(A2, R0, t) or_return
+ internal_mul(t, Q0, B0) or_return
+ internal_sub(A2, A2, t) or_return
/*
While A2 < 0 do Q0 = Q0 - 1, A2 = A2 + B.
*/
for internal_cmp(A2, 0) == -1 {
- internal_decr(Q0) or_return;
- internal_add(A2, A2, b) or_return;
+ internal_decr(Q0) or_return
+ internal_add(A2, A2, b) or_return
}
/*
Return q = (Q1*beta^k) + Q0, r = A2.
*/
- internal_shl_digit(Q1, k) or_return;
- internal_add(quotient, Q1, Q0) or_return;
+ internal_shl_digit(Q1, k) or_return
+ internal_add(quotient, Q1, Q0) or_return
- return internal_copy(remainder, A2);
+ return internal_copy(remainder, A2)
}
_private_int_div_recursive :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- A, B, Q, Q1, R, A_div, A_mod := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(A, B, Q, Q1, R, A_div, A_mod);
+ A, B, Q, Q1, R, A_div, A_mod := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(A, B, Q, Q1, R, A_div, A_mod)
- internal_init_multi(A, B, Q, Q1, R, A_div, A_mod) or_return;
+ internal_init_multi(A, B, Q, Q1, R, A_div, A_mod) or_return
/*
Most significant bit of a limb.
Assumes _DIGIT_MAX < (sizeof(DIGIT) * sizeof(u8)).
*/
- msb := (_DIGIT_MAX + DIGIT(1)) >> 1;
- sigma := 0;
- msb_b := b.digit[b.used - 1];
+ msb := (_DIGIT_MAX + DIGIT(1)) >> 1
+ sigma := 0
+ msb_b := b.digit[b.used - 1]
for msb_b < msb {
- sigma += 1;
- msb_b <<= 1;
+ sigma += 1
+ msb_b <<= 1
}
/*
Use that sigma to normalize B.
*/
- internal_shl(B, b, sigma) or_return;
- internal_shl(A, a, sigma) or_return;
+ internal_shl(B, b, sigma) or_return
+ internal_shl(A, a, sigma) or_return
/*
Fix the sign.
*/
- neg := a.sign != b.sign;
- A.sign = .Zero_or_Positive; B.sign = .Zero_or_Positive;
+ neg := a.sign != b.sign
+ A.sign = .Zero_or_Positive; B.sign = .Zero_or_Positive
/*
If the magnitude of "A" is not more more than twice that of "B" we can work
on them directly, otherwise we need to work at "A" in chunks.
*/
- n := B.used;
- m := A.used - B.used;
+ n := B.used
+ m := A.used - B.used
/*
Q = 0. We already ensured that when we called `internal_init_multi`.
@@ -1317,56 +1317,56 @@ _private_int_div_recursive :: proc(quotient, remainder, a, b: ^Int, allocator := /*
(q, r) = RecursiveDivRem(A / (beta^(m-n)), B)
*/
- j := (m - n) * _DIGIT_BITS;
- internal_shrmod(A_div, A_mod, A, j) or_return;
- _private_div_recursion(Q1, R, A_div, B) or_return;
+ j := (m - n) * _DIGIT_BITS
+ internal_shrmod(A_div, A_mod, A, j) or_return
+ _private_div_recursion(Q1, R, A_div, B) or_return
/*
Q = (Q*beta!(n)) + q
*/
- internal_shl(Q, Q, n * _DIGIT_BITS) or_return;
- internal_add(Q, Q, Q1) or_return;
+ internal_shl(Q, Q, n * _DIGIT_BITS) or_return
+ internal_add(Q, Q, Q1) or_return
/*
A = (r * beta^(m-n)) + (A % beta^(m-n))
*/
- internal_shl(R, R, (m - n) * _DIGIT_BITS) or_return;
- internal_add(A, R, A_mod) or_return;
+ internal_shl(R, R, (m - n) * _DIGIT_BITS) or_return
+ internal_add(A, R, A_mod) or_return
/*
m = m - n
*/
- m -= n;
+ m -= n
}
/*
(q, r) = RecursiveDivRem(A, B)
*/
- _private_div_recursion(Q1, R, A, B) or_return;
+ _private_div_recursion(Q1, R, A, B) or_return
/*
Q = (Q * beta^m) + q, R = r
*/
- internal_shl(Q, Q, m * _DIGIT_BITS) or_return;
- internal_add(Q, Q, Q1) or_return;
+ internal_shl(Q, Q, m * _DIGIT_BITS) or_return
+ internal_add(Q, Q, Q1) or_return
/*
Get sign before writing to dest.
*/
- R.sign = .Zero_or_Positive if internal_is_zero(Q) else a.sign;
+ R.sign = .Zero_or_Positive if internal_is_zero(Q) else a.sign
if quotient != nil {
- swap(quotient, Q);
- quotient.sign = .Negative if neg else .Zero_or_Positive;
+ swap(quotient, Q)
+ quotient.sign = .Negative if neg else .Zero_or_Positive
}
if remainder != nil {
/*
De-normalize the remainder.
*/
- internal_shrmod(R, nil, R, sigma) or_return;
- swap(remainder, R);
+ internal_shrmod(R, nil, R, sigma) or_return
+ swap(remainder, R)
}
- return nil;
+ return nil
}
/*
@@ -1375,53 +1375,53 @@ _private_int_div_recursive :: proc(quotient, remainder, a, b: ^Int, allocator := @(deprecated="Use `_int_div_school`, it's 3.5x faster.")
_private_int_div_small :: proc(quotient, remainder, numerator, denominator: ^Int) -> (err: Error) {
- ta, tb, tq, q := &Int{}, &Int{}, &Int{}, &Int{};
- c: int;
- defer internal_destroy(ta, tb, tq, q);
+ ta, tb, tq, q := &Int{}, &Int{}, &Int{}, &Int{}
+ c: int
+ defer internal_destroy(ta, tb, tq, q)
for {
- internal_one(tq) or_return;
+ internal_one(tq) or_return
- num_bits, _ := count_bits(numerator);
- den_bits, _ := count_bits(denominator);
- n := num_bits - den_bits;
+ num_bits, _ := count_bits(numerator)
+ den_bits, _ := count_bits(denominator)
+ n := num_bits - den_bits
- abs(ta, numerator) or_return;
- abs(tb, denominator) or_return;
- shl(tb, tb, n) or_return;
- shl(tq, tq, n) or_return;
+ abs(ta, numerator) or_return
+ abs(tb, denominator) or_return
+ shl(tb, tb, n) or_return
+ shl(tq, tq, n) or_return
for n >= 0 {
if c, _ = cmp_mag(ta, tb); c == 0 || c == 1 {
// ta -= tb
- sub(ta, ta, tb) or_return;
+ sub(ta, ta, tb) or_return
// q += tq
- add( q, q, tq) or_return;
+ add( q, q, tq) or_return
}
- shr1(tb, tb) or_return;
- shr1(tq, tq) or_return;
+ shr1(tb, tb) or_return
+ shr1(tq, tq) or_return
- n -= 1;
+ n -= 1
}
/*
Now q == quotient and ta == remainder.
*/
- neg := numerator.sign != denominator.sign;
+ neg := numerator.sign != denominator.sign
if quotient != nil {
- swap(quotient, q);
- z, _ := is_zero(quotient);
- quotient.sign = .Negative if neg && !z else .Zero_or_Positive;
+ swap(quotient, q)
+ z, _ := is_zero(quotient)
+ quotient.sign = .Negative if neg && !z else .Zero_or_Positive
}
if remainder != nil {
- swap(remainder, ta);
- z, _ := is_zero(numerator);
- remainder.sign = .Zero_or_Positive if z else numerator.sign;
+ swap(remainder, ta)
+ z, _ := is_zero(numerator)
+ remainder.sign = .Zero_or_Positive if z else numerator.sign
}
- break;
+ break
}
- return err;
+ return err
}
@@ -1430,65 +1430,65 @@ _private_int_div_small :: proc(quotient, remainder, numerator, denominator: ^Int Binary split factorial algo due to: http://www.luschny.de/math/factorial/binarysplitfact.html
*/
_private_int_factorial_binary_split :: proc(res: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- inner, outer, start, stop, temp := &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(inner, outer, start, stop, temp);
+ inner, outer, start, stop, temp := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(inner, outer, start, stop, temp)
- internal_one(inner, false) or_return;
- internal_one(outer, false) or_return;
+ internal_one(inner, false) or_return
+ internal_one(outer, false) or_return
- bits_used := int(_DIGIT_TYPE_BITS - intrinsics.count_leading_zeros(n));
+ bits_used := int(_DIGIT_TYPE_BITS - intrinsics.count_leading_zeros(n))
for i := bits_used; i >= 0; i -= 1 {
- start := (n >> (uint(i) + 1)) + 1 | 1;
- stop := (n >> uint(i)) + 1 | 1;
- _private_int_recursive_product(temp, start, stop, 0) or_return;
- internal_mul(inner, inner, temp) or_return;
- internal_mul(outer, outer, inner) or_return;
+ start := (n >> (uint(i) + 1)) + 1 | 1
+ stop := (n >> uint(i)) + 1 | 1
+ _private_int_recursive_product(temp, start, stop, 0) or_return
+ internal_mul(inner, inner, temp) or_return
+ internal_mul(outer, outer, inner) or_return
}
- shift := n - intrinsics.count_ones(n);
+ shift := n - intrinsics.count_ones(n)
- return internal_shl(res, outer, int(shift));
+ return internal_shl(res, outer, int(shift))
}
/*
Recursive product used by binary split factorial algorithm.
*/
_private_int_recursive_product :: proc(res: ^Int, start, stop: int, level := int(0), allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
- t1, t2 := &Int{}, &Int{};
- defer internal_destroy(t1, t2);
+ t1, t2 := &Int{}, &Int{}
+ defer internal_destroy(t1, t2)
if level > FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS {
- return .Max_Iterations_Reached;
+ return .Max_Iterations_Reached
}
- num_factors := (stop - start) >> 1;
+ num_factors := (stop - start) >> 1
if num_factors == 2 {
- internal_set(t1, start, false) or_return;
+ internal_set(t1, start, false) or_return
when true {
- internal_grow(t2, t1.used + 1, false) or_return;
- internal_add(t2, t1, 2) or_return;
+ internal_grow(t2, t1.used + 1, false) or_return
+ internal_add(t2, t1, 2) or_return
} else {
- internal_add(t2, t1, 2) or_return;
+ internal_add(t2, t1, 2) or_return
}
- return internal_mul(res, t1, t2);
+ return internal_mul(res, t1, t2)
}
if num_factors > 1 {
- mid := (start + num_factors) | 1;
- _private_int_recursive_product(t1, start, mid, level + 1) or_return;
- _private_int_recursive_product(t2, mid, stop, level + 1) or_return;
- return internal_mul(res, t1, t2);
+ mid := (start + num_factors) | 1
+ _private_int_recursive_product(t1, start, mid, level + 1) or_return
+ _private_int_recursive_product(t2, mid, stop, level + 1) or_return
+ return internal_mul(res, t1, t2)
}
if num_factors == 1 {
- return #force_inline internal_set(res, start, true);
+ return #force_inline internal_set(res, start, true)
}
- return #force_inline internal_one(res, true);
+ return #force_inline internal_one(res, true)
}
/*
@@ -1507,10 +1507,10 @@ _private_int_recursive_product :: proc(res: ^Int, start, stop: int, level := int If neither result is wanted, we have nothing to do.
*/
_private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
+ context.allocator = allocator
if res_gcd == nil && res_lcm == nil {
- return nil;
+ return nil
}
/*
@@ -1521,76 +1521,76 @@ _private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context. GCD(0, 0) and LCM(0, 0) are both 0.
*/
if res_gcd != nil {
- internal_zero(res_gcd) or_return;
+ internal_zero(res_gcd) or_return
}
if res_lcm != nil {
- internal_zero(res_lcm) or_return;
+ internal_zero(res_lcm) or_return
}
- return nil;
+ return nil
} else if a.used == 0 {
/*
We can early out with GCD = B and LCM = 0
*/
if res_gcd != nil {
- internal_abs(res_gcd, b) or_return;
+ internal_abs(res_gcd, b) or_return
}
if res_lcm != nil {
- internal_zero(res_lcm) or_return;
+ internal_zero(res_lcm) or_return
}
- return nil;
+ return nil
} else if b.used == 0 {
/*
We can early out with GCD = A and LCM = 0
*/
if res_gcd != nil {
- internal_abs(res_gcd, a) or_return;
+ internal_abs(res_gcd, a) or_return
}
if res_lcm != nil {
- internal_zero(res_lcm) or_return;
+ internal_zero(res_lcm) or_return
}
- return nil;
+ return nil
}
- temp_gcd_res := &Int{};
- defer internal_destroy(temp_gcd_res);
+ temp_gcd_res := &Int{}
+ defer internal_destroy(temp_gcd_res)
/*
If neither `a` or `b` was zero, we need to compute `gcd`.
Get copies of `a` and `b` we can modify.
*/
- u, v := &Int{}, &Int{};
- defer internal_destroy(u, v);
- internal_copy(u, a) or_return;
- internal_copy(v, b) or_return;
+ u, v := &Int{}, &Int{}
+ defer internal_destroy(u, v)
+ internal_copy(u, a) or_return
+ internal_copy(v, b) or_return
/*
Must be positive for the remainder of the algorithm.
*/
- u.sign = .Zero_or_Positive; v.sign = .Zero_or_Positive;
+ u.sign = .Zero_or_Positive; v.sign = .Zero_or_Positive
/*
B1. Find the common power of two for `u` and `v`.
*/
- u_lsb, _ := internal_count_lsb(u);
- v_lsb, _ := internal_count_lsb(v);
- k := min(u_lsb, v_lsb);
+ u_lsb, _ := internal_count_lsb(u)
+ v_lsb, _ := internal_count_lsb(v)
+ k := min(u_lsb, v_lsb)
if k > 0 {
/*
Divide the power of two out.
*/
- internal_shr(u, u, k) or_return;
- internal_shr(v, v, k) or_return;
+ internal_shr(u, u, k) or_return
+ internal_shr(v, v, k) or_return
}
/*
Divide any remaining factors of two out.
*/
if u_lsb != k {
- internal_shr(u, u, u_lsb - k) or_return;
+ internal_shr(u, u, u_lsb - k) or_return
}
if v_lsb != k {
- internal_shr(v, v, v_lsb - k) or_return;
+ internal_shr(v, v, v_lsb - k) or_return
}
for v.used != 0 {
@@ -1601,34 +1601,34 @@ _private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context. /*
Swap `u` and `v` to make sure `v` is >= `u`.
*/
- internal_swap(u, v);
+ internal_swap(u, v)
}
/*
Subtract smallest from largest.
*/
- internal_sub(v, v, u) or_return;
+ internal_sub(v, v, u) or_return
/*
Divide out all factors of two.
*/
- b, _ := internal_count_lsb(v);
- internal_shr(v, v, b) or_return;
+ b, _ := internal_count_lsb(v)
+ internal_shr(v, v, b) or_return
}
/*
Multiply by 2**k which we divided out at the beginning.
*/
- internal_shl(temp_gcd_res, u, k) or_return;
- temp_gcd_res.sign = .Zero_or_Positive;
+ internal_shl(temp_gcd_res, u, k) or_return
+ temp_gcd_res.sign = .Zero_or_Positive
/*
We've computed `gcd`, either the long way, or because one of the inputs was zero.
If we don't want `lcm`, we're done.
*/
if res_lcm == nil {
- internal_swap(temp_gcd_res, res_gcd);
- return nil;
+ internal_swap(temp_gcd_res, res_gcd)
+ return nil
}
/*
@@ -1639,25 +1639,25 @@ _private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context. /*
Store quotient in `t2` such that `t2 * b` is the LCM.
*/
- internal_div(res_lcm, a, temp_gcd_res) or_return;
- err = internal_mul(res_lcm, res_lcm, b);
+ internal_div(res_lcm, a, temp_gcd_res) or_return
+ err = internal_mul(res_lcm, res_lcm, b)
} else {
/*
Store quotient in `t2` such that `t2 * a` is the LCM.
*/
- internal_div(res_lcm, a, temp_gcd_res) or_return;
- err = internal_mul(res_lcm, res_lcm, b);
+ internal_div(res_lcm, a, temp_gcd_res) or_return
+ err = internal_mul(res_lcm, res_lcm, b)
}
if res_gcd != nil {
- internal_swap(temp_gcd_res, res_gcd);
+ internal_swap(temp_gcd_res, res_gcd)
}
/*
Fix the sign to positive and return.
*/
- res_lcm.sign = .Zero_or_Positive;
- return err;
+ res_lcm.sign = .Zero_or_Positive
+ return err
}
/*
@@ -1665,24 +1665,24 @@ _private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context. Assumes `a` not to be `nil` and to have been initialized.
*/
_private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) -> (res: int, err: Error) {
- bracket_low, bracket_high, bracket_mid, t, bi_base := &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
- defer internal_destroy(bracket_low, bracket_high, bracket_mid, t, bi_base);
+ bracket_low, bracket_high, bracket_mid, t, bi_base := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(bracket_low, bracket_high, bracket_mid, t, bi_base)
- ic := #force_inline internal_cmp(a, base);
+ ic := #force_inline internal_cmp(a, base)
if ic == -1 || ic == 0 {
- return 1 if ic == 0 else 0, nil;
+ return 1 if ic == 0 else 0, nil
}
defer if err != nil {
- res = -1;
+ res = -1
}
- internal_set(bi_base, base, true, allocator) or_return;
- internal_clear(bracket_mid, false, allocator) or_return;
- internal_clear(t, false, allocator) or_return;
- internal_one(bracket_low, false, allocator) or_return;
- internal_set(bracket_high, base, false, allocator) or_return;
+ internal_set(bi_base, base, true, allocator) or_return
+ internal_clear(bracket_mid, false, allocator) or_return
+ internal_clear(t, false, allocator) or_return
+ internal_one(bracket_low, false, allocator) or_return
+ internal_set(bracket_high, base, false, allocator) or_return
- low := 0; high := 1;
+ low := 0; high := 1
/*
A kind of Giant-step/baby-step algorithm.
@@ -1696,39 +1696,39 @@ _private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) - Iterate until `a` is bracketed between low + high.
*/
if #force_inline internal_cmp(bracket_high, a) != -1 {
- break;
+ break
}
- low = high;
+ low = high
#force_inline internal_copy(bracket_low, bracket_high) or_return;
- high <<= 1;
+ high <<= 1
#force_inline internal_sqr(bracket_high, bracket_high) or_return;
}
for (high - low) > 1 {
- mid := (high + low) >> 1;
+ mid := (high + low) >> 1
#force_inline internal_pow(t, bi_base, mid - low) or_return;
#force_inline internal_mul(bracket_mid, bracket_low, t) or_return;
- mc := #force_inline internal_cmp(a, bracket_mid);
+ mc := #force_inline internal_cmp(a, bracket_mid)
switch mc {
case -1:
- high = mid;
- internal_swap(bracket_mid, bracket_high);
+ high = mid
+ internal_swap(bracket_mid, bracket_high)
case 0:
- return mid, nil;
+ return mid, nil
case 1:
- low = mid;
- internal_swap(bracket_mid, bracket_low);
+ low = mid
+ internal_swap(bracket_mid, bracket_low)
}
}
- fc := #force_inline internal_cmp(bracket_high, a);
- res = high if fc == 0 else low;
+ fc := #force_inline internal_cmp(bracket_high, a)
+ res = high if fc == 0 else low
- return;
+ return
}
@@ -1741,37 +1741,37 @@ _private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) - Based on Algorithm 14.32 on pp.601 of HAC.
*/
_private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {
- context.allocator = allocator;
- W: [_WARRAY]_WORD = ---;
+ context.allocator = allocator
+ W: [_WARRAY]_WORD = ---
if x.used > _WARRAY { return .Invalid_Argument; }
/*
Get old used count.
*/
- old_used := x.used;
+ old_used := x.used
/*
Grow `x` as required.
*/
- internal_grow(x, n.used + 1) or_return;
+ internal_grow(x, n.used + 1) or_return
/*
First we have to get the digits of the input into an array of double precision words W[...]
Copy the digits of `x` into W[0..`x.used` - 1]
*/
- ix: int;
+ ix: int
for ix = 0; ix < x.used; ix += 1 {
- W[ix] = _WORD(x.digit[ix]);
+ W[ix] = _WORD(x.digit[ix])
}
/*
Zero the high words of W[a->used..m->used*2].
*/
- zero_upper := (n.used * 2) + 1;
+ zero_upper := (n.used * 2) + 1
if ix < zero_upper {
for ix = x.used; ix < zero_upper; ix += 1 {
- W[ix] = {};
+ W[ix] = {}
}
}
@@ -1786,7 +1786,7 @@ _private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := co by casting the value down to a DIGIT. Note this requires
that W[ix-1] have the carry cleared (see after the inner loop)
*/
- mu := ((W[ix] & _WORD(_MASK)) * _WORD(rho)) & _WORD(_MASK);
+ mu := ((W[ix] & _WORD(_MASK)) * _WORD(rho)) & _WORD(_MASK)
/*
`a = a + mu * m * b**i`
@@ -1805,13 +1805,13 @@ _private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := co first m->used words of W[] have the carries fixed.
*/
for iy := 0; iy < n.used; iy += 1 {
- W[ix + iy] += mu * _WORD(n.digit[iy]);
+ W[ix + iy] += mu * _WORD(n.digit[iy])
}
/*
Now fix carry for next digit, W[ix+1].
*/
- W[ix + 1] += (W[ix] >> _DIGIT_BITS);
+ W[ix + 1] += (W[ix] >> _DIGIT_BITS)
}
/*
@@ -1820,7 +1820,7 @@ _private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := co */
for ; ix < n.used * 2; ix += 1 {
- W[ix + 1] += (W[ix] >> _DIGIT_BITS);
+ W[ix + 1] += (W[ix] >> _DIGIT_BITS)
}
/* copy out, A = A/b**n
@@ -1831,69 +1831,69 @@ _private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := co */
for ix = 0; ix < (n.used + 1); ix += 1 {
- x.digit[ix] = DIGIT(W[n.used + ix] & _WORD(_MASK));
+ x.digit[ix] = DIGIT(W[n.used + ix] & _WORD(_MASK))
}
/*
Set the max used.
*/
- x.used = n.used + 1;
+ x.used = n.used + 1
/*
Zero old_used digits, if the input a was larger than m->used+1 we'll have to clear the digits.
*/
- internal_zero_unused(x, old_used);
- internal_clamp(x);
+ internal_zero_unused(x, old_used)
+ internal_clamp(x)
/*
if A >= m then A = A - m
*/
if internal_cmp_mag(x, n) != -1 {
- return internal_sub(x, x, n);
+ return internal_sub(x, x, n)
}
- return nil;
+ return nil
}
/*
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 internal_destroy(x, y, u, v, A, B, C, D);
+ context.allocator = allocator
+ x, y, u, v, A, B, C, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(x, y, u, v, A, B, C, D)
/*
`b` cannot be negative.
*/
if b.sign == .Negative || internal_is_zero(b) {
- return .Invalid_Argument;
+ return .Invalid_Argument
}
/*
init temps.
*/
- internal_init_multi(x, y, u, v, A, B, C, D) or_return;
+ internal_init_multi(x, y, u, v, A, B, C, D) or_return
/*
`x` = `a` % `b`, `y` = `b`
*/
- internal_mod(x, a, b) or_return;
- internal_copy(y, b) or_return;
+ internal_mod(x, a, b) or_return
+ internal_copy(y, b) or_return
/*
2. [modified] if x,y are both even then return an error!
*/
if internal_is_even(x) && internal_is_even(y) {
- return .Invalid_Argument;
+ return .Invalid_Argument
}
/*
3. u=x, v=y, A=1, B=0, C=0, D=1
*/
- internal_copy(u, x) or_return;
- internal_copy(v, y) or_return;
- internal_one(A) or_return;
- internal_one(D) or_return;
+ internal_copy(u, x) or_return
+ internal_copy(v, y) or_return
+ internal_one(A) or_return
+ internal_one(D) or_return
for {
/*
@@ -1903,7 +1903,7 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator /*
4.1 `u` = `u` / 2
*/
- internal_int_shr1(u, u) or_return;
+ internal_int_shr1(u, u) or_return
/*
4.2 if `A` or `B` is odd then:
@@ -1912,14 +1912,14 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator /*
`A` = (`A`+`y`) / 2, `B` = (`B`-`x`) / 2
*/
- internal_add(A, A, y) or_return;
- internal_add(B, B, x) or_return;
+ internal_add(A, A, y) or_return
+ internal_add(B, B, x) or_return
}
/*
`A` = `A` / 2, `B` = `B` / 2
*/
- internal_int_shr1(A, A) or_return;
- internal_int_shr1(B, B) or_return;
+ internal_int_shr1(A, A) or_return
+ internal_int_shr1(B, B) or_return
}
/*
@@ -1929,7 +1929,7 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator /*
5.1 `v` = `v` / 2
*/
- internal_int_shr1(v, v) or_return;
+ internal_int_shr1(v, v) or_return
/*
5.2 if `C` or `D` is odd then:
@@ -1938,14 +1938,14 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator /*
`C` = (`C`+`y`) / 2, `D` = (`D`-`x`) / 2
*/
- internal_add(C, C, y) or_return;
- internal_add(D, D, x) or_return;
+ internal_add(C, C, y) or_return
+ internal_add(D, D, x) or_return
}
/*
`C` = `C` / 2, `D` = `D` / 2
*/
- internal_int_shr1(C, C) or_return;
- internal_int_shr1(D, D) or_return;
+ internal_int_shr1(C, C) or_return
+ internal_int_shr1(D, D) or_return
}
/*
@@ -1955,21 +1955,21 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator /*
`u` = `u` - `v`, `A` = `A` - `C`, `B` = `B` - `D`
*/
- internal_sub(u, u, v) or_return;
- internal_sub(A, A, C) or_return;
- internal_sub(B, B, D) or_return;
+ internal_sub(u, u, v) or_return
+ internal_sub(A, A, C) or_return
+ internal_sub(B, B, D) or_return
} else {
/* v - v - u, C = C - A, D = D - B */
- internal_sub(v, v, u) or_return;
- internal_sub(C, C, A) or_return;
- internal_sub(D, D, B) or_return;
+ internal_sub(v, v, u) or_return
+ internal_sub(C, C, A) or_return
+ internal_sub(D, D, B) or_return
}
/*
If not zero goto step 4
*/
if internal_is_zero(u) {
- break;
+ break
}
}
@@ -1981,29 +1981,29 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator If `v` != `1` then there is no inverse.
*/
if internal_cmp(v, 1) != 0 {
- return .Invalid_Argument;
+ return .Invalid_Argument
}
/*
If its too low.
*/
if internal_cmp(C, 0) == -1 {
- internal_add(C, C, b) or_return;
+ internal_add(C, C, b) or_return
}
/*
Too big.
*/
if internal_cmp(C, 0) != -1 {
- internal_sub(C, C, b) or_return;
+ internal_sub(C, C, b) or_return
}
/*
`C` is now the inverse.
*/
- swap(dest, C);
+ swap(dest, C)
- return;
+ return
}
/*
@@ -2013,11 +2013,11 @@ _private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator 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 internal_destroy(x, y, u, v, B, D);
+ context.allocator = allocator
+ x, y, u, v, B, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+ defer internal_destroy(x, y, u, v, B, D)
- sign: Sign;
+ sign: Sign
/*
2. [modified] `b` must be odd.
@@ -2027,17 +2027,17 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
Init all our temps.
*/
- internal_init_multi(x, y, u, v, B, D) or_return;
+ internal_init_multi(x, y, u, v, B, D) or_return
/*
`x` == modulus, `y` == value to invert.
*/
- internal_copy(x, b) or_return;
+ internal_copy(x, b) or_return
/*
We need `y` = `|a|`.
*/
- internal_mod(y, a, b) or_return;
+ internal_mod(y, a, b) or_return
/*
If one of `x`, `y` is zero return an error!
@@ -2047,10 +2047,10 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
3. `u` = `x`, `v` = `y`, `A` = 1, `B` = 0, `C` = 0, `D` = 1
*/
- internal_copy(u, x) or_return;
- internal_copy(v, y) or_return;
+ internal_copy(u, x) or_return
+ internal_copy(v, y) or_return
- internal_one(D) or_return;
+ internal_one(D) or_return
for {
/*
@@ -2060,7 +2060,7 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
4.1 `u` = `u` / 2
*/
- internal_int_shr1(u, u) or_return;
+ internal_int_shr1(u, u) or_return
/*
4.2 if `B` is odd then:
@@ -2069,13 +2069,13 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
`B` = (`B` - `x`) / 2
*/
- internal_sub(B, B, x) or_return;
+ internal_sub(B, B, x) or_return
}
/*
`B` = `B` / 2
*/
- internal_int_shr1(B, B) or_return;
+ internal_int_shr1(B, B) or_return
}
/*
@@ -2085,7 +2085,7 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
5.1 `v` = `v` / 2
*/
- internal_int_shr1(v, v) or_return;
+ internal_int_shr1(v, v) or_return
/*
5.2 if `D` is odd then:
@@ -2094,12 +2094,12 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
`D` = (`D` - `x`) / 2
*/
- internal_sub(D, D, x) or_return;
+ internal_sub(D, D, x) or_return
}
/*
`D` = `D` / 2
*/
- internal_int_shr1(D, D) or_return;
+ internal_int_shr1(D, D) or_return
}
/*
@@ -2109,14 +2109,14 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc /*
`u` = `u` - `v`, `B` = `B` - `D`
*/
- internal_sub(u, u, v) or_return;
- internal_sub(B, B, D) or_return;
+ internal_sub(u, u, v) or_return
+ internal_sub(B, B, D) or_return
} else {
/*
`v` - `v` - `u`, `D` = `D` - `B`
*/
- internal_sub(v, v, u) or_return;
- internal_sub(D, D, B) or_return;
+ internal_sub(v, v, u) or_return
+ internal_sub(D, D, B) or_return
}
/*
@@ -2133,27 +2133,27 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc if `v` != 1 then there is no inverse
*/
if internal_cmp(v, 1) != 0 {
- return .Invalid_Argument;
+ return .Invalid_Argument
}
/*
`b` is now the inverse.
*/
- sign = a.sign;
+ sign = a.sign
for internal_int_is_negative(D) {
- internal_add(D, D, b) or_return;
+ internal_add(D, D, b) or_return
}
/*
Too big.
*/
for internal_cmp_mag(D, b) != -1 {
- internal_sub(D, D, b) or_return;
+ internal_sub(D, D, b) or_return
}
- swap(dest, D);
- dest.sign = sign;
- return nil;
+ swap(dest, D)
+ dest.sign = sign
+ return nil
}
@@ -2163,14 +2163,14 @@ _private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.alloc Also assumes `base` is a power of two.
*/
_private_log_power_of_two :: proc(a: ^Int, base: DIGIT) -> (log: int, err: Error) {
- base := base;
- y: int;
+ base := base
+ y: int
for y = 0; base & 1 == 0; {
- y += 1;
- base >>= 1;
+ y += 1
+ base >>= 1
}
- log = internal_count_bits(a);
- return (log - 1) / y, err;
+ log = internal_count_bits(a)
+ return (log - 1) / y, err
}
/*
@@ -2178,17 +2178,17 @@ _private_log_power_of_two :: proc(a: ^Int, base: DIGIT) -> (log: int, err: Error Assumes `src` and `dest` to not be `nil` and have been initialized.
*/
_private_copy_digits :: proc(dest, src: ^Int, digits: int, offset := int(0)) -> (err: Error) {
- digits := digits;
+ digits := digits
/*
If dest == src, do nothing
*/
if dest == src {
- return nil;
+ return nil
}
- digits = min(digits, len(src.digit), len(dest.digit));
- mem.copy_non_overlapping(&dest.digit[0], &src.digit[offset], size_of(DIGIT) * digits);
- return nil;
+ digits = min(digits, len(src.digit), len(dest.digit))
+ mem.copy_non_overlapping(&dest.digit[0], &src.digit[offset], size_of(DIGIT) * digits)
+ return nil
}
/*
@@ -2208,7 +2208,7 @@ _private_int_rem_128 := [?]DIGIT{ 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
-};
+}
#assert(128 * size_of(DIGIT) == size_of(_private_int_rem_128));
_private_int_rem_105 := [?]DIGIT{
@@ -2219,7 +2219,7 @@ _private_int_rem_105 := [?]DIGIT{ 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,
-};
+}
#assert(105 * size_of(DIGIT) == size_of(_private_int_rem_105));
_private_prime_table := [?]DIGIT{
@@ -2258,7 +2258,7 @@ _private_prime_table := [?]DIGIT{ 0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653,
-};
+}
#assert(256 * size_of(DIGIT) == size_of(_private_prime_table));
when MATH_BIG_FORCE_64_BIT || (!MATH_BIG_FORCE_32_BIT && size_of(rawptr) == 8) {
@@ -2298,7 +2298,7 @@ when MATH_BIG_FORCE_64_BIT || (!MATH_BIG_FORCE_32_BIT && size_of(rawptr) == 8) { /* f(32): */ 263_130_836_933_693_530_167_218_012_160_000_000,
/* f(33): */ 8_683_317_618_811_886_495_518_194_401_280_000_000,
/* f(34): */ 295_232_799_039_604_140_847_618_609_643_520_000_000,
- };
+ }
} else {
_factorial_table := [21]_WORD{
/* f(00): */ 1,
@@ -2322,7 +2322,7 @@ when MATH_BIG_FORCE_64_BIT || (!MATH_BIG_FORCE_32_BIT && size_of(rawptr) == 8) { /* f(18): */ 6_402_373_705_728_000,
/* f(19): */ 121_645_100_408_832_000,
/* f(20): */ 2_432_902_008_176_640_000,
- };
+ }
};
/*
|