vm: adding bignum_gcd primitive.

diff --git a/basis/stack-checker/known-words/known-words.factor b/basis/stack-checker/known-words/known-words.factor
index 99a5e7ace8432775d0a68449a35b13b8119c324c..ae760a7e47dd9d33cda2261a2cc5ed244ead9b4d 100644 (file)
--- 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 fc9e3670c7924874ccaee7882410f7b897aa40cd..c5dd0ee3085c00b0c807688c82146a43ed7ea8a4 100755 (executable)
--- 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 93bdaf8a1dd9c07743d69335dddef1a6229e1519..87a483db5513235295be3f3780096dc6296c215d 100755 (executable)
--- 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 ff0d86a681076f835d029bc7e0cc1061b7534a26..e8e5ce011682bfb42d4c5bb30b083d0b6ab0b6eb 100644 (file)
--- 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 ce8ac6eaf4bdf722ddddecaea8904f7c61891bbb..45218aff01aba703b0de3a72c71c529f702b428a 100755 (executable)
--- a/vm/math.cpp
+++ b/vm/math.cpp
@@ -147,6 +147,12 @@ void factor_vm::primitive_bignum_mod()
        ctx->push(tag<bignum>(bignum_remainder(x,y)));
 }
 
+void factor_vm::primitive_bignum_gcd()
+{
+       POP_BIGNUMS(x,y);
+       ctx->push(tag<bignum>(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 7884cf2e30c9442aebb3360793debbbe0e8b6474..226cc3541af9e54427a0c0334bedf4fe30d13e1c 100644 (file)
--- 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 01f0a2afc209b66ff5aafd504b60b4c7578cc6be..112b65b3413a5dda83050ae2a86ab7d5c5521bfe 100755 (executable)
--- 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();