From 22c26ff3f5959d4135898b439154ab8537d45335 Mon Sep 17 00:00:00 2001 From: John Benediktsson Date: Thu, 5 Apr 2012 09:17:35 -0700 Subject: [PATCH] vm: adding bignum_gcd primitive. --- .../known-words/known-words.factor | 1 + core/bootstrap/primitives.factor | 1 + vm/bignum.cpp | 147 ++++++++++++++++++ vm/bignumint.hpp | 6 + vm/math.cpp | 6 + vm/primitives.hpp | 1 + vm/vm.hpp | 2 + 7 files changed, 164 insertions(+) diff --git a/basis/stack-checker/known-words/known-words.factor b/basis/stack-checker/known-words/known-words.factor index 99a5e7ace8..ae760a7e47 100644 --- a/basis/stack-checker/known-words/known-words.factor +++ b/basis/stack-checker/known-words/known-words.factor @@ -333,6 +333,7 @@ M: object infer-call* \ call bad-macro-input ; \ bignum-bitxor { bignum bignum } { bignum } define-primitive \ bignum-bitxor make-foldable \ bignum-log2 { bignum } { bignum } define-primitive \ bignum-log2 make-foldable \ bignum-mod { bignum bignum } { bignum } define-primitive \ bignum-mod make-foldable +\ bignum-gcd { bignum bignum } { bignum } define-primitive \ bignum-gcd make-foldable \ bignum-shift { bignum fixnum } { bignum } define-primitive \ bignum-shift make-foldable \ bignum/i { bignum bignum } { bignum } define-primitive \ bignum/i make-foldable \ bignum/mod { bignum bignum } { bignum bignum } define-primitive \ bignum/mod make-foldable diff --git a/core/bootstrap/primitives.factor b/core/bootstrap/primitives.factor index fc9e3670c7..c5dd0ee308 100755 --- a/core/bootstrap/primitives.factor +++ b/core/bootstrap/primitives.factor @@ -491,6 +491,7 @@ tuple { "bignum-bitxor" "math.private" "primitive_bignum_xor" ( x y -- z ) } { "bignum-log2" "math.private" "primitive_bignum_log2" ( x -- n ) } { "bignum-mod" "math.private" "primitive_bignum_mod" ( x y -- z ) } + { "bignum-gcd" "math.private" "primitive_bignum_gcd" ( x y -- z ) } { "bignum-shift" "math.private" "primitive_bignum_shift" ( x y -- z ) } { "bignum/i" "math.private" "primitive_bignum_divint" ( x y -- z ) } { "bignum/mod" "math.private" "primitive_bignum_divmod" ( x y -- z w ) } diff --git a/vm/bignum.cpp b/vm/bignum.cpp index 93bdaf8a1d..87a483db55 100755 --- a/vm/bignum.cpp +++ b/vm/bignum.cpp @@ -1709,4 +1709,151 @@ int factor_vm::bignum_unsigned_logbitp(int shift, bignum * bignum) return (digit & mask) ? 1 : 0; } +# + +/* Allocates memory */ +bignum * factor_vm::bignum_gcd(bignum * a, bignum * b) +{ + bignum * d; + bignum_digit_type x, y; + bignum_twodigit_type q, s, t, A, B, C, D; + int nbits, k; + bignum_length_type size_a, size_b; + bignum_digit_type * scan_a, * scan_b, * scan_c, * scan_d, * a_end, * b_end; + + /* clone the bignums so we can modify them in-place */ + scan_a = BIGNUM_START_PTR (a); + size_a = BIGNUM_LENGTH (a); + a_end = scan_a + size_a; + d = allot_bignum (size_a, 0); + scan_d = BIGNUM_START_PTR (d); + while (scan_a < a_end) + (*scan_d++) = (*scan_a++); + a = d; + scan_b = BIGNUM_START_PTR (b); + size_b = BIGNUM_LENGTH (b); + b_end = scan_b + size_b; + d = allot_bignum (size_b, 0); + scan_d = BIGNUM_START_PTR (d); + while (scan_b < b_end) + (*scan_d++) = (*scan_b++); + b = d; + + /* Initial reduction: make sure that 0 <= b <= a. */ + BIGNUM_SET_NEGATIVE_P (a, 0); + BIGNUM_SET_NEGATIVE_P (b, 0); + if (bignum_compare(a, b) == bignum_comparison_less) { + d = a; + a = b; + b = d; + } + + while ((size_a = BIGNUM_LENGTH (a)) > 1) { + nbits = log2 (BIGNUM_REF (a, size_a-1)); + size_b = BIGNUM_LENGTH (b); + x = ((BIGNUM_REF (a, size_a-1) << (BIGNUM_DIGIT_LENGTH - nbits)) | + (BIGNUM_REF (a, size_a-2) >> nbits)); + y = ((size_b >= size_a - 1 ? BIGNUM_REF (b, size_a-2) >> nbits : 0) | + (size_b >= size_a ? BIGNUM_REF (b, size_a-1) << (BIGNUM_DIGIT_LENGTH - nbits) : 0)); + + /* inner loop of Lehmer's algorithm; */ + A = 1; B = 0; C = 0; D = 1; + for (k=0 ;; k++) { + if (y - C == 0) + break; + q = (x + (A - 1)) / (y - C); + s = B + (q * D); + t = x - (q * y); + if (s > t) + break; + x = y; + y = t; + t = A + (q * C); + A = D; B = C; C = s; D = t; + } + + if (k == 0) { + /* no progress; do a Euclidean step */ + if (BIGNUM_LENGTH (b) == 0) { + return a; + } + d = bignum_remainder (a, b); + if (d == BIGNUM_OUT_OF_BAND) { + return d; + } + a = b; + b = d; + continue; + } + + /* + a, b = A*b - B*a, D*a - C*b if k is odd + a, b = A*a - B*b, D*b - C*a if k is even + */ + scan_a = BIGNUM_START_PTR (a); + scan_b = BIGNUM_START_PTR (b); + scan_c = BIGNUM_START_PTR (a); + scan_d = BIGNUM_START_PTR (b); + a_end = scan_a + size_a; + b_end = scan_b + size_b; + s = 0; + t = 0; + if (k & 1) { + while (scan_b < b_end) { + s += (A * *scan_b) - (B * *scan_a); + t += (D * *scan_a++) - (C * *scan_b++); + *scan_c++ = (bignum_digit_type) (s & BIGNUM_DIGIT_MASK); + *scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); + s >>= BIGNUM_DIGIT_LENGTH; + t >>= BIGNUM_DIGIT_LENGTH; + } + while (scan_a < a_end) { + s -= (B * *scan_a); + t += (D * *scan_a++); + *scan_c++ = (bignum_digit_type) (s & BIGNUM_DIGIT_MASK); + *scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); + s >>= BIGNUM_DIGIT_LENGTH; + t >>= BIGNUM_DIGIT_LENGTH; + } + } + else { + while (scan_b < b_end) { + s += (A * *scan_a) - (B * *scan_b); + t += (D * *scan_b++) - (C * *scan_a++); + *scan_c++ = (bignum_digit_type) (s & BIGNUM_DIGIT_MASK); + *scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); + s >>= BIGNUM_DIGIT_LENGTH; + t >>= BIGNUM_DIGIT_LENGTH; + } + while (scan_a < a_end) { + s += (A * *scan_a); + t -= (C * *scan_a++); + *scan_c++ = (bignum_digit_type) (s & BIGNUM_DIGIT_MASK); + *scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK); + s >>= BIGNUM_DIGIT_LENGTH; + t >>= BIGNUM_DIGIT_LENGTH; + } + } + BIGNUM_ASSERT (s == 0); + BIGNUM_ASSERT (t == 0); + + a = bignum_trim (a); + b = bignum_trim (b); + } + + /* a fits into a fixnum, so b must too */ + fixnum xx = bignum_to_fixnum (a); + fixnum yy = bignum_to_fixnum (b); + fixnum tt; + + /* usual Euclidean algorithm for longs */ + while (yy != 0) { + tt = yy; + yy = xx % yy; + xx = tt; + } + + return fixnum_to_bignum (xx); +} + } diff --git a/vm/bignumint.hpp b/vm/bignumint.hpp index ff0d86a681..e8e5ce0116 100644 --- a/vm/bignumint.hpp +++ b/vm/bignumint.hpp @@ -48,6 +48,12 @@ namespace factor typedef fixnum bignum_digit_type; typedef fixnum bignum_length_type; +#ifdef FACTOR_64 +typedef __int128_t bignum_twodigit_type; +#else +typedef s64 bignum_twodigit_type; +#endif + /* BIGNUM_TO_POINTER casts a bignum object to a digit array pointer. */ #define BIGNUM_TO_POINTER(bignum) ((bignum_digit_type *)(bignum + 1)) diff --git a/vm/math.cpp b/vm/math.cpp index ce8ac6eaf4..45218aff01 100755 --- a/vm/math.cpp +++ b/vm/math.cpp @@ -147,6 +147,12 @@ void factor_vm::primitive_bignum_mod() ctx->push(tag(bignum_remainder(x,y))); } +void factor_vm::primitive_bignum_gcd() +{ + POP_BIGNUMS(x,y); + ctx->push(tag(bignum_gcd(x,y))); +} + void factor_vm::primitive_bignum_and() { POP_BIGNUMS(x,y); diff --git a/vm/primitives.hpp b/vm/primitives.hpp index 7884cf2e30..226cc3541a 100644 --- a/vm/primitives.hpp +++ b/vm/primitives.hpp @@ -21,6 +21,7 @@ namespace factor _(bignum_lesseq) \ _(bignum_log2) \ _(bignum_mod) \ + _(bignum_gcd) \ _(bignum_multiply) \ _(bignum_not) \ _(bignum_or) \ diff --git a/vm/vm.hpp b/vm/vm.hpp index 01f0a2afc2..112b65b341 100755 --- a/vm/vm.hpp +++ b/vm/vm.hpp @@ -291,6 +291,7 @@ struct factor_vm bignum *bignum_integer_length(bignum * x); int bignum_logbitp(int shift, bignum * arg); int bignum_unsigned_logbitp(int shift, bignum * bignum); + bignum *bignum_gcd(bignum * a, bignum * b); //data heap void init_card_decks(); @@ -489,6 +490,7 @@ struct factor_vm void primitive_bignum_divint(); void primitive_bignum_divmod(); void primitive_bignum_mod(); + void primitive_bignum_gcd(); void primitive_bignum_and(); void primitive_bignum_or(); void primitive_bignum_xor();