/* Allocates memory */
bignum * factor_vm::bignum_gcd(bignum * a, bignum * b)
{
- GC_BIGNUM(a);
- GC_BIGNUM(b);
- bignum * d;
- bignum_length_type size_a, size_b;
- bignum_digit_type * scan_a, * scan_b, * scan_d, * a_end, * b_end;
-
- if (BIGNUM_NEGATIVE_P (a)) {
- size_a = BIGNUM_LENGTH (a);
- d = allot_bignum (size_a, 0);
- scan_d = BIGNUM_START_PTR (d);
- scan_a = BIGNUM_START_PTR (a);
- a_end = scan_a + size_a;
- while (scan_a < a_end)
- (*scan_d++) = (*scan_a++);
- a = d;
- }
-
- if (BIGNUM_NEGATIVE_P (b)) {
- size_b = BIGNUM_LENGTH (b);
- d = allot_bignum (size_b, 0);
- scan_d = BIGNUM_START_PTR (d);
- scan_b = BIGNUM_START_PTR (b);
- b_end = scan_b + size_b;
- while (scan_b < b_end)
- (*scan_d++) = (*scan_b++);
- b = d;
- }
-
- if (bignum_compare(a, b) == bignum_comparison_less) {
- std::swap(a, b);
- }
-
- while (BIGNUM_LENGTH (b) != 0) {
- d = bignum_remainder (a, b);
- GC_BIGNUM(d);
- if (d == BIGNUM_OUT_OF_BAND) {
- return d;
- }
- a = b;
- b = d;
- }
-
- return a;
+ GC_BIGNUM(a);
+ GC_BIGNUM(b);
+ bignum * d;
+ bignum_length_type size_a, size_b;
+ bignum_digit_type * scan_a, * scan_b, * scan_d, * a_end, * b_end;
+
+ if (BIGNUM_NEGATIVE_P (a)) {
+ size_a = BIGNUM_LENGTH (a);
+ d = allot_bignum (size_a, 0);
+ scan_d = BIGNUM_START_PTR (d);
+ scan_a = BIGNUM_START_PTR (a);
+ a_end = scan_a + size_a;
+ while (scan_a < a_end)
+ (*scan_d++) = (*scan_a++);
+ a = d;
+ }
+
+ if (BIGNUM_NEGATIVE_P (b)) {
+ size_b = BIGNUM_LENGTH (b);
+ d = allot_bignum (size_b, 0);
+ scan_d = BIGNUM_START_PTR (d);
+ scan_b = BIGNUM_START_PTR (b);
+ b_end = scan_b + size_b;
+ while (scan_b < b_end)
+ (*scan_d++) = (*scan_b++);
+ b = d;
+ }
+
+ if (bignum_compare(a, b) == bignum_comparison_less) {
+ std::swap(a, b);
+ }
+
+ while (BIGNUM_LENGTH (b) != 0) {
+ d = bignum_remainder (a, b);
+ GC_BIGNUM(d);
+ if (d == BIGNUM_OUT_OF_BAND) {
+ return d;
+ }
+ a = b;
+ b = d;
+ }
+
+ return a;
}
#else
/* Allocates memory */
bignum * factor_vm::bignum_gcd(bignum * a, bignum * b)
{
- GC_BIGNUM(a);
- GC_BIGNUM(b);
- bignum *c, *d, *e, *f;
- bignum_twodigit_type x, y, q, s, t, A, B, C, D;
- int nbits, k;
- bignum_length_type size_a, size_b, size_c;
- bignum_digit_type *scan_a, *scan_b, *scan_c, *scan_d;
- bignum_digit_type *a_end, *b_end, *c_end;
-
- /* clone the bignums so we can modify them in-place */
- size_a = BIGNUM_LENGTH (a);
- c = allot_bignum (size_a, 0);
- scan_a = BIGNUM_START_PTR (a);
- a_end = scan_a + size_a;
- GC_BIGNUM(c);
- scan_c = BIGNUM_START_PTR (c);
- while (scan_a < a_end)
- (*scan_c++) = (*scan_a++);
- a = c;
- size_b = BIGNUM_LENGTH (b);
- d = allot_bignum (size_b, 0);
- scan_b = BIGNUM_START_PTR (b);
- b_end = scan_b + size_b;
- GC_BIGNUM(d);
- 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. */
- if (bignum_compare(a, b) == bignum_comparison_less) {
- std::swap(a, b);
- std::swap(size_a, size_b);
- }
-
- while (size_a > 1) {
- nbits = log2 (BIGNUM_REF (a, size_a-1));
- 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));
+ GC_BIGNUM(a);
+ GC_BIGNUM(b);
+ bignum *c, *d, *e, *f;
+ bignum_twodigit_type x, y, q, s, t, A, B, C, D;
+ int nbits, k;
+ bignum_length_type size_a, size_b, size_c;
+ bignum_digit_type *scan_a, *scan_b, *scan_c, *scan_d;
+ bignum_digit_type *a_end, *b_end, *c_end;
+
+ /* clone the bignums so we can modify them in-place */
+ size_a = BIGNUM_LENGTH (a);
+ c = allot_bignum (size_a, 0);
+ scan_a = BIGNUM_START_PTR (a);
+ a_end = scan_a + size_a;
+ GC_BIGNUM(c);
+ scan_c = BIGNUM_START_PTR (c);
+ while (scan_a < a_end)
+ (*scan_c++) = (*scan_a++);
+ a = c;
+ size_b = BIGNUM_LENGTH (b);
+ d = allot_bignum (size_b, 0);
+ scan_b = BIGNUM_START_PTR (b);
+ b_end = scan_b + size_b;
+ GC_BIGNUM(d);
+ 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. */
+ if (bignum_compare(a, b) == bignum_comparison_less) {
+ std::swap(a, b);
+ std::swap(size_a, size_b);
+ }
+
+ while (size_a > 1) {
+ nbits = log2 (BIGNUM_REF (a, size_a-1));
+ 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 (size_b == 0) {
- return bignum_trim (a);
- }
- e = bignum_trim (a);
- GC_BIGNUM(e);
- f = bignum_trim (b);
- GC_BIGNUM(f);
- c = bignum_remainder (e, f);
- GC_BIGNUM (c);
- if (c == BIGNUM_OUT_OF_BAND) {
- return c;
- }
-
- // copy 'b' to 'a'
- scan_a = BIGNUM_START_PTR (a);
- scan_b = BIGNUM_START_PTR (b);
- a_end = scan_a + size_a;
- b_end = scan_b + size_b;
- while (scan_b < b_end) *(scan_a++) = *(scan_b++);
- while (scan_a < a_end) *(scan_a++) = 0;
- size_a = size_b;
-
- // copy 'c' to 'b'
- scan_b = BIGNUM_START_PTR (b);
- scan_c = BIGNUM_START_PTR (c);
- size_c = BIGNUM_LENGTH (c);
- c_end = scan_c + size_c;
- while (scan_c < c_end) *(scan_b++) = *(scan_c++);
- while (scan_b < b_end) *(scan_b++) = 0;
- size_b = size_c;
-
- 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 = scan_a;
- scan_d = scan_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);
-
- // update size_a and size_b to remove any zeros at end
- while (size_a > 0 && *(--scan_a) == 0) size_a--;
- while (size_b > 0 && *(--scan_b) == 0) size_b--;
-
- BIGNUM_ASSERT (size_a >= size_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);
+ 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 (size_b == 0) {
+ return bignum_trim (a);
+ }
+ e = bignum_trim (a);
+ GC_BIGNUM(e);
+ f = bignum_trim (b);
+ GC_BIGNUM(f);
+ c = bignum_remainder (e, f);
+ GC_BIGNUM (c);
+ if (c == BIGNUM_OUT_OF_BAND) {
+ return c;
+ }
+
+ // copy 'b' to 'a'
+ scan_a = BIGNUM_START_PTR (a);
+ scan_b = BIGNUM_START_PTR (b);
+ a_end = scan_a + size_a;
+ b_end = scan_b + size_b;
+ while (scan_b < b_end) *(scan_a++) = *(scan_b++);
+ while (scan_a < a_end) *(scan_a++) = 0;
+ size_a = size_b;
+
+ // copy 'c' to 'b'
+ scan_b = BIGNUM_START_PTR (b);
+ scan_c = BIGNUM_START_PTR (c);
+ size_c = BIGNUM_LENGTH (c);
+ c_end = scan_c + size_c;
+ while (scan_c < c_end) *(scan_b++) = *(scan_c++);
+ while (scan_b < b_end) *(scan_b++) = 0;
+ size_b = size_c;
+
+ 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 = scan_a;
+ scan_d = scan_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);
+
+ // update size_a and size_b to remove any zeros at end
+ while (size_a > 0 && *(--scan_a) == 0) size_a--;
+ while (size_b > 0 && *(--scan_b) == 0) size_b--;
+
+ BIGNUM_ASSERT (size_a >= size_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);
}
#endif
projects / factor.git / commitdiff