Lines Matching refs:r
74 typedef int bf_op2_func_t(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec,
159 slimb_t r; in sat_add() local
160 r = a + b; in sat_add()
162 if (((a ^ r) & (b ^ r)) < 0) in sat_add()
163 r = (a >> (LIMB_BITS - 1)) ^ (((limb_t)1 << (LIMB_BITS - 1)) - 1); in sat_add()
164 return r; in sat_add()
184 void bf_init(bf_context_t *s, bf_t *r) in bf_init() argument
186 r->ctx = s; in bf_init()
187 r->sign = 0; in bf_init()
188 r->expn = BF_EXP_ZERO; in bf_init()
189 r->len = 0; in bf_init()
190 r->tab = NULL; in bf_init()
194 int bf_resize(bf_t *r, limb_t len) in bf_resize() argument
198 if (len != r->len) { in bf_resize()
199 tab = bf_realloc(r->ctx, r->tab, len * sizeof(limb_t)); in bf_resize()
202 r->tab = tab; in bf_resize()
203 r->len = len; in bf_resize()
209 int bf_set_ui(bf_t *r, uint64_t a) in bf_set_ui() argument
211 r->sign = 0; in bf_set_ui()
213 r->expn = BF_EXP_ZERO; in bf_set_ui()
214 bf_resize(r, 0); /* cannot fail */ in bf_set_ui()
223 if (bf_resize(r, 1)) in bf_set_ui()
226 r->tab[0] = a << shift; in bf_set_ui()
227 r->expn = LIMB_BITS - shift; in bf_set_ui()
233 if (bf_resize(r, 2)) in bf_set_ui()
238 r->tab[0] = a0 << shift; in bf_set_ui()
239 r->tab[1] = (a1 << shift) | (a0 >> (LIMB_BITS - shift)); in bf_set_ui()
240 r->expn = 2 * LIMB_BITS - shift; in bf_set_ui()
245 bf_set_nan(r); in bf_set_ui()
250 int bf_set_si(bf_t *r, int64_t a) in bf_set_si() argument
255 ret = bf_set_ui(r, -a); in bf_set_si()
256 r->sign = 1; in bf_set_si()
258 ret = bf_set_ui(r, a); in bf_set_si()
263 void bf_set_nan(bf_t *r) in bf_set_nan() argument
265 bf_resize(r, 0); /* cannot fail */ in bf_set_nan()
266 r->expn = BF_EXP_NAN; in bf_set_nan()
267 r->sign = 0; in bf_set_nan()
270 void bf_set_zero(bf_t *r, int is_neg) in bf_set_zero() argument
272 bf_resize(r, 0); /* cannot fail */ in bf_set_zero()
273 r->expn = BF_EXP_ZERO; in bf_set_zero()
274 r->sign = is_neg; in bf_set_zero()
277 void bf_set_inf(bf_t *r, int is_neg) in bf_set_inf() argument
279 bf_resize(r, 0); /* cannot fail */ in bf_set_inf()
280 r->expn = BF_EXP_INF; in bf_set_inf()
281 r->sign = is_neg; in bf_set_inf()
285 int bf_set(bf_t *r, const bf_t *a) in bf_set() argument
287 if (r == a) in bf_set()
289 if (bf_resize(r, a->len)) { in bf_set()
290 bf_set_nan(r); in bf_set()
293 r->sign = a->sign; in bf_set()
294 r->expn = a->expn; in bf_set()
295 memcpy(r->tab, a->tab, a->len * sizeof(limb_t)); in bf_set()
300 void bf_move(bf_t *r, bf_t *a) in bf_move() argument
302 bf_context_t *s = r->ctx; in bf_move()
303 if (r == a) in bf_move()
305 bf_free(s, r->tab); in bf_move()
306 *r = *a; in bf_move()
373 static inline limb_t scan_bit_nz(const bf_t *r, slimb_t bit_pos) in scan_bit_nz() argument
381 v = r->tab[pos] & limb_mask(0, bit_pos & (LIMB_BITS - 1)); in scan_bit_nz()
386 if (r->tab[pos] != 0) in scan_bit_nz()
395 static int bf_get_rnd_add(int *pret, const bf_t *r, limb_t l, in bf_get_rnd_add() argument
405 bit0 = scan_bit_nz(r, l * LIMB_BITS - 1 - bf_max(0, prec + 1)); in bf_get_rnd_add()
409 bit1 = get_bit(r->tab, l, l * LIMB_BITS - 1 - prec); in bf_get_rnd_add()
423 get_bit(r->tab, l, l * LIMB_BITS - 1 - (prec - 1)); in bf_get_rnd_add()
429 if (r->sign == (rnd_mode == BF_RNDD)) in bf_get_rnd_add()
448 static int bf_set_overflow(bf_t *r, int sign, limb_t prec, bf_flags_t flags) in bf_set_overflow() argument
460 bf_set_inf(r, sign); in bf_set_overflow()
464 if (bf_resize(r, l)) { in bf_set_overflow()
465 bf_set_nan(r); in bf_set_overflow()
468 r->tab[0] = limb_mask((-prec) & (LIMB_BITS - 1), in bf_set_overflow()
471 r->tab[i] = (limb_t)-1; in bf_set_overflow()
473 r->expn = e_max; in bf_set_overflow()
474 r->sign = sign; in bf_set_overflow()
484 static int __bf_round(bf_t *r, limb_t prec1, bf_flags_t flags, limb_t l, in __bf_round() argument
499 prec = r->expn + prec1; in __bf_round()
502 } else if (unlikely(r->expn < e_min) && (flags & BF_FLAG_SUBNORMAL)) { in __bf_round()
506 prec = prec1 - (e_min - r->expn); in __bf_round()
513 add_one = bf_get_rnd_add(&ret, r, l, prec, rnd_mode); in __bf_round()
517 bf_resize(r, 1); /* cannot fail */ in __bf_round()
518 r->tab[0] = (limb_t)1 << (LIMB_BITS - 1); in __bf_round()
519 r->expn += 1 - prec; in __bf_round()
534 v = r->tab[i] + carry; in __bf_round()
536 r->tab[i] = v; in __bf_round()
544 a = r->tab[i]; in __bf_round()
545 r->tab[i] = (a >> 1) | (v << (LIMB_BITS - 1)); in __bf_round()
548 r->expn++; in __bf_round()
553 if (unlikely(r->expn < e_min)) { in __bf_round()
561 bf_set_zero(r, r->sign); in __bf_round()
567 if (unlikely(r->expn > e_max)) in __bf_round()
568 return bf_set_overflow(r, r->sign, prec1, flags); in __bf_round()
576 r->tab[i] &= limb_mask(shift, LIMB_BITS - 1); in __bf_round()
581 while (r->tab[i] == 0) in __bf_round()
585 memmove(r->tab, r->tab + i, l * sizeof(limb_t)); in __bf_round()
587 bf_resize(r, l); /* cannot fail */ in __bf_round()
592 int bf_normalize_and_round(bf_t *r, limb_t prec1, bf_flags_t flags) in bf_normalize_and_round() argument
599 l = r->len; in bf_normalize_and_round()
600 while (l > 0 && r->tab[l - 1] == 0) in bf_normalize_and_round()
604 r->expn = BF_EXP_ZERO; in bf_normalize_and_round()
605 bf_resize(r, 0); /* cannot fail */ in bf_normalize_and_round()
608 r->expn -= (r->len - l) * LIMB_BITS; in bf_normalize_and_round()
610 v = r->tab[l - 1]; in bf_normalize_and_round()
615 a = r->tab[i]; in bf_normalize_and_round()
616 r->tab[i] = (a << shift) | (v >> (LIMB_BITS - shift)); in bf_normalize_and_round()
619 r->expn -= shift; in bf_normalize_and_round()
621 ret = __bf_round(r, prec1, flags, l, 0); in bf_normalize_and_round()
667 int bf_round(bf_t *r, limb_t prec, bf_flags_t flags) in bf_round() argument
669 if (r->len == 0) in bf_round()
671 return __bf_round(r, prec, flags, r->len, 0); in bf_round()
868 static int bf_add_internal(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_add_internal() argument
887 bf_set_zero(r, (flags & BF_RND_MASK) == BF_RNDD); in bf_add_internal()
894 bf_set_nan(r); in bf_add_internal()
897 bf_set_nan(r); in bf_add_internal()
900 bf_set_inf(r, a_sign); in bf_add_internal()
904 bf_set(r, a); in bf_add_internal()
905 r->sign = a_sign; in bf_add_internal()
912 r->sign = a_sign; in bf_add_internal()
913 r->expn = a->expn; in bf_add_internal()
930 if (bf_resize(r, r_len)) in bf_add_internal()
991 r->tab[i] = u; in bf_add_internal()
994 r->tab[0] |= (z != 0); in bf_add_internal()
998 if (bf_resize(r, r_len + 1)) in bf_add_internal()
1000 r->tab[r_len] = 1; in bf_add_internal()
1001 r->expn += LIMB_BITS; in bf_add_internal()
1004 ret = bf_normalize_and_round(r, prec, flags); in bf_add_internal()
1008 bf_set_nan(r); in bf_add_internal()
1012 static int __bf_add(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in __bf_add() argument
1015 return bf_add_internal(r, a, b, prec, flags, 0); in __bf_add()
1018 static int __bf_sub(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in __bf_sub() argument
1021 return bf_add_internal(r, a, b, prec, flags, 1); in __bf_sub()
1165 limb_t i, r; in mp_mul_basecase() local
1169 r = mp_add_mul1(result + i, op1, op1_size, op2[i]); in mp_mul_basecase()
1170 result[i + op1_size] = r; in mp_mul_basecase()
1182 bf_t r_s, *r = &r_s; in mp_mul() local
1183 r->tab = result; in mp_mul()
1185 if (fft_mul(s, r, (limb_t *)op1, op1_size, in mp_mul()
1227 limb_t n1m, n_adj, q, r, ah; in udiv1norm() local
1239 r = (limb_t)a + (ah & d); in udiv1norm()
1240 *pr = r; in udiv1norm()
1246 limb_t b, limb_t r) in mp_div1norm() argument
1254 tabr[i] = udiv1norm(&r, r, taba[i], b, b_inv); in mp_div1norm()
1259 a1 = ((dlimb_t)r << LIMB_BITS) | taba[i]; in mp_div1norm()
1261 r = a1 % b; in mp_div1norm()
1264 return r; in mp_div1norm()
1278 limb_t r, a, c, q, v, b1, b1_inv, n, dummy_r; in mp_divnorm() local
1319 r = al % b1; in mp_divnorm()
1321 r = mp_sub_mul1(taba + i, tabb, nb, q); in mp_divnorm()
1324 a = v - r; in mp_divnorm()
1524 int bf_mul(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_mul() argument
1538 bf_set_nan(r); in bf_mul()
1543 bf_set_nan(r); in bf_mul()
1546 bf_set_inf(r, r_sign); in bf_mul()
1550 bf_set_zero(r, r_sign); in bf_mul()
1573 if (r == a) in bf_mul()
1575 if (r == b) in bf_mul()
1577 if (fft_mul(r->ctx, r, a_tab, a_len, b_tab, b_len, mul_flags)) in bf_mul()
1582 if (r == a || r == b) { in bf_mul()
1583 bf_init(r->ctx, &tmp); in bf_mul()
1584 r1 = r; in bf_mul()
1585 r = &tmp; in bf_mul()
1587 if (bf_resize(r, a_len + b_len)) { in bf_mul()
1589 bf_set_nan(r); in bf_mul()
1593 mp_mul_basecase(r->tab, a_tab, a_len, b_tab, b_len); in bf_mul()
1595 r->sign = r_sign; in bf_mul()
1596 r->expn = a->expn + b->expn; in bf_mul()
1597 ret = bf_normalize_and_round(r, prec, flags); in bf_mul()
1599 if (r == &tmp) in bf_mul()
1606 int bf_mul_2exp(bf_t *r, slimb_t e, limb_t prec, bf_flags_t flags) in bf_mul_2exp() argument
1609 if (r->len == 0) in bf_mul_2exp()
1614 r->expn += e; in bf_mul_2exp()
1615 return __bf_round(r, prec, flags, r->len, 0); in bf_mul_2exp()
1638 static void bf_tdivremu(bf_t *q, bf_t *r, in bf_tdivremu() argument
1643 bf_set(r, a); in bf_tdivremu()
1647 bf_mul(r, q, b, BF_PREC_INF, BF_RNDZ); in bf_tdivremu()
1648 bf_sub(r, a, r, BF_PREC_INF, BF_RNDZ); in bf_tdivremu()
1652 static int __bf_div(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in __bf_div() argument
1655 bf_context_t *s = r->ctx; in __bf_div()
1662 bf_set_nan(r); in __bf_div()
1665 bf_set_nan(r); in __bf_div()
1668 bf_set_inf(r, r_sign); in __bf_div()
1671 bf_set_zero(r, r_sign); in __bf_div()
1676 bf_set_nan(r); in __bf_div()
1679 bf_set_zero(r, r_sign); in __bf_div()
1683 bf_set_inf(r, r_sign); in __bf_div()
1703 if (bf_resize(r, n + 1)) in __bf_div()
1705 if (mp_divnorm(s, r->tab, taba, na, b->tab, nb)) { in __bf_div()
1712 r->tab[0] |= 1; in __bf_div()
1713 bf_free(r->ctx, taba); in __bf_div()
1714 r->expn = a->expn - b->expn + LIMB_BITS; in __bf_div()
1715 r->sign = r_sign; in __bf_div()
1716 ret = bf_normalize_and_round(r, prec, flags); in __bf_div()
1720 bf_set_nan(r); in __bf_div()
1732 int bf_divrem(bf_t *q, bf_t *r, const bf_t *a, const bf_t *b, in bf_divrem() argument
1741 assert(r != a && r != b); in bf_divrem()
1742 assert(q != r); in bf_divrem()
1747 bf_set_nan(r); in bf_divrem()
1750 bf_set_nan(r); in bf_divrem()
1753 bf_set(r, a); in bf_divrem()
1754 return bf_round(r, prec, flags); in bf_divrem()
1792 bf_tdivremu(q, r, a1, b1); in bf_divrem()
1793 if (bf_is_nan(q) || bf_is_nan(r)) in bf_divrem()
1796 if (r->len != 0) { in bf_divrem()
1800 res = bf_cmpu(r, b1); in bf_divrem()
1811 ret |= bf_sub(r, r, b1, BF_PREC_INF, BF_RNDZ); in bf_divrem()
1817 r->sign ^= a->sign; in bf_divrem()
1819 return bf_round(r, prec, flags); in bf_divrem()
1822 bf_set_nan(r); in bf_divrem()
1826 int bf_rem(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_rem() argument
1832 bf_init(r->ctx, q); in bf_rem()
1833 ret = bf_divrem(q, r, a, b, prec, flags, rnd_mode); in bf_rem()
1847 int bf_remquo(slimb_t *pq, bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_remquo() argument
1853 bf_init(r->ctx, q); in bf_remquo()
1854 ret = bf_divrem(q, r, a, b, prec, flags, rnd_mode); in bf_remquo()
1868 static limb_t mp_mod1(const limb_t *tab, limb_t n, limb_t m, limb_t r) in mp_mod1() argument
1874 t = ((dlimb_t)r << LIMB_BITS) | tab[i]; in mp_mod1()
1875 r = t % m; in mp_mod1()
1877 return r; in mp_mod1()
1889 limb_t s1, r1, s, r, q, u, num; in mp_sqrtrem1() local
1904 r = (u << 8) | ((a >> (LIMB_BITS - 32)) & 0xff); in mp_sqrtrem1()
1905 r -= q * q; in mp_sqrtrem1()
1906 if ((slimb_t)r < 0) { in mp_sqrtrem1()
1908 r += 2 * s + 1; in mp_sqrtrem1()
1913 r1 = r; in mp_sqrtrem1()
1919 r = (u << 16) | ((a >> (LIMB_BITS - 64)) & 0xffff); in mp_sqrtrem1()
1920 r -= q * q; in mp_sqrtrem1()
1921 if ((slimb_t)r < 0) { in mp_sqrtrem1()
1923 r += 2 * s + 1; in mp_sqrtrem1()
1926 *pr = r; in mp_sqrtrem1()
1933 limb_t s, r; in bf_isqrt() local
1939 s = mp_sqrtrem1(&r, a << k); in bf_isqrt()
1947 dlimb_t r, num; in mp_sqrtrem2() local
1958 r = ((dlimb_t)u << l) | (a0 & (((limb_t)1 << l) - 1)); in mp_sqrtrem2()
1960 r -= (dlimb_t)1 << LIMB_BITS; /* special case when q=2^l */ in mp_sqrtrem2()
1962 r -= q * q; in mp_sqrtrem2()
1963 if ((slimb_t)(r >> LIMB_BITS) < 0) { in mp_sqrtrem2()
1965 r += 2 * (dlimb_t)s + 1; in mp_sqrtrem2()
1968 taba[0] = r; in mp_sqrtrem2()
1969 return r >> LIMB_BITS; in mp_sqrtrem2()
2073 int bf_sqrtrem(bf_t *r, bf_t *rem1, const bf_t *a) in bf_sqrtrem() argument
2079 bf_set_nan(r); in bf_sqrtrem()
2083 bf_set(r, a); in bf_sqrtrem()
2090 bf_set_nan(r); in bf_sqrtrem()
2097 bf_sqrt(r, a, (a->expn + 1) / 2, BF_RNDZ); in bf_sqrtrem()
2098 bf_rint(r, BF_RNDZ); in bf_sqrtrem()
2104 bf_init(r->ctx, rem); in bf_sqrtrem()
2107 bf_mul(rem, r, r, BF_PREC_INF, BF_RNDZ); in bf_sqrtrem()
2126 int bf_sqrt(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_sqrt() argument
2131 assert(r != a); in bf_sqrt()
2135 bf_set_nan(r); in bf_sqrt()
2139 bf_set(r, a); in bf_sqrt()
2144 bf_set_nan(r); in bf_sqrt()
2154 if (bf_resize(r, n)) in bf_sqrt()
2167 if (mp_sqrtrem(s, r->tab, a1, n)) { in bf_sqrt()
2179 r->tab[0] |= 1; in bf_sqrt()
2180 r->sign = 0; in bf_sqrt()
2181 r->expn = (a->expn + 1) >> 1; in bf_sqrt()
2182 ret = bf_round(r, prec, flags); in bf_sqrt()
2186 bf_set_nan(r); in bf_sqrt()
2190 static no_inline int bf_op2(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_op2() argument
2196 if (r == a || r == b) { in bf_op2()
2197 bf_init(r->ctx, &tmp); in bf_op2()
2199 bf_move(r, &tmp); in bf_op2()
2201 ret = func(r, a, b, prec, flags); in bf_op2()
2206 int bf_add(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_add() argument
2209 return bf_op2(r, a, b, prec, flags, __bf_add); in bf_add()
2212 int bf_sub(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_sub() argument
2215 return bf_op2(r, a, b, prec, flags, __bf_sub); in bf_sub()
2218 int bf_div(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, in bf_div() argument
2221 return bf_op2(r, a, b, prec, flags, __bf_div); in bf_div()
2224 int bf_mul_ui(bf_t *r, const bf_t *a, uint64_t b1, limb_t prec, in bf_mul_ui() argument
2229 bf_init(r->ctx, &b); in bf_mul_ui()
2231 ret |= bf_mul(r, a, &b, prec, flags); in bf_mul_ui()
2236 int bf_mul_si(bf_t *r, const bf_t *a, int64_t b1, limb_t prec, in bf_mul_si() argument
2241 bf_init(r->ctx, &b); in bf_mul_si()
2243 ret |= bf_mul(r, a, &b, prec, flags); in bf_mul_si()
2248 int bf_add_si(bf_t *r, const bf_t *a, int64_t b1, limb_t prec, in bf_add_si() argument
2254 bf_init(r->ctx, &b); in bf_add_si()
2256 ret |= bf_add(r, a, &b, prec, flags); in bf_add_si()
2261 static int bf_pow_ui(bf_t *r, const bf_t *a, limb_t b, limb_t prec, in bf_pow_ui() argument
2266 assert(r != a); in bf_pow_ui()
2268 return bf_set_ui(r, 1); in bf_pow_ui()
2269 ret = bf_set(r, a); in bf_pow_ui()
2272 ret |= bf_mul(r, r, r, prec, flags); in bf_pow_ui()
2274 ret |= bf_mul(r, r, a, prec, flags); in bf_pow_ui()
2279 static int bf_pow_ui_ui(bf_t *r, limb_t a1, limb_t b, in bf_pow_ui_ui() argument
2288 ret = bf_set_ui(r, mp_pow_dec[b]); in bf_pow_ui_ui()
2290 bf_init(r->ctx, &a); in bf_pow_ui_ui()
2292 ret |= bf_pow_ui(r, &a, b, prec, flags); in bf_pow_ui_ui()
2299 int bf_rint(bf_t *r, int rnd_mode) in bf_rint() argument
2301 return bf_round(r, 0, rnd_mode | BF_FLAG_RADPNT_PREC); in bf_rint()
2322 static int bf_logic_op(bf_t *r, const bf_t *a1, const bf_t *b1, int op) in bf_logic_op() argument
2330 assert(r != a1 && r != b1); in bf_logic_op()
2344 bf_init(r->ctx, a); in bf_logic_op()
2355 bf_init(r->ctx, b); in bf_logic_op()
2376 if (bf_resize(r, l)) in bf_logic_op()
2386 r->tab[i] = bf_logic_op1(v1, v2, op) ^ r_mask; in bf_logic_op()
2388 r->expn = l * LIMB_BITS; in bf_logic_op()
2389 r->sign = r_sign; in bf_logic_op()
2390 bf_normalize_and_round(r, BF_PREC_INF, BF_RNDZ); /* cannot fail */ in bf_logic_op()
2392 if (bf_add_si(r, r, -1, BF_PREC_INF, BF_RNDZ)) in bf_logic_op()
2403 bf_set_nan(r); in bf_logic_op()
2409 int bf_logic_or(bf_t *r, const bf_t *a, const bf_t *b) in bf_logic_or() argument
2411 return bf_logic_op(r, a, b, BF_LOGIC_OR); in bf_logic_or()
2415 int bf_logic_xor(bf_t *r, const bf_t *a, const bf_t *b) in bf_logic_xor() argument
2417 return bf_logic_op(r, a, b, BF_LOGIC_XOR); in bf_logic_xor()
2421 int bf_logic_and(bf_t *r, const bf_t *a, const bf_t *b) in bf_logic_and() argument
2423 return bf_logic_op(r, a, b, BF_LOGIC_AND); in bf_logic_and()
2693 static int bf_integer_from_radix_rec(bf_t *r, const limb_t *tab, in bf_integer_from_radix_rec() argument
2699 ret = bf_set_ui(r, tab[0]); in bf_integer_from_radix_rec()
2713 ret = bf_integer_from_radix_rec(r, tab + n2, n1, level + 1, n0, in bf_integer_from_radix_rec()
2717 ret = bf_mul(r, r, B, BF_PREC_INF, BF_RNDZ); in bf_integer_from_radix_rec()
2720 bf_init(r->ctx, T); in bf_integer_from_radix_rec()
2724 ret = bf_add(r, r, T, BF_PREC_INF, BF_RNDZ); in bf_integer_from_radix_rec()
2732 static int bf_integer_from_radix(bf_t *r, const limb_t *tab, in bf_integer_from_radix() argument
2735 bf_context_t *s = r->ctx; in bf_integer_from_radix()
2746 bf_init(r->ctx, &pow_tab[i]); in bf_integer_from_radix()
2747 ret = bf_integer_from_radix_rec(r, tab, n, 0, n, radixl, pow_tab); in bf_integer_from_radix()
2756 int bf_mul_pow_radix(bf_t *r, const bf_t *T, limb_t radix, in bf_mul_pow_radix() argument
2764 return bf_set(r, T); in bf_mul_pow_radix()
2766 ret = bf_set(r, T); in bf_mul_pow_radix()
2767 ret |= bf_round(r, prec, flags); in bf_mul_pow_radix()
2777 bf_init(r->ctx, B); in bf_mul_pow_radix()
2782 ret |= bf_div(r, T, B, T->len * LIMB_BITS, BF_RNDN); in bf_mul_pow_radix()
2784 ret |= bf_mul(r, T, B, BF_PREC_INF, BF_RNDN); in bf_mul_pow_radix()
2798 ret |= bf_div(r, T, B, prec1 + extra_bits, BF_RNDN | BF_FLAG_EXT_EXP); in bf_mul_pow_radix()
2800 ret |= bf_mul(r, T, B, prec1 + extra_bits, BF_RNDN | BF_FLAG_EXT_EXP); in bf_mul_pow_radix()
2804 !bf_can_round(r, prec, flags & BF_RND_MASK, prec1) && in bf_mul_pow_radix()
2811 ret = bf_round(r, prec, flags) | (ret & BF_ST_INEXACT); in bf_mul_pow_radix()
2878 static int bf_atof_internal(bf_t *r, slimb_t *pexponent, in bf_atof_internal() argument
2893 bf_set_nan(r); in bf_atof_internal()
2928 bf_set_nan(r); in bf_atof_internal()
2936 bf_set_inf(r, is_neg); in bf_atof_internal()
2947 a = r; in bf_atof_internal()
2951 bf_init(r->ctx, a); in bf_atof_internal()
2954 a = r; in bf_atof_internal()
3028 bf_set_nan(r); in bf_atof_internal()
3041 bf_set_nan(r); in bf_atof_internal()
3069 bf_set_zero(r, is_neg); in bf_atof_internal()
3072 bf_set_inf(r, is_neg); in bf_atof_internal()
3099 bf_set_zero(r, is_neg); in bf_atof_internal()
3105 bf_init(r->ctx, T); in bf_atof_internal()
3107 bf_set_nan(r); in bf_atof_internal()
3114 ret = bf_set(r, T); in bf_atof_internal()
3116 ret = bf_mul_pow_radix(r, T, radix, expn, prec, flags); in bf_atof_internal()
3137 int bf_atof2(bf_t *r, slimb_t *pexponent, in bf_atof2() argument
3141 return bf_atof_internal(r, pexponent, str, pnext, radix, prec, flags, in bf_atof2()
3145 int bf_atof(bf_t *r, const char *str, const char **pnext, int radix, in bf_atof() argument
3149 return bf_atof_internal(r, &dummy_exp, str, pnext, radix, prec, flags, FALSE); in bf_atof()
3459 static int bf_integer_to_radix(bf_t *r, const bf_t *a, limb_t radixl) in bf_integer_to_radix() argument
3461 bf_context_t *s = r->ctx; in bf_integer_to_radix()
3466 r_len = r->len; in bf_integer_to_radix()
3472 bf_init(r->ctx, &pow_tab[i]); in bf_integer_to_radix()
3474 ret = bf_integer_to_radix_rec(pow_tab, r->tab, a, r_len, 0, r_len, radixl, in bf_integer_to_radix()
3487 static int bf_convert_to_radix(bf_t *r, slimb_t *pE, in bf_convert_to_radix() argument
3499 return bf_set(r, a); in bf_convert_to_radix()
3525 ret = bf_pow_ui_ui(r, radix, e, prec + extra_bits, in bf_convert_to_radix()
3528 ret |= bf_mul(r, r, a, prec + extra_bits, in bf_convert_to_radix()
3531 ret |= bf_div(r, a, r, prec + extra_bits, in bf_convert_to_radix()
3538 !bf_can_round(r, r->expn, rnd_mode, prec)) { in bf_convert_to_radix()
3543 ret = bf_rint(r, rnd_mode); in bf_convert_to_radix()
3553 bf_init(r->ctx, B); in bf_convert_to_radix()
3559 res = bf_cmpu(r, B); in bf_convert_to_radix()
4218 typedef int ZivFunc(bf_t *r, const bf_t *a, limb_t prec, void *opaque);
4220 static int bf_ziv_rounding(bf_t *r, const bf_t *a, in bf_ziv_rounding() argument
4230 f(r, a, prec, opaque); in bf_ziv_rounding()
4236 ret = f(r, a, prec1, opaque); in bf_ziv_rounding()
4248 if (bf_can_round(r, prec, rnd_mode, prec1)) { in bf_ziv_rounding()
4256 if (r->len == 0) in bf_ziv_rounding()
4259 return __bf_round(r, prec, flags, r->len, ret); in bf_ziv_rounding()
4263 static int bf_add_epsilon(bf_t *r, const bf_t *a, slimb_t e, int e_sign, in bf_add_epsilon() argument
4273 ret = bf_add(r, r, T, prec, flags); in bf_add_epsilon()
4280 static int bf_exp_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_exp_internal() argument
4282 bf_context_t *s = r->ctx; in bf_exp_internal()
4286 assert(r != a); in bf_exp_internal()
4327 bf_set_ui(r, 1); in bf_exp_internal()
4331 bf_mul(r, r, U, prec1, BF_RNDN); in bf_exp_internal()
4332 bf_add_si(r, r, 1, prec1, BF_RNDN); in bf_exp_internal()
4340 bf_mul(r, r, r, prec1, BF_RNDN | BF_FLAG_EXT_EXP); in bf_exp_internal()
4344 bf_mul_2exp(r, n, BF_PREC_INF, BF_RNDZ | BF_FLAG_EXT_EXP); in bf_exp_internal()
4350 static int check_exp_underflow_overflow(bf_context_t *s, bf_t *r, in check_exp_underflow_overflow() argument
4375 return bf_set_overflow(r, 0, prec, flags); in check_exp_underflow_overflow()
4388 bf_set_ui(r, 1); in check_exp_underflow_overflow()
4389 r->expn = e_min; in check_exp_underflow_overflow()
4391 bf_set_zero(r, 0); in check_exp_underflow_overflow()
4400 int bf_exp(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_exp() argument
4402 bf_context_t *s = r->ctx; in bf_exp()
4404 assert(r != a); in bf_exp()
4407 bf_set_nan(r); in bf_exp()
4410 bf_set_zero(r, 0); in bf_exp()
4412 bf_set_inf(r, 0); in bf_exp()
4414 bf_set_ui(r, 1); in bf_exp()
4419 ret = check_exp_underflow_overflow(s, r, a, a, prec, flags); in bf_exp()
4424 bf_set_ui(r, 1); in bf_exp()
4425 return bf_add_epsilon(r, r, -(prec + 2), a->sign, prec, flags); in bf_exp()
4428 return bf_ziv_rounding(r, a, prec, flags, bf_exp_internal, NULL); in bf_exp()
4431 static int bf_log_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_log_internal() argument
4433 bf_context_t *s = r->ctx; in bf_log_internal()
4439 assert(r != a); in bf_log_internal()
4501 bf_set_ui(r, 0); in bf_log_internal()
4506 bf_add(r, r, U, prec1, BF_RNDN); in bf_log_internal()
4507 bf_mul(r, r, Y2, prec1, BF_RNDN); in bf_log_internal()
4509 bf_add_si(r, r, 1, prec1, BF_RNDN); in bf_log_internal()
4510 bf_mul(r, r, Y, prec1, BF_RNDN); in bf_log_internal()
4519 bf_mul_2exp(r, K + 1, BF_PREC_INF, BF_RNDZ); in bf_log_internal()
4524 bf_add(r, r, T, prec1, BF_RNDN); in bf_log_internal()
4530 int bf_log(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_log() argument
4532 bf_context_t *s = r->ctx; in bf_log()
4535 assert(r != a); in bf_log()
4538 bf_set_nan(r); in bf_log()
4542 bf_set_nan(r); in bf_log()
4545 bf_set_inf(r, 0); in bf_log()
4549 bf_set_inf(r, 1); in bf_log()
4554 bf_set_nan(r); in bf_log()
4560 bf_set_zero(r, 0); in bf_log()
4566 return bf_ziv_rounding(r, a, prec, flags, bf_log_internal, NULL); in bf_log()
4570 static int bf_pow_generic(bf_t *r, const bf_t *x, limb_t prec, void *opaque) in bf_pow_generic() argument
4572 bf_context_t *s = r->ctx; in bf_pow_generic()
4583 bf_set_nan(r); in bf_pow_generic()
4585 bf_exp_internal(r, T, prec1, NULL); /* no overflow/underlow test needed */ in bf_pow_generic()
4591 static int bf_pow_int(bf_t *r, const bf_t *x, limb_t prec, void *opaque) in bf_pow_int() argument
4593 bf_context_t *s = r->ctx; in bf_pow_int()
4605 ret = bf_pow_ui(r, x, y1 < 0 ? -y1 : y1, prec1, BF_RNDN | BF_FLAG_EXT_EXP); in bf_pow_int()
4609 ret |= bf_div(r, T, r, prec1, BF_RNDN | BF_FLAG_EXT_EXP); in bf_pow_int()
4618 static BOOL check_exact_power2n(bf_t *r, const bf_t *x, slimb_t n) in check_exact_power2n() argument
4620 bf_context_t *s = r->ctx; in check_exact_power2n()
4647 bf_set(T, r); in check_exact_power2n()
4648 if (bf_sqrtrem(r, NULL, T) != 0) in check_exact_power2n()
4651 r->expn += er; in check_exact_power2n()
4656 int bf_pow(bf_t *r, const bf_t *x, const bf_t *y, limb_t prec, bf_flags_t flags) in bf_pow() argument
4658 bf_context_t *s = r->ctx; in bf_pow()
4668 bf_set_ui(r, 1); in bf_pow()
4670 bf_set_nan(r); in bf_pow()
4673 bf_set_ui(r, 1); in bf_pow()
4674 cmp_x_abs_1 = bf_cmpu(x, r); in bf_pow()
4677 bf_set_nan(r); in bf_pow()
4683 bf_set_nan(r); in bf_pow()
4686 bf_set_zero(r, 0); in bf_pow()
4688 bf_set_inf(r, 0); in bf_pow()
4694 bf_set_inf(r, y_is_odd & x->sign); in bf_pow()
4700 bf_set_zero(r, y_is_odd & x->sign); in bf_pow()
4713 bf_set_nan(r); in bf_pow()
4728 bf_set_ui(r, 1); in bf_pow()
4729 if (bf_cmp_eq(T, r)) { in bf_pow()
4748 ret = check_exp_underflow_overflow(s, r, al, ah, prec, flags); in bf_pow()
4763 bf_set_ui(r, 1); in bf_pow()
4764 ret = bf_mul_2exp(r, e, prec, flags); in bf_pow()
4773 return bf_set_overflow(r, 0, BF_PREC_INF, flags); in bf_pow()
4775 ret = bf_pow_ui(r, T, y1, BF_PREC_INF, BF_RNDZ); in bf_pow()
4789 bf_mul_si(r, y, T_bits - 1, LIMB_BITS, BF_RNDZ); in bf_pow()
4790 bf_get_limb(&e, r, 0); in bf_pow()
4794 ret = bf_ziv_rounding(r, T, prec, flags, bf_pow_int, (void *)y); in bf_pow()
4799 if (y_emin < 0 && check_exact_power2n(r, T, -y_emin)) { in bf_pow()
4804 bf_print_str("r", r); in bf_pow()
4806 bf_set(T, r); in bf_pow()
4817 ret = bf_ziv_rounding(r, T, prec, flags, bf_pow_generic, (void *)y); in bf_pow()
4822 r->sign = r_sign; in bf_pow()
4827 static void bf_sqrt_sin(bf_t *r, const bf_t *x, limb_t prec1) in bf_sqrt_sin() argument
4829 bf_context_t *s = r->ctx; in bf_sqrt_sin()
4833 bf_mul(r, T, T, prec1, BF_RNDN); in bf_sqrt_sin()
4835 bf_add(T, T, r, prec1, BF_RNDN); in bf_sqrt_sin()
4837 bf_sqrt(r, T, prec1, BF_RNDF); in bf_sqrt_sin()
4846 bf_t r_s, *r = &r_s; in bf_sincos() local
4854 bf_init(s1, r); in bf_sincos()
4893 bf_set_ui(r, 1); in bf_sincos()
4898 bf_mul(r, r, U, prec1, BF_RNDN); in bf_sincos()
4899 bf_neg(r); in bf_sincos()
4901 bf_add_si(r, r, 1, prec1, BF_RNDN); in bf_sincos()
4909 bf_mul(T, r, r, prec1, BF_RNDN); in bf_sincos()
4910 bf_mul_2exp(r, 1, BF_PREC_INF, BF_RNDZ); in bf_sincos()
4911 bf_add(r, r, T, prec1, BF_RNDN); in bf_sincos()
4912 bf_mul_2exp(r, 1, BF_PREC_INF, BF_RNDZ); in bf_sincos()
4918 bf_add_si(c, r, 1, prec1, BF_RNDN); in bf_sincos()
4920 bf_sqrt_sin(c, r, prec1); in bf_sincos()
4927 bf_sqrt_sin(s, r, prec1); in bf_sincos()
4930 bf_add_si(s, r, 1, prec1, BF_RNDN); in bf_sincos()
4934 bf_delete(r); in bf_sincos()
4938 static int bf_cos_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_cos_internal() argument
4940 return bf_sincos(NULL, r, a, prec); in bf_cos_internal()
4943 int bf_cos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_cos() argument
4947 bf_set_nan(r); in bf_cos()
4950 bf_set_nan(r); in bf_cos()
4953 bf_set_ui(r, 1); in bf_cos()
4964 bf_set_ui(r, 1); in bf_cos()
4965 return bf_add_epsilon(r, r, e, 1, prec, flags); in bf_cos()
4969 return bf_ziv_rounding(r, a, prec, flags, bf_cos_internal, NULL); in bf_cos()
4972 static int bf_sin_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_sin_internal() argument
4974 return bf_sincos(r, NULL, a, prec); in bf_sin_internal()
4977 int bf_sin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_sin() argument
4981 bf_set_nan(r); in bf_sin()
4984 bf_set_nan(r); in bf_sin()
4987 bf_set_zero(r, a->sign); in bf_sin()
4998 bf_set(r, a); in bf_sin()
4999 return bf_add_epsilon(r, r, e, 1 - a->sign, prec, flags); in bf_sin()
5003 return bf_ziv_rounding(r, a, prec, flags, bf_sin_internal, NULL); in bf_sin()
5006 static int bf_tan_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_tan_internal() argument
5008 bf_context_t *s = r->ctx; in bf_tan_internal()
5015 bf_sincos(r, T, a, prec1); in bf_tan_internal()
5016 bf_div(r, r, T, prec1, BF_RNDF); in bf_tan_internal()
5021 int bf_tan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_tan() argument
5023 assert(r != a); in bf_tan()
5026 bf_set_nan(r); in bf_tan()
5029 bf_set_nan(r); in bf_tan()
5032 bf_set_zero(r, a->sign); in bf_tan()
5043 bf_set(r, a); in bf_tan()
5044 return bf_add_epsilon(r, r, e, a->sign, prec, flags); in bf_tan()
5048 return bf_ziv_rounding(r, a, prec, flags, bf_tan_internal, NULL); in bf_tan()
5053 static int bf_atan_internal(bf_t *r, const bf_t *a, limb_t prec, in bf_atan_internal() argument
5056 bf_context_t *s = r->ctx; in bf_atan_internal()
5100 bf_set_ui(r, 0); in bf_atan_internal()
5105 bf_neg(r); in bf_atan_internal()
5106 bf_add(r, r, U, prec1, BF_RNDN); in bf_atan_internal()
5107 bf_mul(r, r, X2, prec1, BF_RNDN); in bf_atan_internal()
5109 bf_neg(r); in bf_atan_internal()
5110 bf_add_si(r, r, 1, prec1, BF_RNDN); in bf_atan_internal()
5111 bf_mul(r, r, T, prec1, BF_RNDN); in bf_atan_internal()
5114 bf_mul_2exp(r, K, BF_PREC_INF, BF_RNDZ); in bf_atan_internal()
5123 bf_neg(r); in bf_atan_internal()
5132 bf_add(r, T, r, prec1, BF_RNDN); in bf_atan_internal()
5139 int bf_atan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_atan() argument
5141 bf_context_t *s = r->ctx; in bf_atan()
5147 bf_set_nan(r); in bf_atan()
5151 bf_const_pi_signed(r, a->sign, prec, flags); in bf_atan()
5152 bf_mul_2exp(r, -1, BF_PREC_INF, BF_RNDZ); in bf_atan()
5155 bf_set_zero(r, a->sign); in bf_atan()
5166 bf_const_pi_signed(r, a->sign, prec, flags); in bf_atan()
5167 bf_mul_2exp(r, -2, BF_PREC_INF, BF_RNDZ); in bf_atan()
5177 bf_set(r, a); in bf_atan()
5178 return bf_add_epsilon(r, r, e, 1 - a->sign, prec, flags); in bf_atan()
5182 return bf_ziv_rounding(r, a, prec, flags, bf_atan_internal, (void *)FALSE); in bf_atan()
5185 static int bf_atan2_internal(bf_t *r, const bf_t *y, limb_t prec, void *opaque) in bf_atan2_internal() argument
5187 bf_context_t *s = r->ctx; in bf_atan2_internal()
5194 bf_set_nan(r); in bf_atan2_internal()
5209 ret = bf_atan(r, T, prec1, BF_RNDF); in bf_atan2_internal()
5215 bf_add(r, r, T, prec1, BF_RNDN); in bf_atan2_internal()
5223 int bf_atan2(bf_t *r, const bf_t *y, const bf_t *x, in bf_atan2() argument
5226 return bf_ziv_rounding(r, y, prec, flags, bf_atan2_internal, (void *)x); in bf_atan2()
5229 static int bf_asin_internal(bf_t *r, const bf_t *a, limb_t prec, void *opaque) in bf_asin_internal() argument
5231 bf_context_t *s = r->ctx; in bf_asin_internal()
5251 bf_sqrt(r, T, prec1, BF_RNDN); in bf_asin_internal()
5252 bf_div(T, a, r, prec1, BF_RNDN); in bf_asin_internal()
5255 bf_atan_internal(r, T, prec1, (void *)(intptr_t)is_acos); in bf_asin_internal()
5260 int bf_asin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_asin() argument
5262 bf_context_t *s = r->ctx; in bf_asin()
5268 bf_set_nan(r); in bf_asin()
5271 bf_set_nan(r); in bf_asin()
5274 bf_set_zero(r, a->sign); in bf_asin()
5283 bf_set_nan(r); in bf_asin()
5293 bf_set(r, a); in bf_asin()
5294 return bf_add_epsilon(r, r, e, a->sign, prec, flags); in bf_asin()
5298 return bf_ziv_rounding(r, a, prec, flags, bf_asin_internal, (void *)FALSE); in bf_asin()
5301 int bf_acos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags) in bf_acos() argument
5303 bf_context_t *s = r->ctx; in bf_acos()
5309 bf_set_nan(r); in bf_acos()
5312 bf_set_nan(r); in bf_acos()
5315 bf_const_pi(r, prec, flags); in bf_acos()
5316 bf_mul_2exp(r, -1, BF_PREC_INF, BF_RNDZ); in bf_acos()
5325 bf_set_nan(r); in bf_acos()
5328 bf_set_zero(r, 0); in bf_acos()
5332 return bf_ziv_rounding(r, a, prec, flags, bf_asin_internal, (void *)TRUE); in bf_acos()
5365 #define divdq(q, r, a1, a0, b) \ argument
5371 r = __t % __b; \
5384 #define divdq(q, r, a1, a0, b) \ argument
5390 r = __t % __b; \
5414 #define divdq_base(q, r, a1, a0)\ argument
5430 r = __a0;\
5437 #define divdq_base(q, r, a1, a0)\ argument
5444 r = a0 - q * __b;\
5445 __t1 = (r >= __b);\
5448 r -= __b;\
5467 limb_t q, r, m1; in fast_udiv_init() local
5472 divdq(q, r, ((limb_t)1 << l) - d, 0, d); in fast_udiv_init()
5473 (void)r; in fast_udiv_init()
5563 #define fast_shr_rem_dec(q, r, a, shift) q = fast_shr_dec(a, shift), r = a - q * mp_pow_dec[shift] argument
5649 limb_t t0, t1, r; in mp_mul1_dec() local
5654 divdq_base(l, r, t1, t0); in mp_mul1_dec()
5655 tabr[i] = r; in mp_mul1_dec()
5666 limb_t l, t0, t1, r; in mp_add_mul1_dec() local
5673 divdq_base(l, r, t1, t0); in mp_add_mul1_dec()
5674 tabr[i] = r; in mp_add_mul1_dec()
5686 limb_t l, t0, t1, r, a, v, c; in mp_sub_mul1_dec() local
5693 divdq_base(l, r, t1, t0); in mp_sub_mul1_dec()
5695 a = v - r; in mp_sub_mul1_dec()
5712 limb_t r; in mp_mul_basecase_dec() local
5717 r = mp_add_mul1_dec(result + i, op1, op1_size, op2[i]); in mp_mul_basecase_dec()
5718 result[i + op1_size] = r; in mp_mul_basecase_dec()
5725 limb_t b, limb_t r) in mp_div1_dec() argument
5737 if (r) in mp_div1_dec()
5738 r = base_div2; in mp_div1_dec()
5741 tabr[i] = (t0 >> 1) + r; in mp_div1_dec()
5742 r = 0; in mp_div1_dec()
5744 r = base_div2; in mp_div1_dec()
5746 if (r) in mp_div1_dec()
5747 r = 1; in mp_div1_dec()
5757 muldq(t1, t0, r, base); in mp_div1_dec()
5759 q = udiv1norm(&r, t1, t0, b, b_inv); in mp_div1_dec()
5767 muldq(t1, t0, r, base); in mp_div1_dec()
5771 q = udiv1norm(&r, t1, t0, b, b_inv); in mp_div1_dec()
5772 r >>= shift; in mp_div1_dec()
5778 muldq(t1, t0, r, base); in mp_div1_dec()
5780 divdq(q, r, t1, t0, b); in mp_div1_dec()
5784 return r; in mp_div1_dec()
5836 limb_t r, mult, t0, t1, a, c, q, v, *tabb; in mp_div_dec() local
5846 r = tabb1[nb - 1]; in mp_div_dec()
5847 assert(r != 0); in mp_div_dec()
5849 if (r >= BF_DEC_BASE / 2) { in mp_div_dec()
5866 mult = base / (r + 1); in mp_div_dec()
5891 divdq(q, r, t1, t0, tabb[nb - 1]); in mp_div_dec()
5895 r = mp_sub_mul1_dec(taba + i, tabb, nb, q); in mp_div_dec()
5900 a = v - r; in mp_div_dec()
5941 limb_t l, a, q, r; in mp_shr_dec() local
5947 fast_shr_rem_dec(q, r, a, shift); in mp_shr_dec()
5949 l = r; in mp_shr_dec()
5959 limb_t l, a, q, r; in mp_shl_dec() local
5965 fast_shr_rem_dec(q, r, a, LIMB_DIGITS - shift); in mp_shl_dec()
5966 tab_r[i] = r * mp_pow_dec[shift] + l; in mp_shl_dec()
5975 dlimb_t a, b, r; in mp_sqrtrem2_dec() local
5987 r = a - (dlimb_t)s * (dlimb_t)s; in mp_sqrtrem2_dec()
5988 divdq_base(r1, r0, r >> LIMB_BITS, r); in mp_sqrtrem2_dec()
6295 static inline limb_t scan_digit_nz(const bfdec_t *r, slimb_t bit_pos) in scan_digit_nz() argument
6305 fast_shr_rem_dec(q, v, r->tab[pos], shift + 1); in scan_digit_nz()
6311 if (r->tab[pos] != 0) in scan_digit_nz()
6358 static int bfdec_get_rnd_add(int *pret, const bfdec_t *r, limb_t l, in bfdec_get_rnd_add() argument
6369 digit0 = scan_digit_nz(r, l * LIMB_DIGITS - 1 - bf_max(0, prec + 1)); in bfdec_get_rnd_add()
6373 digit1 = get_digit(r->tab, l, l * LIMB_DIGITS - 1 - prec); in bfdec_get_rnd_add()
6387 get_digit(r->tab, l, l * LIMB_DIGITS - 1 - (prec - 1)) & 1; in bfdec_get_rnd_add()
6395 if (r->sign == (rnd_mode == BF_RNDD)) in bfdec_get_rnd_add()
6419 static int __bfdec_round(bfdec_t *r, limb_t prec1, bf_flags_t flags, limb_t l) in __bfdec_round() argument
6432 prec = r->expn + prec1; in __bfdec_round()
6435 } else if (unlikely(r->expn < e_min) && (flags & BF_FLAG_SUBNORMAL)) { in __bfdec_round()
6439 prec = prec1 - (e_min - r->expn); in __bfdec_round()
6447 add_one = bfdec_get_rnd_add(&ret, r, l, prec, rnd_mode); in __bfdec_round()
6451 bfdec_resize(r, 1); /* cannot fail because r is non zero */ in __bfdec_round()
6452 r->tab[0] = BF_DEC_BASE / 10; in __bfdec_round()
6453 r->expn += 1 - prec; in __bfdec_round()
6466 carry = mp_add_ui_dec(r->tab + pos, carry, l - pos); in __bfdec_round()
6469 mp_shr_dec(r->tab + pos, r->tab + pos, l - pos, 1, 1); in __bfdec_round()
6470 r->expn++; in __bfdec_round()
6475 if (unlikely(r->expn < e_min)) { in __bfdec_round()
6482 bfdec_set_zero(r, r->sign); in __bfdec_round()
6489 if (unlikely(r->expn > e_max)) { in __bfdec_round()
6490 bfdec_set_inf(r, r->sign); in __bfdec_round()
6501 r->tab[i] = fast_shr_dec(r->tab[i], shift) * in __bfdec_round()
6508 while (r->tab[i] == 0) in __bfdec_round()
6512 memmove(r->tab, r->tab + i, l * sizeof(limb_t)); in __bfdec_round()
6514 bfdec_resize(r, l); /* cannot fail */ in __bfdec_round()
6519 int bfdec_round(bfdec_t *r, limb_t prec, bf_flags_t flags) in bfdec_round() argument
6521 if (r->len == 0) in bfdec_round()
6523 return __bfdec_round(r, prec, flags, r->len); in bfdec_round()
6527 int bfdec_normalize_and_round(bfdec_t *r, limb_t prec1, bf_flags_t flags) in bfdec_normalize_and_round() argument
6533 l = r->len; in bfdec_normalize_and_round()
6534 while (l > 0 && r->tab[l - 1] == 0) in bfdec_normalize_and_round()
6538 r->expn = BF_EXP_ZERO; in bfdec_normalize_and_round()
6539 bfdec_resize(r, 0); /* cannot fail */ in bfdec_normalize_and_round()
6542 r->expn -= (r->len - l) * LIMB_DIGITS; in bfdec_normalize_and_round()
6544 v = r->tab[l - 1]; in bfdec_normalize_and_round()
6547 mp_shl_dec(r->tab, r->tab, l, shift, 0); in bfdec_normalize_and_round()
6548 r->expn -= shift; in bfdec_normalize_and_round()
6550 ret = __bfdec_round(r, prec1, flags, l); in bfdec_normalize_and_round()
6556 int bfdec_set_ui(bfdec_t *r, uint64_t v) in bfdec_set_ui() argument
6560 if (bfdec_resize(r, 3)) in bfdec_set_ui()
6562 r->tab[0] = v % BF_DEC_BASE; in bfdec_set_ui()
6564 r->tab[1] = v % BF_DEC_BASE; in bfdec_set_ui()
6565 r->tab[2] = v / BF_DEC_BASE; in bfdec_set_ui()
6566 r->expn = 3 * LIMB_DIGITS; in bfdec_set_ui()
6570 if (bfdec_resize(r, 2)) in bfdec_set_ui()
6572 r->tab[0] = v % BF_DEC_BASE; in bfdec_set_ui()
6573 r->tab[1] = v / BF_DEC_BASE; in bfdec_set_ui()
6574 r->expn = 2 * LIMB_DIGITS; in bfdec_set_ui()
6576 if (bfdec_resize(r, 1)) in bfdec_set_ui()
6578 r->tab[0] = v; in bfdec_set_ui()
6579 r->expn = LIMB_DIGITS; in bfdec_set_ui()
6581 r->sign = 0; in bfdec_set_ui()
6582 return bfdec_normalize_and_round(r, BF_PREC_INF, 0); in bfdec_set_ui()
6584 bfdec_set_nan(r); in bfdec_set_ui()
6588 int bfdec_set_si(bfdec_t *r, int64_t v) in bfdec_set_si() argument
6592 ret = bfdec_set_ui(r, -v); in bfdec_set_si()
6593 r->sign = 1; in bfdec_set_si()
6595 ret = bfdec_set_ui(r, v); in bfdec_set_si()
6600 static int bfdec_add_internal(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, bf_flags… in bfdec_add_internal() argument
6602 bf_context_t *s = r->ctx; in bfdec_add_internal()
6619 bfdec_set_zero(r, (flags & BF_RND_MASK) == BF_RNDD); in bfdec_add_internal()
6626 bfdec_set_nan(r); in bfdec_add_internal()
6630 bfdec_set_nan(r); in bfdec_add_internal()
6633 bfdec_set_inf(r, a_sign); in bfdec_add_internal()
6637 if (bfdec_set(r, a)) in bfdec_add_internal()
6639 r->sign = a_sign; in bfdec_add_internal()
6654 if (bfdec_resize(r, r_len)) in bfdec_add_internal()
6656 r->sign = a_sign; in bfdec_add_internal()
6657 r->expn = a->expn; in bfdec_add_internal()
6661 r->tab[i] = 0; in bfdec_add_internal()
6663 r->tab[a_offset + i] = a->tab[i]; in bfdec_add_internal()
6680 carry = mp_sub_dec(r->tab + b_offset, r->tab + b_offset, in bfdec_add_internal()
6683 carry = mp_sub_ui_dec(r->tab + b_offset + b1_len, carry, in bfdec_add_internal()
6688 carry = mp_add_dec(r->tab + b_offset, r->tab + b_offset, in bfdec_add_internal()
6691 carry = mp_add_ui_dec(r->tab + b_offset + b1_len, carry, in bfdec_add_internal()
6695 if (bfdec_resize(r, r_len + 1)) { in bfdec_add_internal()
6700 r->tab[r_len] = 1; in bfdec_add_internal()
6701 r->expn += LIMB_DIGITS; in bfdec_add_internal()
6707 ret = bfdec_normalize_and_round(r, prec, flags); in bfdec_add_internal()
6711 bfdec_set_nan(r); in bfdec_add_internal()
6715 static int __bfdec_add(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in __bfdec_add() argument
6718 return bfdec_add_internal(r, a, b, prec, flags, 0); in __bfdec_add()
6721 static int __bfdec_sub(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in __bfdec_sub() argument
6724 return bfdec_add_internal(r, a, b, prec, flags, 1); in __bfdec_sub()
6727 int bfdec_add(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in bfdec_add() argument
6730 return bf_op2((bf_t *)r, (bf_t *)a, (bf_t *)b, prec, flags, in bfdec_add()
6734 int bfdec_sub(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in bfdec_sub() argument
6737 return bf_op2((bf_t *)r, (bf_t *)a, (bf_t *)b, prec, flags, in bfdec_sub()
6741 int bfdec_mul(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in bfdec_mul() argument
6755 bfdec_set_nan(r); in bfdec_mul()
6760 bfdec_set_nan(r); in bfdec_mul()
6763 bfdec_set_inf(r, r_sign); in bfdec_mul()
6767 bfdec_set_zero(r, r_sign); in bfdec_mul()
6780 if (r == a || r == b) { in bfdec_mul()
6781 bfdec_init(r->ctx, &tmp); in bfdec_mul()
6782 r1 = r; in bfdec_mul()
6783 r = &tmp; in bfdec_mul()
6785 if (bfdec_resize(r, a_len + b_len)) { in bfdec_mul()
6786 bfdec_set_nan(r); in bfdec_mul()
6790 mp_mul_basecase_dec(r->tab, a_tab, a_len, b_tab, b_len); in bfdec_mul()
6791 r->sign = r_sign; in bfdec_mul()
6792 r->expn = a->expn + b->expn; in bfdec_mul()
6793 ret = bfdec_normalize_and_round(r, prec, flags); in bfdec_mul()
6795 if (r == &tmp) in bfdec_mul()
6801 int bfdec_mul_si(bfdec_t *r, const bfdec_t *a, int64_t b1, limb_t prec, in bfdec_mul_si() argument
6806 bfdec_init(r->ctx, &b); in bfdec_mul_si()
6808 ret |= bfdec_mul(r, a, &b, prec, flags); in bfdec_mul_si()
6813 int bfdec_add_si(bfdec_t *r, const bfdec_t *a, int64_t b1, limb_t prec, in bfdec_add_si() argument
6819 bfdec_init(r->ctx, &b); in bfdec_add_si()
6821 ret |= bfdec_add(r, a, &b, prec, flags); in bfdec_add_si()
6826 static int __bfdec_div(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, in __bfdec_div() argument
6835 bfdec_set_nan(r); in __bfdec_div()
6838 bfdec_set_nan(r); in __bfdec_div()
6841 bfdec_set_inf(r, r_sign); in __bfdec_div()
6844 bfdec_set_zero(r, r_sign); in __bfdec_div()
6849 bfdec_set_nan(r); in __bfdec_div()
6852 bfdec_set_zero(r, r_sign); in __bfdec_div()
6856 bfdec_set_inf(r, r_sign); in __bfdec_div()
6882 taba = bf_malloc(r->ctx, (na + 1) * sizeof(limb_t)); in __bfdec_div()
6888 if (bfdec_resize(r, n + 1)) in __bfdec_div()
6890 if (mp_div_dec(r->ctx, r->tab, taba, na, b->tab, nb)) { in __bfdec_div()
6892 bf_free(r->ctx, taba); in __bfdec_div()
6900 bf_free(r->ctx, taba); in __bfdec_div()
6903 bfdec_set_nan(r); in __bfdec_div()
6906 r->tab[0] |= 1; in __bfdec_div()
6909 r->expn = a->expn - b->expn + LIMB_DIGITS; in __bfdec_div()
6910 r->sign = r_sign; in __bfdec_div()
6911 ret = bfdec_normalize_and_round(r, prec, flags); in __bfdec_div()
6915 bfdec_set_nan(r); in __bfdec_div()
6919 int bfdec_div(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in bfdec_div() argument
6922 return bf_op2((bf_t *)r, (bf_t *)a, (bf_t *)b, prec, flags, in bfdec_div()
6928 static void bfdec_tdivremu(bf_context_t *s, bfdec_t *q, bfdec_t *r, in bfdec_tdivremu() argument
6933 bfdec_set(r, a); in bfdec_tdivremu()
6936 bfdec_mul(r, q, b, BF_PREC_INF, BF_RNDZ); in bfdec_tdivremu()
6937 bfdec_sub(r, a, r, BF_PREC_INF, BF_RNDZ); in bfdec_tdivremu()
6949 int bfdec_divrem(bfdec_t *q, bfdec_t *r, const bfdec_t *a, const bfdec_t *b, in bfdec_divrem() argument
6960 assert(r != a && r != b); in bfdec_divrem()
6961 assert(q != r); in bfdec_divrem()
6966 bfdec_set_nan(r); in bfdec_divrem()
6969 bfdec_set_nan(r); in bfdec_divrem()
6972 bfdec_set(r, a); in bfdec_divrem()
6973 return bfdec_round(r, prec, flags); in bfdec_divrem()
7013 bfdec_tdivremu(s, q, r, a1, b1); in bfdec_divrem()
7014 if (bfdec_is_nan(q) || bfdec_is_nan(r)) in bfdec_divrem()
7019 if (r->len != 0) { in bfdec_divrem()
7022 if (bfdec_set(r1, r)) in bfdec_divrem()
7039 res |= bfdec_sub(r, r, b1, BF_PREC_INF, BF_RNDZ); in bfdec_divrem()
7045 r->sign ^= a->sign; in bfdec_divrem()
7047 return bfdec_round(r, prec, flags); in bfdec_divrem()
7050 bfdec_set_nan(r); in bfdec_divrem()
7054 int bfdec_rem(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec, in bfdec_rem() argument
7060 bfdec_init(r->ctx, q); in bfdec_rem()
7061 ret = bfdec_divrem(q, r, a, b, prec, flags, rnd_mode); in bfdec_rem()
7067 int bfdec_rint(bfdec_t *r, int rnd_mode) in bfdec_rint() argument
7069 return bfdec_round(r, 0, rnd_mode | BF_FLAG_RADPNT_PREC); in bfdec_rint()
7072 int bfdec_sqrt(bfdec_t *r, const bfdec_t *a, limb_t prec, bf_flags_t flags) in bfdec_sqrt() argument
7080 assert(r != a); in bfdec_sqrt()
7084 bfdec_set_nan(r); in bfdec_sqrt()
7088 bfdec_set(r, a); in bfdec_sqrt()
7093 bfdec_set_nan(r); in bfdec_sqrt()
7104 if (bfdec_resize(r, n)) in bfdec_sqrt()
7129 if (mp_sqrtrem_dec(s, r->tab, a1, n)) { in bfdec_sqrt()
7134 mp_div1_dec(r->tab, r->tab, n, 1 << k, 0); in bfdec_sqrt()
7143 r->tab[0] |= 1; in bfdec_sqrt()
7144 r->sign = 0; in bfdec_sqrt()
7145 r->expn = (a->expn + 1) >> 1; in bfdec_sqrt()
7146 ret = bfdec_round(r, prec, flags); in bfdec_sqrt()
7150 bfdec_set_nan(r); in bfdec_sqrt()
7204 int bfdec_pow_ui(bfdec_t *r, const bfdec_t *a, limb_t b) in bfdec_pow_ui() argument
7208 assert(r != a); in bfdec_pow_ui()
7210 return bfdec_set_ui(r, 1); in bfdec_pow_ui()
7211 ret = bfdec_set(r, a); in bfdec_pow_ui()
7214 ret |= bfdec_mul(r, r, r, BF_PREC_INF, BF_RNDZ); in bfdec_pow_ui()
7216 ret |= bfdec_mul(r, r, a, BF_PREC_INF, BF_RNDZ); in bfdec_pow_ui()
7226 int bfdec_atof(bfdec_t *r, const char *str, const char **pnext, in bfdec_atof() argument
7230 return bf_atof_internal((bf_t *)r, &dummy_exp, str, pnext, 10, prec, in bfdec_atof()
7370 limb_t r; in add_mod() local
7371 r = a + b; in add_mod()
7372 if (r >= m) in add_mod()
7373 r -= m; in add_mod()
7374 return r; in add_mod()
7380 limb_t r; in sub_mod() local
7381 r = a - b; in sub_mod()
7382 if (r > a) in sub_mod()
7383 r += m; in sub_mod()
7384 return r; in sub_mod()
7390 static inline limb_t mod_fast(dlimb_t r, in mod_fast() argument
7395 a1 = r >> NTT_MOD_LOG2_MIN; in mod_fast()
7398 r = r - (dlimb_t)q * m - m * 2; in mod_fast()
7399 r1 = r >> LIMB_BITS; in mod_fast()
7401 r += m & t0; in mod_fast()
7402 r0 = r; in mod_fast()
7403 r1 = r >> LIMB_BITS; in mod_fast()
7413 dlimb_t r; in mul_mod_fast() local
7414 r = (dlimb_t)a * (dlimb_t)b; in mul_mod_fast()
7415 return mod_fast(r, m, m_inv); in mul_mod_fast()
7432 limb_t r, q; in mul_mod_fast2() local
7435 r = a * b - q * m; in mul_mod_fast2()
7436 if (r >= m) in mul_mod_fast2()
7437 r -= m; in mul_mod_fast2()
7438 return r; in mul_mod_fast2()
7447 limb_t r, q; in mul_mod_fast3() local
7450 r = a * b - q * m; in mul_mod_fast3()
7451 return r; in mul_mod_fast3()
7485 static inline __m256d ntt_mod1(__m256d r, __m256d m) in ntt_mod1() argument
7487 return _mm256_blendv_pd(r, r + m, r); in ntt_mod1()
7491 static inline __m256d ntt_mod(__m256d r, __m256d mf, __m256d m2f) in ntt_mod() argument
7493 return _mm256_blendv_pd(r, r + m2f, r) - mf; in ntt_mod()
7500 __m256d r, q, ab1, ab0, qm0, qm1; in ntt_mul_mod() local
7506 r = (ab1 - qm1) + (ab0 - qm0); in ntt_mul_mod()
7507 return r; in ntt_mul_mod()
7981 limb_t base_mask1, a0, a1, a2, r, m, m_inv; in limb_to_ntt() local
8023 r = mod_fast(a, m, m_inv); in limb_to_ntt()
8025 b = ((dlimb_t)r << (LIMB_BITS - NTT_MOD_LOG2_MAX + NTT_MOD_LOG2_MIN)) | a0; in limb_to_ntt()
8026 r = mod_fast(b, m, m_inv); in limb_to_ntt()
8028 tabr[i + j * fft_len] = int_to_ntt_limb(r, m); in limb_to_ntt()
8049 limb_t u[NB_MODS], carry[NB_MODS], fft_len, base_mask1, r; in ntt_to_limb() local
8095 r = (int64_t)y[j].d[p]; in ntt_to_limb()
8097 t = (dlimb_t)u[k] * mods[j] + r; in ntt_to_limb()
8098 r = t >> LIMB_BITS; in ntt_to_limb()
8101 u[l] = r; in ntt_to_limb()
8107 r = (int64_t)y[0].d[p]; in ntt_to_limb()
8109 t = (dlimb_t)u[k] * mods[j] + r + carry[k]; in ntt_to_limb()
8110 r = t >> LIMB_BITS; in ntt_to_limb()
8113 u[l] = r + carry[l]; in ntt_to_limb()
8152 limb_t y[NB_MODS], u[NB_MODS], carry[NB_MODS], fft_len, base_mask1, r; in ntt_to_limb() local
8197 r = y[j]; in ntt_to_limb()
8199 t = (dlimb_t)u[k] * mods[j] + r; in ntt_to_limb()
8200 r = t >> LIMB_BITS; in ntt_to_limb()
8203 u[l] = r; in ntt_to_limb()
8208 r = y[0]; in ntt_to_limb()
8210 t = (dlimb_t)u[k] * mods[j] + r + carry[k]; in ntt_to_limb()
8211 r = t >> LIMB_BITS; in ntt_to_limb()
8214 u[l] = r + carry[l]; in ntt_to_limb()