diff options
| -rw-r--r-- | core/math/big/example.odin | 19 | ||||
| -rw-r--r-- | core/math/big/internal.odin | 18 | ||||
| -rw-r--r-- | core/math/big/private.odin | 66 |
3 files changed, 91 insertions, 12 deletions
diff --git a/core/math/big/example.odin b/core/math/big/example.odin index 3d51ac400..1d9a4cb3c 100644 --- a/core/math/big/example.odin +++ b/core/math/big/example.odin @@ -80,16 +80,27 @@ demo :: proc() { err: Error; bs: string; - if err = factorial(a, 500); err != nil { fmt.printf("factorial err: %v\n", err); return; } + // if err = factorial(a, 850); err != nil { fmt.printf("factorial err: %v\n", err); return; } + + foo := "615037959146039477924633848896619112832171971562900618409305032006863881436080"; + if err = atoi(a, foo, 10); err != nil { return; } + print("a: ", a, 10, true, true, true); + fmt.println(); + { SCOPED_TIMING(.sqr); - if err = sqr(b, a); err != nil { fmt.printf("sqr err: %v\n", err); return; } + if err = sqr(b, a); err != nil { fmt.printf("sqr err: %v\n", err); return; } } + fmt.println(); + print("b _sqr_karatsuba: ", b); + fmt.println(); - bs, err = itoa(b, 10); + bs, err = itoa(b, 16); defer delete(bs); - assert(bs[:50] == "14887338741396604108836218987068397819515734169330"); + if bs[:50] != "1C367982F3050A8A3C62A8A7906D165438B54B287AF3F15D36" { + fmt.println("sqr failed"); + } } main :: proc() { diff --git a/core/math/big/internal.odin b/core/math/big/internal.odin index 9b267176d..00b1ed7bf 100644 --- a/core/math/big/internal.odin +++ b/core/math/big/internal.odin @@ -36,7 +36,7 @@ import "core:mem" import "core:intrinsics" import rnd "core:math/rand" -//import "core:fmt" +// import "core:fmt" /* Low-level addition, unsigned. Handbook of Applied Cryptography, algorithm 14.7. @@ -627,20 +627,22 @@ internal_int_mul :: proc(dest, src, multiplier: ^Int, allocator := context.alloc Do we need to square? */ if src.used >= SQR_TOOM_CUTOFF { - /* Use Toom-Cook? */ - // err = s_mp_sqr_toom(a, c); + /* + Use Toom-Cook? + */ // fmt.printf("_private_int_sqr_toom: %v\n", src.used); - err = #force_inline _private_int_sqr(dest, src); + err = #force_inline _private_int_sqr_karatsuba(dest, src); } else if src.used >= SQR_KARATSUBA_CUTOFF { - /* Karatsuba? */ - // err = s_mp_sqr_karatsuba(a, c); - // fmt.printf("_private_int_sqr_karatsuba: %v\n", src.used); - err = #force_inline _private_int_sqr(dest, src); + /* + Karatsuba? + */ + err = #force_inline _private_int_sqr_karatsuba(dest, src); } else if ((src.used * 2) + 1) < _WARRAY && src.used < (_MAX_COMBA / 2) { /* Fast comba? */ err = #force_inline _private_int_sqr_comba(dest, src); + //err = #force_inline _private_int_sqr(dest, src); } else { err = #force_inline _private_int_sqr(dest, src); } diff --git a/core/math/big/private.odin b/core/math/big/private.odin index f49cbc51a..3d8497c72 100644 --- a/core/math/big/private.odin +++ b/core/math/big/private.odin @@ -355,6 +355,72 @@ _private_int_sqr_comba :: proc(dest, src: ^Int, allocator := context.allocator) }
/*
+ 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.
+*/
+_private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
+ context.allocator = allocator;
+
+ 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;
+
+ /*
+ Init temps.
+ */
+ if err = internal_grow(x0, B); err != nil { return err; }
+ if err = internal_grow(x1, src.used - B); err != nil { return err; }
+ if err = internal_grow(t1, src.used * 2); err != nil { return err; }
+ if err = internal_grow(t2, src.used * 2); err != nil { return err; }
+ if err = internal_grow(x0x0, B * 2 ); err != nil { return err; }
+ if err = internal_grow(x1x1, (src.used - B) * 2); err != nil { return err; }
+
+ /*
+ Now shift the digits.
+ */
+ x0.used = B;
+ x1.used = src.used - B;
+
+ internal_copy_digits(x0, src, x0.used);
+ #force_inline mem.copy_non_overlapping(&x1.digit[0], &src.digit[B], size_of(DIGIT) * x1.used);
+ internal_clamp(x0);
+
+ /*
+ Now calc the products x0*x0 and x1*x1.
+ */
+ if err = internal_sqr(x0x0, x0); err != nil { return err; }
+ if err = internal_sqr(x1x1, x1); err != nil { return err; }
+
+ /*
+ Now calc (x1+x0)^2
+ */
+ if err = internal_add(t1, x0, x1); err != nil { return err; }
+ if err = internal_sqr(t1, t1); err != nil { return err; }
+
+ /*
+ Add x0y0
+ */
+ if err = internal_add(t2, x0x0, x1x1); err != nil { return err; }
+ if err = internal_sub(t1, t1, t2); err != nil { return err; }
+
+ /*
+ Shift by B.
+ */
+ if err = internal_shl_digit(t1, B); err != nil { return err; }
+ if err = internal_shl_digit(x1x1, B * 2); err != nil { return err; }
+ if err = internal_add(t1, t1, x0x0); err != nil { return err; }
+ if err = internal_add(dest, t1, x1x1); err != nil { return err; }
+
+ return 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) {
|