JavaScript-QuickJS
view release on metacpan or search on metacpan
quickjs/libbf.c view on Meta::CPAN
if (get_bit(a->tab, a->len, bit_pos) != bit)
return TRUE;
bit_pos--;
n--;
}
return FALSE;
}
/* Cannot fail with BF_ST_MEM_ERROR. */
int bf_round(bf_t *r, limb_t prec, bf_flags_t flags)
{
if (r->len == 0)
return 0;
return __bf_round(r, prec, flags, r->len, 0);
}
/* for debugging */
static __maybe_unused void dump_limbs(const char *str, const limb_t *tab, limb_t n)
{
limb_t i;
printf("%s: len=%" PRId_LIMB "\n", str, n);
for(i = 0; i < n; i++) {
printf("%" PRId_LIMB ": " FMT_LIMB "\n",
i, tab[i]);
}
}
void mp_print_str(const char *str, const limb_t *tab, limb_t n)
{
slimb_t i;
printf("%s= 0x", str);
for(i = n - 1; i >= 0; i--) {
if (i != (n - 1))
printf("_");
printf(FMT_LIMB, tab[i]);
}
printf("\n");
}
static __maybe_unused void mp_print_str_h(const char *str,
const limb_t *tab, limb_t n,
limb_t high)
{
slimb_t i;
printf("%s= 0x", str);
printf(FMT_LIMB, high);
for(i = n - 1; i >= 0; i--) {
printf("_");
printf(FMT_LIMB, tab[i]);
}
printf("\n");
}
/* for debugging */
void bf_print_str(const char *str, const bf_t *a)
{
slimb_t i;
printf("%s=", str);
if (a->expn == BF_EXP_NAN) {
printf("NaN");
} else {
if (a->sign)
putchar('-');
if (a->expn == BF_EXP_ZERO) {
putchar('0');
} else if (a->expn == BF_EXP_INF) {
printf("Inf");
} else {
printf("0x0.");
for(i = a->len - 1; i >= 0; i--)
printf(FMT_LIMB, a->tab[i]);
printf("p%" PRId_LIMB, a->expn);
}
}
printf("\n");
}
/* compare the absolute value of 'a' and 'b'. Return < 0 if a < b, 0
if a = b and > 0 otherwise. */
int bf_cmpu(const bf_t *a, const bf_t *b)
{
slimb_t i;
limb_t len, v1, v2;
if (a->expn != b->expn) {
if (a->expn < b->expn)
return -1;
else
return 1;
}
len = bf_max(a->len, b->len);
for(i = len - 1; i >= 0; i--) {
v1 = get_limbz(a, a->len - len + i);
v2 = get_limbz(b, b->len - len + i);
if (v1 != v2) {
if (v1 < v2)
return -1;
else
return 1;
}
}
return 0;
}
/* Full order: -0 < 0, NaN == NaN and NaN is larger than all other numbers */
int bf_cmp_full(const bf_t *a, const bf_t *b)
{
int res;
if (a->expn == BF_EXP_NAN || b->expn == BF_EXP_NAN) {
if (a->expn == b->expn)
res = 0;
else if (a->expn == BF_EXP_NAN)
res = 1;
else
res = -1;
} else if (a->sign != b->sign) {
res = 1 - 2 * a->sign;
} else {
res = bf_cmpu(a, b);
if (a->sign)
res = -res;
}
return res;
}
/* Standard floating point comparison: return 2 if one of the operands
is NaN (unordered) or -1, 0, 1 depending on the ordering assuming
-0 == +0 */
int bf_cmp(const bf_t *a, const bf_t *b)
{
int res;
if (a->expn == BF_EXP_NAN || b->expn == BF_EXP_NAN) {
res = 2;
} else if (a->sign != b->sign) {
if (a->expn == BF_EXP_ZERO && b->expn == BF_EXP_ZERO)
res = 0;
else
res = 1 - 2 * a->sign;
} else {
res = bf_cmpu(a, b);
if (a->sign)
res = -res;
}
return res;
}
/* Compute the number of bits 'n' matching the pattern:
a= X1000..0
b= X0111..1
When computing a-b, the result will have at least n leading zero
bits.
Precondition: a > b and a.expn - b.expn = 0 or 1
*/
static limb_t count_cancelled_bits(const bf_t *a, const bf_t *b)
{
slimb_t bit_offset, b_offset, n;
int p, p1;
limb_t v1, v2, mask;
bit_offset = a->len * LIMB_BITS - 1;
b_offset = (b->len - a->len) * LIMB_BITS - (LIMB_BITS - 1) +
a->expn - b->expn;
n = 0;
/* first search the equals bits */
for(;;) {
v1 = get_limbz(a, bit_offset >> LIMB_LOG2_BITS);
v2 = get_bits(b->tab, b->len, bit_offset + b_offset);
// printf("v1=" FMT_LIMB " v2=" FMT_LIMB "\n", v1, v2);
if (v1 != v2)
break;
n += LIMB_BITS;
bit_offset -= LIMB_BITS;
}
/* find the position of the first different bit */
p = clz(v1 ^ v2) + 1;
n += p;
/* then search for '0' in a and '1' in b */
p = LIMB_BITS - p;
if (p > 0) {
/* search in the trailing p bits of v1 and v2 */
mask = limb_mask(0, p - 1);
p1 = bf_min(clz(v1 & mask), clz((~v2) & mask)) - (LIMB_BITS - p);
n += p1;
if (p1 != p)
goto done;
}
bit_offset -= LIMB_BITS;
for(;;) {
v1 = get_limbz(a, bit_offset >> LIMB_LOG2_BITS);
v2 = get_bits(b->tab, b->len, bit_offset + b_offset);
// printf("v1=" FMT_LIMB " v2=" FMT_LIMB "\n", v1, v2);
if (v1 != 0 || v2 != -1) {
/* different: count the matching bits */
p1 = bf_min(clz(v1), clz(~v2));
n += p1;
break;
}
n += LIMB_BITS;
bit_offset -= LIMB_BITS;
}
done:
return n;
}
static int bf_add_internal(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec,
bf_flags_t flags, int b_neg)
{
const bf_t *tmp;
int is_sub, ret, cmp_res, a_sign, b_sign;
a_sign = a->sign;
b_sign = b->sign ^ b_neg;
is_sub = a_sign ^ b_sign;
cmp_res = bf_cmpu(a, b);
if (cmp_res < 0) {
tmp = a;
a = b;
b = tmp;
a_sign = b_sign; /* b_sign is never used later */
}
/* abs(a) >= abs(b) */
if (cmp_res == 0 && is_sub && a->expn < BF_EXP_INF) {
/* zero result */
bf_set_zero(r, (flags & BF_RND_MASK) == BF_RNDD);
ret = 0;
} else if (a->len == 0 || b->len == 0) {
ret = 0;
if (a->expn >= BF_EXP_INF) {
if (a->expn == BF_EXP_NAN) {
/* at least one operand is NaN */
bf_set_nan(r);
} else if (b->expn == BF_EXP_INF && is_sub) {
/* infinities with different signs */
bf_set_nan(r);
ret = BF_ST_INVALID_OP;
} else {
bf_set_inf(r, a_sign);
}
} else {
/* at least one zero and not subtract */
bf_set(r, a);
r->sign = a_sign;
goto renorm;
}
} else {
slimb_t d, a_offset, b_bit_offset, i, cancelled_bits;
limb_t carry, v1, v2, u, r_len, carry1, precl, tot_len, z, sub_mask;
r->sign = a_sign;
r->expn = a->expn;
d = a->expn - b->expn;
/* must add more precision for the leading cancelled bits in
subtraction */
if (is_sub) {
if (d <= 1)
cancelled_bits = count_cancelled_bits(a, b);
else
cancelled_bits = 1;
} else {
cancelled_bits = 0;
}
/* add two extra bits for rounding */
precl = (cancelled_bits + prec + 2 + LIMB_BITS - 1) / LIMB_BITS;
tot_len = bf_max(a->len, b->len + (d + LIMB_BITS - 1) / LIMB_BITS);
r_len = bf_min(precl, tot_len);
if (bf_resize(r, r_len))
goto fail;
a_offset = a->len - r_len;
b_bit_offset = (b->len - r_len) * LIMB_BITS + d;
/* compute the bits before for the rounding */
carry = is_sub;
z = 0;
sub_mask = -is_sub;
i = r_len - tot_len;
while (i < 0) {
slimb_t ap, bp;
BOOL inflag;
ap = a_offset + i;
bp = b_bit_offset + i * LIMB_BITS;
inflag = FALSE;
if (ap >= 0 && ap < a->len) {
v1 = a->tab[ap];
inflag = TRUE;
} else {
v1 = 0;
}
if (bp + LIMB_BITS > 0 && bp < (slimb_t)(b->len * LIMB_BITS)) {
quickjs/libbf.c view on Meta::CPAN
}
if (bf_integer_to_radix(a, a1, radixl)) {
dbuf_set_error(s);
goto done;
}
radix_bits = 0;
pos = n;
pos_incr = 1;
first_buf_pos = pos * digits_per_limb - n_digits;
}
buf_pos = digits_per_limb;
i = 0;
while (i < n_digits) {
if (buf_pos == digits_per_limb) {
pos -= pos_incr;
if (radix_bits == 0) {
v = get_limbz(a, pos);
limb_to_a(buf, v, radix, digits_per_limb);
} else {
v = get_bits(a->tab, a->len, pos);
limb_to_a2(buf, v, radix_bits, digits_per_limb);
}
buf_pos = first_buf_pos;
first_buf_pos = 0;
}
if (i < dot_pos) {
l = dot_pos;
} else {
if (i == dot_pos)
dbuf_putc(s, '.');
l = n_digits;
}
l = bf_min(digits_per_limb - buf_pos, l - i);
dbuf_put(s, (uint8_t *)(buf + buf_pos), l);
buf_pos += l;
i += l;
}
done:
if (a != a1)
bf_delete(a);
}
static void *bf_dbuf_realloc(void *opaque, void *ptr, size_t size)
{
bf_context_t *s = opaque;
return bf_realloc(s, ptr, size);
}
/* return the length in bytes. A trailing '\0' is added */
static char *bf_ftoa_internal(size_t *plen, const bf_t *a2, int radix,
limb_t prec, bf_flags_t flags, BOOL is_dec)
{
bf_context_t *ctx = a2->ctx;
DynBuf s_s, *s = &s_s;
int radix_bits;
// bf_print_str("ftoa", a2);
// printf("radix=%d\n", radix);
dbuf_init2(s, ctx, bf_dbuf_realloc);
if (a2->expn == BF_EXP_NAN) {
dbuf_putstr(s, "NaN");
} else {
if (a2->sign)
dbuf_putc(s, '-');
if (a2->expn == BF_EXP_INF) {
if (flags & BF_FTOA_JS_QUIRKS)
dbuf_putstr(s, "Infinity");
else
dbuf_putstr(s, "Inf");
} else {
int fmt, ret;
slimb_t n_digits, n, i, n_max, n1;
bf_t a1_s, *a1 = &a1_s;
if ((radix & (radix - 1)) != 0)
radix_bits = 0;
else
radix_bits = ceil_log2(radix);
fmt = flags & BF_FTOA_FORMAT_MASK;
bf_init(ctx, a1);
if (fmt == BF_FTOA_FORMAT_FRAC) {
if (is_dec || radix_bits != 0) {
if (bf_set(a1, a2))
goto fail1;
#ifdef USE_BF_DEC
if (is_dec) {
if (bfdec_round((bfdec_t *)a1, prec, (flags & BF_RND_MASK) | BF_FLAG_RADPNT_PREC) & BF_ST_MEM_ERROR)
goto fail1;
n = a1->expn;
} else
#endif
{
if (bf_round(a1, prec * radix_bits, (flags & BF_RND_MASK) | BF_FLAG_RADPNT_PREC) & BF_ST_MEM_ERROR)
goto fail1;
n = ceil_div(a1->expn, radix_bits);
}
if (flags & BF_FTOA_ADD_PREFIX) {
if (radix == 16)
dbuf_putstr(s, "0x");
else if (radix == 8)
dbuf_putstr(s, "0o");
else if (radix == 2)
dbuf_putstr(s, "0b");
}
if (a1->expn == BF_EXP_ZERO) {
dbuf_putstr(s, "0");
if (prec > 0) {
dbuf_putstr(s, ".");
for(i = 0; i < prec; i++) {
dbuf_putc(s, '0');
}
}
} else {
n_digits = prec + n;
if (n <= 0) {
/* 0.x */
dbuf_putstr(s, "0.");
for(i = 0; i < -n; i++) {
dbuf_putc(s, '0');
}
quickjs/libbf.c view on Meta::CPAN
bf_context_t *s = r->ctx;
bf_t T_s, *T = &T_s;
slimb_t e, i, er;
limb_t v;
/* x = m*2^e with m odd integer */
e = bf_get_exp_min(x);
/* fast check on the exponent */
if (n > (LIMB_BITS - 1)) {
if (e != 0)
return FALSE;
er = 0;
} else {
if ((e & (((limb_t)1 << n) - 1)) != 0)
return FALSE;
er = e >> n;
}
/* every perfect odd square = 1 modulo 8 */
v = get_bits(x->tab, x->len, x->len * LIMB_BITS - x->expn + e);
if ((v & 7) != 1)
return FALSE;
bf_init(s, T);
bf_set(T, x);
T->expn -= e;
for(i = 0; i < n; i++) {
if (i != 0)
bf_set(T, r);
if (bf_sqrtrem(r, NULL, T) != 0)
return FALSE;
}
r->expn += er;
return TRUE;
}
/* prec = BF_PREC_INF is accepted for x and y integers and y >= 0 */
int bf_pow(bf_t *r, const bf_t *x, const bf_t *y, limb_t prec, bf_flags_t flags)
{
bf_context_t *s = r->ctx;
bf_t T_s, *T = &T_s;
bf_t ytmp_s;
BOOL y_is_int, y_is_odd;
int r_sign, ret, rnd_mode;
slimb_t y_emin;
if (x->len == 0 || y->len == 0) {
if (y->expn == BF_EXP_ZERO) {
/* pow(x, 0) = 1 */
bf_set_ui(r, 1);
} else if (x->expn == BF_EXP_NAN) {
bf_set_nan(r);
} else {
int cmp_x_abs_1;
bf_set_ui(r, 1);
cmp_x_abs_1 = bf_cmpu(x, r);
if (cmp_x_abs_1 == 0 && (flags & BF_POW_JS_QUIRKS) &&
(y->expn >= BF_EXP_INF)) {
bf_set_nan(r);
} else if (cmp_x_abs_1 == 0 &&
(!x->sign || y->expn != BF_EXP_NAN)) {
/* pow(1, y) = 1 even if y = NaN */
/* pow(-1, +/-inf) = 1 */
} else if (y->expn == BF_EXP_NAN) {
bf_set_nan(r);
} else if (y->expn == BF_EXP_INF) {
if (y->sign == (cmp_x_abs_1 > 0)) {
bf_set_zero(r, 0);
} else {
bf_set_inf(r, 0);
}
} else {
y_emin = bf_get_exp_min(y);
y_is_odd = (y_emin == 0);
if (y->sign == (x->expn == BF_EXP_ZERO)) {
bf_set_inf(r, y_is_odd & x->sign);
if (y->sign) {
/* pow(0, y) with y < 0 */
return BF_ST_DIVIDE_ZERO;
}
} else {
bf_set_zero(r, y_is_odd & x->sign);
}
}
}
return 0;
}
bf_init(s, T);
bf_set(T, x);
y_emin = bf_get_exp_min(y);
y_is_int = (y_emin >= 0);
rnd_mode = flags & BF_RND_MASK;
if (x->sign) {
if (!y_is_int) {
bf_set_nan(r);
bf_delete(T);
return BF_ST_INVALID_OP;
}
y_is_odd = (y_emin == 0);
r_sign = y_is_odd;
/* change the directed rounding mode if the sign of the result
is changed */
if (r_sign && (rnd_mode == BF_RNDD || rnd_mode == BF_RNDU))
flags ^= 1;
bf_neg(T);
} else {
r_sign = 0;
}
bf_set_ui(r, 1);
if (bf_cmp_eq(T, r)) {
/* abs(x) = 1: nothing more to do */
ret = 0;
} else {
/* check the overflow/underflow cases */
{
bf_t al_s, *al = &al_s;
bf_t ah_s, *ah = &ah_s;
limb_t precl = LIMB_BITS;
bf_init(s, al);
bf_init(s, ah);
quickjs/libbf.c view on Meta::CPAN
case 44: /* 17592186044416-35184372088831 */
return LIMB_DIGITS - 14;
case 45: /* 35184372088832-70368744177663 */
return LIMB_DIGITS - 14;
case 46: /* 70368744177664-140737488355327 */
if (a < 100000000000000)
return LIMB_DIGITS - 14;
else
return LIMB_DIGITS - 15;
case 47: /* 140737488355328-281474976710655 */
return LIMB_DIGITS - 15;
case 48: /* 281474976710656-562949953421311 */
return LIMB_DIGITS - 15;
case 49: /* 562949953421312-1125899906842623 */
if (a < 1000000000000000)
return LIMB_DIGITS - 15;
else
return LIMB_DIGITS - 16;
case 50: /* 1125899906842624-2251799813685247 */
return LIMB_DIGITS - 16;
case 51: /* 2251799813685248-4503599627370495 */
return LIMB_DIGITS - 16;
case 52: /* 4503599627370496-9007199254740991 */
return LIMB_DIGITS - 16;
case 53: /* 9007199254740992-18014398509481983 */
if (a < 10000000000000000)
return LIMB_DIGITS - 16;
else
return LIMB_DIGITS - 17;
case 54: /* 18014398509481984-36028797018963967 */
return LIMB_DIGITS - 17;
case 55: /* 36028797018963968-72057594037927935 */
return LIMB_DIGITS - 17;
case 56: /* 72057594037927936-144115188075855871 */
if (a < 100000000000000000)
return LIMB_DIGITS - 17;
else
return LIMB_DIGITS - 18;
case 57: /* 144115188075855872-288230376151711743 */
return LIMB_DIGITS - 18;
case 58: /* 288230376151711744-576460752303423487 */
return LIMB_DIGITS - 18;
case 59: /* 576460752303423488-1152921504606846975 */
if (a < 1000000000000000000)
return LIMB_DIGITS - 18;
else
return LIMB_DIGITS - 19;
#endif
default:
return 0;
}
}
/* for debugging */
void bfdec_print_str(const char *str, const bfdec_t *a)
{
slimb_t i;
printf("%s=", str);
if (a->expn == BF_EXP_NAN) {
printf("NaN");
} else {
if (a->sign)
putchar('-');
if (a->expn == BF_EXP_ZERO) {
putchar('0');
} else if (a->expn == BF_EXP_INF) {
printf("Inf");
} else {
printf("0.");
for(i = a->len - 1; i >= 0; i--)
printf("%0*" PRIu_LIMB, LIMB_DIGITS, a->tab[i]);
printf("e%" PRId_LIMB, a->expn);
}
}
printf("\n");
}
/* return != 0 if one digit between 0 and bit_pos inclusive is not zero. */
static inline limb_t scan_digit_nz(const bfdec_t *r, slimb_t bit_pos)
{
slimb_t pos;
limb_t v, q;
int shift;
if (bit_pos < 0)
return 0;
pos = (limb_t)bit_pos / LIMB_DIGITS;
shift = (limb_t)bit_pos % LIMB_DIGITS;
fast_shr_rem_dec(q, v, r->tab[pos], shift + 1);
(void)q;
if (v != 0)
return 1;
pos--;
while (pos >= 0) {
if (r->tab[pos] != 0)
return 1;
pos--;
}
return 0;
}
static limb_t get_digit(const limb_t *tab, limb_t len, slimb_t pos)
{
slimb_t i;
int shift;
i = floor_div(pos, LIMB_DIGITS);
if (i < 0 || i >= len)
return 0;
shift = pos - i * LIMB_DIGITS;
return fast_shr_dec(tab[i], shift) % 10;
}
#if 0
static limb_t get_digits(const limb_t *tab, limb_t len, slimb_t pos)
{
limb_t a0, a1;
int shift;
slimb_t i;
i = floor_div(pos, LIMB_DIGITS);
quickjs/libbf.c view on Meta::CPAN
r->tab[2] = v / BF_DEC_BASE;
r->expn = 3 * LIMB_DIGITS;
} else
#endif
if (v >= BF_DEC_BASE) {
if (bfdec_resize(r, 2))
goto fail;
r->tab[0] = v % BF_DEC_BASE;
r->tab[1] = v / BF_DEC_BASE;
r->expn = 2 * LIMB_DIGITS;
} else {
if (bfdec_resize(r, 1))
goto fail;
r->tab[0] = v;
r->expn = LIMB_DIGITS;
}
r->sign = 0;
return bfdec_normalize_and_round(r, BF_PREC_INF, 0);
fail:
bfdec_set_nan(r);
return BF_ST_MEM_ERROR;
}
int bfdec_set_si(bfdec_t *r, int64_t v)
{
int ret;
if (v < 0) {
ret = bfdec_set_ui(r, -v);
r->sign = 1;
} else {
ret = bfdec_set_ui(r, v);
}
return ret;
}
static int bfdec_add_internal(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, bf_flags_t flags, int b_neg)
{
bf_context_t *s = r->ctx;
int is_sub, cmp_res, a_sign, b_sign, ret;
a_sign = a->sign;
b_sign = b->sign ^ b_neg;
is_sub = a_sign ^ b_sign;
cmp_res = bfdec_cmpu(a, b);
if (cmp_res < 0) {
const bfdec_t *tmp;
tmp = a;
a = b;
b = tmp;
a_sign = b_sign; /* b_sign is never used later */
}
/* abs(a) >= abs(b) */
if (cmp_res == 0 && is_sub && a->expn < BF_EXP_INF) {
/* zero result */
bfdec_set_zero(r, (flags & BF_RND_MASK) == BF_RNDD);
ret = 0;
} else if (a->len == 0 || b->len == 0) {
ret = 0;
if (a->expn >= BF_EXP_INF) {
if (a->expn == BF_EXP_NAN) {
/* at least one operand is NaN */
bfdec_set_nan(r);
ret = 0;
} else if (b->expn == BF_EXP_INF && is_sub) {
/* infinities with different signs */
bfdec_set_nan(r);
ret = BF_ST_INVALID_OP;
} else {
bfdec_set_inf(r, a_sign);
}
} else {
/* at least one zero and not subtract */
if (bfdec_set(r, a))
return BF_ST_MEM_ERROR;
r->sign = a_sign;
goto renorm;
}
} else {
slimb_t d, a_offset, b_offset, i, r_len;
limb_t carry;
limb_t *b1_tab;
int b_shift;
mp_size_t b1_len;
d = a->expn - b->expn;
/* XXX: not efficient in time and memory if the precision is
not infinite */
r_len = bf_max(a->len, b->len + (d + LIMB_DIGITS - 1) / LIMB_DIGITS);
if (bfdec_resize(r, r_len))
goto fail;
r->sign = a_sign;
r->expn = a->expn;
a_offset = r_len - a->len;
for(i = 0; i < a_offset; i++)
r->tab[i] = 0;
for(i = 0; i < a->len; i++)
r->tab[a_offset + i] = a->tab[i];
b_shift = d % LIMB_DIGITS;
if (b_shift == 0) {
b1_len = b->len;
b1_tab = (limb_t *)b->tab;
} else {
b1_len = b->len + 1;
b1_tab = bf_malloc(s, sizeof(limb_t) * b1_len);
if (!b1_tab)
goto fail;
b1_tab[0] = mp_shr_dec(b1_tab + 1, b->tab, b->len, b_shift, 0) *
mp_pow_dec[LIMB_DIGITS - b_shift];
}
b_offset = r_len - (b->len + (d + LIMB_DIGITS - 1) / LIMB_DIGITS);
if (is_sub) {
carry = mp_sub_dec(r->tab + b_offset, r->tab + b_offset,
b1_tab, b1_len, 0);
if (carry != 0) {
carry = mp_sub_ui_dec(r->tab + b_offset + b1_len, carry,
r_len - (b_offset + b1_len));
assert(carry == 0);
( run in 1.185 second using v1.01-cache-2.11-cpan-adec679a428 )