vm: adding bignum_gcd primitive.

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

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 ) }

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);
+}
+
 }

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))
 

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);

vm/primitives.hpp

@@ -21,6 +21,7 @@ namespace factor
 	_(bignum_lesseq) \
 	_(bignum_log2) \
 	_(bignum_mod) \
+	_(bignum_gcd) \
 	_(bignum_multiply) \
 	_(bignum_not) \
 	_(bignum_or) \

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();