in GaiaXAndroidQuickJS/quickjs/gxquickjs/libbf.c [4656:4824]
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);
/* compute bounds of log(abs(x)) * y with a low precision */
/* XXX: compute bf_log() once */
/* XXX: add a fast test before this slow test */
bf_log(al, T, precl, BF_RNDD);
bf_log(ah, T, precl, BF_RNDU);
bf_mul(al, al, y, precl, BF_RNDD ^ y->sign);
bf_mul(ah, ah, y, precl, BF_RNDU ^ y->sign);
ret = check_exp_underflow_overflow(s, r, al, ah, prec, flags);
bf_delete(al);
bf_delete(ah);
if (ret)
goto done;
}
if (y_is_int) {
slimb_t T_bits, e;
int_pow:
T_bits = T->expn - bf_get_exp_min(T);
if (T_bits == 1) {
/* pow(2^b, y) = 2^(b*y) */
bf_mul_si(T, y, T->expn - 1, LIMB_BITS, BF_RNDZ);
bf_get_limb(&e, T, 0);
bf_set_ui(r, 1);
ret = bf_mul_2exp(r, e, prec, flags);
} else if (prec == BF_PREC_INF) {
slimb_t y1;
/* specific case for infinite precision (integer case) */
bf_get_limb(&y1, y, 0);
assert(!y->sign);
/* x must be an integer, so abs(x) >= 2 */
if (y1 >= ((slimb_t)1 << BF_EXP_BITS_MAX)) {
bf_delete(T);
return bf_set_overflow(r, 0, BF_PREC_INF, flags);
}
ret = bf_pow_ui(r, T, y1, BF_PREC_INF, BF_RNDZ);
} else {
if (y->expn <= 31) {
/* small enough power: use exponentiation in all cases */
} else if (y->sign) {
/* cannot be exact */
goto general_case;
} else {
if (rnd_mode == BF_RNDF)
goto general_case; /* no need to track exact results */
/* see if the result has a chance to be exact:
if x=a*2^b (a odd), x^y=a^y*2^(b*y)
x^y needs a precision of at least floor_log2(a)*y bits
*/
bf_mul_si(r, y, T_bits - 1, LIMB_BITS, BF_RNDZ);
bf_get_limb(&e, r, 0);
if (prec < e)
goto general_case;
}
ret = bf_ziv_rounding(r, T, prec, flags, bf_pow_int, (void *)y);
}
} else {
if (rnd_mode != BF_RNDF) {
bf_t *y1;
if (y_emin < 0 && check_exact_power2n(r, T, -y_emin)) {
/* the problem is reduced to a power to an integer */
#if 0
printf("\nn=%" PRId64 "\n", -(int64_t)y_emin);
bf_print_str("T", T);
bf_print_str("r", r);
#endif
bf_set(T, r);
y1 = &ytmp_s;
y1->tab = y->tab;
y1->len = y->len;
y1->sign = y->sign;
y1->expn = y->expn - y_emin;
y = y1;
goto int_pow;
}
}
general_case:
ret = bf_ziv_rounding(r, T, prec, flags, bf_pow_generic, (void *)y);
}
}
done:
bf_delete(T);
r->sign = r_sign;
return ret;
}