diff --git a/packages/runtime/src/scr_bigint.c b/packages/runtime/src/scr_bigint.c index dd0a70ff7..5af2eafaa 100644 --- a/packages/runtime/src/scr_bigint.c +++ b/packages/runtime/src/scr_bigint.c @@ -215,30 +215,81 @@ static int bi_cmp_limbs(const bi_limb *a, size_t an, const bi_limb *b, size_t bn return 0; } -/* r = a + b over n limbs; returns the carry. r may alias a or b. */ +/* x + y + carry_in (0 or 1), with the carry out in *carry: the add-with- + * carry builtins compile to one adc chain. */ +#if defined(__has_builtin) +#if __has_builtin(__builtin_addcll) && __has_builtin(__builtin_subcll) && \ + __has_builtin(__builtin_addc) && __has_builtin(__builtin_subc) +#define BI_HAVE_ADDC 1 +#endif +#endif +static inline bi_limb bi_addc(bi_limb x, bi_limb y, bi_limb *carry) { +#if defined(BI_HAVE_ADDC) && BI_LIMB_BITS == 64 + unsigned long long c; + bi_limb s = __builtin_addcll(x, y, *carry, &c); + *carry = c; + return s; +#elif defined(BI_HAVE_ADDC) + unsigned c; + bi_limb s = __builtin_addc(x, y, *carry, &c); + *carry = c; + return s; +#else + bi_limb s = x + y; + bi_limb c1 = s < x; + bi_limb t = s + *carry; + *carry = c1 | (t < s); + return t; +#endif +} + +/* x - y - borrow_in (0 or 1), with the borrow out in *borrow. */ +static inline bi_limb bi_subc(bi_limb x, bi_limb y, bi_limb *borrow) { +#if defined(BI_HAVE_ADDC) && BI_LIMB_BITS == 64 + unsigned long long b; + bi_limb d = __builtin_subcll(x, y, *borrow, &b); + *borrow = b; + return d; +#elif defined(BI_HAVE_ADDC) + unsigned b; + bi_limb d = __builtin_subc(x, y, *borrow, &b); + *borrow = b; + return d; +#else + bi_limb d = x - y; + bi_limb b1 = x < y; + bi_limb t = d - *borrow; + *borrow = b1 | (d < *borrow); + return t; +#endif +} + +/* r = a + b over n limbs; returns the carry. r may alias a or b. Four + * limbs a pass keep the carry in one chain. */ static bi_limb bi_add_n(bi_limb *r, const bi_limb *a, const bi_limb *b, size_t n) { bi_limb carry = 0; - for (size_t i = 0; i < n; i++) { - bi_limb x = a[i], s = x + b[i]; - bi_limb c1 = s < x; - bi_limb t = s + carry; - carry = c1 | (t < s); - r[i] = t; + size_t i = 0; + for (; i + 4 <= n; i += 4) { + r[i] = bi_addc(a[i], b[i], &carry); + r[i + 1] = bi_addc(a[i + 1], b[i + 1], &carry); + r[i + 2] = bi_addc(a[i + 2], b[i + 2], &carry); + r[i + 3] = bi_addc(a[i + 3], b[i + 3], &carry); } + for (; i < n; i++) r[i] = bi_addc(a[i], b[i], &carry); return carry; } /* r = a - b over n limbs; returns the borrow. r may alias a or b. */ static bi_limb bi_sub_n(bi_limb *r, const bi_limb *a, const bi_limb *b, size_t n) { bi_limb borrow = 0; - for (size_t i = 0; i < n; i++) { - bi_limb x = a[i], y = b[i]; - bi_limb d = x - y; - bi_limb b1 = x < y; - bi_limb t = d - borrow; - borrow = b1 | (d < borrow); - r[i] = t; + size_t i = 0; + for (; i + 4 <= n; i += 4) { + r[i] = bi_subc(a[i], b[i], &borrow); + r[i + 1] = bi_subc(a[i + 1], b[i + 1], &borrow); + r[i + 2] = bi_subc(a[i + 2], b[i + 2], &borrow); + r[i + 3] = bi_subc(a[i + 3], b[i + 3], &borrow); } + for (; i < n; i++) r[i] = bi_subc(a[i], b[i], &borrow); return borrow; } @@ -277,9 +328,22 @@ static bi_limb bi_sub_into(bi_limb *r, size_t rn, const bi_limb *a, size_t an) { return borrow; } -/* r = a * m + carry over n limbs; returns the high limb. r may alias a. */ +/* r = a * m + carry over n limbs; returns the high limb. r may alias a. + * Four limbs a pass: the four products, then one carry chain adding each + * product's high limb into the next low limb. */ static bi_limb bi_mul1(bi_limb *r, const bi_limb *a, size_t n, bi_limb m, bi_limb carry) { - for (size_t i = 0; i < n; i++) { + size_t i = 0; + for (; i + 4 <= n; i += 4) { + bi_wide p0 = (bi_wide)a[i] * m, p1 = (bi_wide)a[i + 1] * m; + bi_wide p2 = (bi_wide)a[i + 2] * m, p3 = (bi_wide)a[i + 3] * m; + bi_limb c = 0; + r[i] = bi_addc((bi_limb)p0, carry, &c); + r[i + 1] = bi_addc((bi_limb)p1, (bi_limb)(p0 >> BI_LIMB_BITS), &c); + r[i + 2] = bi_addc((bi_limb)p2, (bi_limb)(p1 >> BI_LIMB_BITS), &c); + r[i + 3] = bi_addc((bi_limb)p3, (bi_limb)(p2 >> BI_LIMB_BITS), &c); + carry = (bi_limb)(p3 >> BI_LIMB_BITS) + c; + } + for (; i < n; i++) { bi_wide p = (bi_wide)a[i] * m + carry; r[i] = (bi_limb)p; carry = (bi_limb)(p >> BI_LIMB_BITS); @@ -287,10 +351,27 @@ static bi_limb bi_mul1(bi_limb *r, const bi_limb *a, size_t n, bi_limb m, bi_lim return carry; } -/* r += a * m over n limbs; returns the high limb. */ +/* r += a * m over n limbs; returns the high limb. Four limbs a pass: one + * carry chain forms the four product limbs, a second adds them to r. */ static bi_limb bi_addmul1(bi_limb *r, const bi_limb *a, size_t n, bi_limb m) { bi_limb carry = 0; - for (size_t i = 0; i < n; i++) { + size_t i = 0; + for (; i + 4 <= n; i += 4) { + bi_wide p0 = (bi_wide)a[i] * m, p1 = (bi_wide)a[i + 1] * m; + bi_wide p2 = (bi_wide)a[i + 2] * m, p3 = (bi_wide)a[i + 3] * m; + bi_limb c = 0, k = 0; + bi_limb t0 = bi_addc((bi_limb)p0, carry, &c); + bi_limb t1 = bi_addc((bi_limb)p1, (bi_limb)(p0 >> BI_LIMB_BITS), &c); + bi_limb t2 = bi_addc((bi_limb)p2, (bi_limb)(p1 >> BI_LIMB_BITS), &c); + bi_limb t3 = bi_addc((bi_limb)p3, (bi_limb)(p2 >> BI_LIMB_BITS), &c); + r[i] = bi_addc(r[i], t0, &k); + r[i + 1] = bi_addc(r[i + 1], t1, &k); + r[i + 2] = bi_addc(r[i + 2], t2, &k); + r[i + 3] = bi_addc(r[i + 3], t3, &k); + /* r + a * m over these limbs stays below B^5, so this cannot wrap */ + carry = (bi_limb)(p3 >> BI_LIMB_BITS) + c + k; + } + for (; i < n; i++) { bi_wide p = (bi_wide)a[i] * m + r[i] + carry; r[i] = (bi_limb)p; carry = (bi_limb)(p >> BI_LIMB_BITS); @@ -301,7 +382,23 @@ static bi_limb bi_addmul1(bi_limb *r, const bi_limb *a, size_t n, bi_limb m) { /* r -= a * m over n limbs; returns the borrow limb. */ static bi_limb bi_submul1(bi_limb *r, const bi_limb *a, size_t n, bi_limb m) { bi_limb carry = 0; - for (size_t i = 0; i < n; i++) { + size_t i = 0; + for (; i + 4 <= n; i += 4) { + bi_wide p0 = (bi_wide)a[i] * m, p1 = (bi_wide)a[i + 1] * m; + bi_wide p2 = (bi_wide)a[i + 2] * m, p3 = (bi_wide)a[i + 3] * m; + bi_limb c = 0, k = 0; + bi_limb t0 = bi_addc((bi_limb)p0, carry, &c); + bi_limb t1 = bi_addc((bi_limb)p1, (bi_limb)(p0 >> BI_LIMB_BITS), &c); + bi_limb t2 = bi_addc((bi_limb)p2, (bi_limb)(p1 >> BI_LIMB_BITS), &c); + bi_limb t3 = bi_addc((bi_limb)p3, (bi_limb)(p2 >> BI_LIMB_BITS), &c); + r[i] = bi_subc(r[i], t0, &k); + r[i + 1] = bi_subc(r[i + 1], t1, &k); + r[i + 2] = bi_subc(r[i + 2], t2, &k); + r[i + 3] = bi_subc(r[i + 3], t3, &k); + /* a * m + carry over these limbs stays below B^5: the borrow fits */ + carry = (bi_limb)(p3 >> BI_LIMB_BITS) + c + k; + } + for (; i < n; i++) { bi_wide p = (bi_wide)a[i] * m + carry; bi_limb lo = (bi_limb)p; carry = (bi_limb)(p >> BI_LIMB_BITS); @@ -543,6 +640,56 @@ static void bi_divmod_knuth(bi_limb *q, bi_limb *r, const bi_limb *a, size_t an, free(un); } +/* a / b for a and b of n >= 2 limbs each, b's top limb nonzero: the + * quotient is one limb. Knuth's estimate comes from the top limbs as + * they read normalized (three of a, two of b) and is the quotient or one + * more, which one pass of a - qhat * b over the operands as they are + * settles, with no normalized copies. r, unless NULL, receives the n + * limbs of a % b. */ +static bi_limb bi_div_qlimb(bi_limb *r, const bi_limb *a, const bi_limb *b, size_t n) { + unsigned s = bi_clz(b[n - 1]); + bi_limb d1 = b[n - 1], d0 = b[n - 2], u2 = 0, u1 = a[n - 1], u0 = a[n - 2]; + if (s) { + unsigned t = BI_LIMB_BITS - s; + bi_limb b3 = n >= 3 ? b[n - 3] : 0, a3 = n >= 3 ? a[n - 3] : 0; + d1 = (d1 << s) | (d0 >> t); + d0 = (d0 << s) | (b3 >> t); + u2 = u1 >> t; + u1 = (u1 << s) | (u0 >> t); + u0 = (u0 << s) | (a3 >> t); + } + bi_limb qhat, rhat; + bool rhat_big; /* rhat >= B: the estimate needs no further test */ + if (u2 >= d1) { + qhat = BI_LIMB_MAX; + rhat = u1 + d1; + rhat_big = rhat < d1; + } else { + bi_divisor top = bi_divisor_of(d1); + qhat = bi_div21(u2, u1, top.d, top.v, &rhat); + rhat_big = false; + } + while (!rhat_big && (bi_wide)qhat * d0 > (((bi_wide)rhat << BI_LIMB_BITS) | u0)) { + qhat--; + rhat += d1; + rhat_big = rhat < d1; + } + /* a - qhat * b: negative when the product's high limb or the borrow + * reaches past a's n limbs */ + bi_limb carry = 0, borrow = 0; + for (size_t i = 0; i < n; i++) { + bi_wide p = (bi_wide)b[i] * qhat + carry; + carry = (bi_limb)(p >> BI_LIMB_BITS); + bi_limb d = bi_subc(a[i], (bi_limb)p, &borrow); + if (r) r[i] = d; + } + if (carry || borrow) { + qhat--; + if (r) bi_add_n(r, r, b, n); + } + return qhat; +} + /* ---- comparisons and hashing ----------------------------------------- */ static int bi_cmp_abs(const ScrBigInt *a, const ScrBigInt *b) { @@ -709,6 +856,13 @@ static void bi_divmod_abs(const ScrBigInt *a, const ScrBigInt *b, ScrBigInt **qo if (rout) *rout = bi_u64_signed(rem, 1); return; } + if (an == bn) { + ScrBigInt *r = rout ? bi_alloc(bn) : NULL; + bi_limb q = bi_div_qlimb(r ? r->limb : NULL, a->limb, b->limb, bn); + if (qout) *qout = bi_u64_signed(q, 1); + if (rout) *rout = bi_finish(r, bn, 1); + return; + } size_t qn = an - bn + 1; ScrBigInt *q = bi_alloc(qn); ScrBigInt *r = rout ? bi_alloc(bn) : NULL;