/* * Tiny arbitrary precision floating point library * * Copyright (c) 2017-2020 Fabrice Bellard * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. */#include <stdlib.h>#include <stdio.h>#include <inttypes.h>#include <math.h>#include <string.h>#include <assert.h> #ifdef __AVX2__#include <immintrin.h>#endif #include "cutils.h"#include "libbf.h" /* enable it to check the multiplication result *///#define USE_MUL_CHECK/* enable it to use FFT/NTT multiplication */#define USE_FFT_MUL/* enable decimal floating point support */#define USE_BF_DEC //#define inline __attribute__((always_inline)) #ifdef __AVX2__#define FFT_MUL_THRESHOLD 100 /* in limbs of the smallest factor */#else#define FFT_MUL_THRESHOLD 100 /* in limbs of the smallest factor */#endif /* XXX: adjust */#define DIVNORM_LARGE_THRESHOLD 50#define UDIV1NORM_THRESHOLD 3 #if LIMB_BITS == 64#define FMT_LIMB1 "%" PRIx64 #define FMT_LIMB "%016" PRIx64 #define PRId_LIMB PRId64#define PRIu_LIMB PRIu64 #else #define FMT_LIMB1 "%x"#define FMT_LIMB "%08x"#define PRId_LIMB "d"#define PRIu_LIMB "u" #endif typedef intptr_t mp_size_t; typedef int bf_op2_func_t(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags); #ifdef USE_FFT_MUL #define FFT_MUL_R_OVERLAP_A (1 << 0)#define FFT_MUL_R_OVERLAP_B (1 << 1)#define FFT_MUL_R_NORESIZE (1 << 2) static no_inline int fft_mul(bf_context_t *s, bf_t *res, limb_t *a_tab, limb_t a_len, limb_t *b_tab, limb_t b_len, int mul_flags);static void fft_clear_cache(bf_context_t *s);#endif#ifdef USE_BF_DECstatic limb_t get_digit(const limb_t *tab, limb_t len, slimb_t pos);#endif /* could leading zeros */static inline int clz(limb_t a){ if (a == 0) { return LIMB_BITS; } else {#if LIMB_BITS == 64 return clz64(a);#else return clz32(a);#endif }} static inline int ctz(limb_t a){ if (a == 0) { return LIMB_BITS; } else {#if LIMB_BITS == 64 return ctz64(a);#else return ctz32(a);#endif }} static inline int ceil_log2(limb_t a){ if (a <= 1) return 0; else return LIMB_BITS - clz(a - 1);} /* b must be >= 1 */static inline slimb_t ceil_div(slimb_t a, slimb_t b){ if (a >= 0) return (a + b - 1) / b; else return a / b;} /* b must be >= 1 */static inline slimb_t floor_div(slimb_t a, slimb_t b){ if (a >= 0) { return a / b; } else { return (a - b + 1) / b; }} /* return r = a modulo b (0 <= r <= b - 1. b must be >= 1 */static inline limb_t smod(slimb_t a, slimb_t b){ a = a % (slimb_t)b; if (a < 0) a += b; return a;} /* signed addition with saturation */static inline slimb_t sat_add(slimb_t a, slimb_t b){ slimb_t r; r = a + b; /* overflow ? */ if (((a ^ r) & (b ^ r)) < 0) r = (a >> (LIMB_BITS - 1)) ^ (((limb_t)1 << (LIMB_BITS - 1)) - 1); return r;} #define malloc(s) malloc_is_forbidden(s)#define free(p) free_is_forbidden(p)#define realloc(p, s) realloc_is_forbidden(p, s) void bf_context_init(bf_context_t *s, bf_realloc_func_t *realloc_func, void *realloc_opaque){ memset(s, 0, sizeof(*s)); s->realloc_func = realloc_func; s->realloc_opaque = realloc_opaque;} void bf_context_end(bf_context_t *s){ bf_clear_cache(s);} void bf_init(bf_context_t *s, bf_t *r){ r->ctx = s; r->sign = 0; r->expn = BF_EXP_ZERO; r->len = 0; r->tab = NULL;} /* return 0 if OK, -1 if alloc error */int bf_resize(bf_t *r, limb_t len){ limb_t *tab; if (len != r->len) { tab = bf_realloc(r->ctx, r->tab, len * sizeof(limb_t)); if (!tab && len != 0) return -1; r->tab = tab; r->len = len; } return 0;} /* return 0 or BF_ST_MEM_ERROR */int bf_set_ui(bf_t *r, uint64_t a){ r->sign = 0; if (a == 0) { r->expn = BF_EXP_ZERO; bf_resize(r, 0); /* cannot fail */ } #if LIMB_BITS == 32 else if (a <= 0xffffffff)#else else#endif { int shift; if (bf_resize(r, 1)) goto fail; shift = clz(a); r->tab[0] = a << shift; r->expn = LIMB_BITS - shift; }#if LIMB_BITS == 32 else { uint32_t a1, a0; int shift; if (bf_resize(r, 2)) goto fail; a0 = a; a1 = a >> 32; shift = clz(a1); r->tab[0] = a0 << shift; r->tab[1] = (a1 << shift) | (a0 >> (LIMB_BITS - shift)); r->expn = 2 * LIMB_BITS - shift; }#endif return 0; fail: bf_set_nan(r); return BF_ST_MEM_ERROR;} /* return 0 or BF_ST_MEM_ERROR */int bf_set_si(bf_t *r, int64_t a){ int ret; if (a < 0) { ret = bf_set_ui(r, -a); r->sign = 1; } else { ret = bf_set_ui(r, a); } return ret;} void bf_set_nan(bf_t *r){ bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_NAN; r->sign = 0;} void bf_set_zero(bf_t *r, int is_neg){ bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_ZERO; r->sign = is_neg;} void bf_set_inf(bf_t *r, int is_neg){ bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_INF; r->sign = is_neg;} /* return 0 or BF_ST_MEM_ERROR */int bf_set(bf_t *r, const bf_t *a){ if (r == a) return 0; if (bf_resize(r, a->len)) { bf_set_nan(r); return BF_ST_MEM_ERROR; } r->sign = a->sign; r->expn = a->expn; memcpy(r->tab, a->tab, a->len * sizeof(limb_t)); return 0;} /* equivalent to bf_set(r, a); bf_delete(a) */void bf_move(bf_t *r, bf_t *a){ bf_context_t *s = r->ctx; if (r == a) return; bf_free(s, r->tab); *r = *a;} static limb_t get_limbz(const bf_t *a, limb_t idx){ if (idx >= a->len) return 0; else return a->tab[idx];} /* get LIMB_BITS at bit position 'pos' in tab */static inline limb_t get_bits(const limb_t *tab, limb_t len, slimb_t pos){ limb_t i, a0, a1; int p; i = pos >> LIMB_LOG2_BITS; p = pos & (LIMB_BITS - 1); if (i < len) a0 = tab[i]; else a0 = 0; if (p == 0) { return a0; } else { i++; if (i < len) a1 = tab[i]; else a1 = 0; return (a0 >> p) | (a1 << (LIMB_BITS - p)); }} static inline limb_t get_bit(const limb_t *tab, limb_t len, slimb_t pos){ slimb_t i; i = pos >> LIMB_LOG2_BITS; if (i < 0 || i >= len) return 0; return (tab[i] >> (pos & (LIMB_BITS - 1))) & 1;} static inline limb_t limb_mask(int start, int last){ limb_t v; int n; n = last - start + 1; if (n == LIMB_BITS) v = -1; else v = (((limb_t)1 << n) - 1) << start; return v;} static limb_t mp_scan_nz(const limb_t *tab, mp_size_t n){ mp_size_t i; for(i = 0; i < n; i++) { if (tab[i] != 0) return 1; } return 0;} /* return != 0 if one bit between 0 and bit_pos inclusive is not zero. */static inline limb_t scan_bit_nz(const bf_t *r, slimb_t bit_pos){ slimb_t pos; limb_t v; pos = bit_pos >> LIMB_LOG2_BITS; if (pos < 0) return 0; v = r->tab[pos] & limb_mask(0, bit_pos & (LIMB_BITS - 1)); if (v != 0) return 1; pos--; while (pos >= 0) { if (r->tab[pos] != 0) return 1; pos--; } return 0;} /* return the addend for rounding. Note that prec can be <= 0 (for BF_FLAG_RADPNT_PREC) */static int bf_get_rnd_add(int *pret, const bf_t *r, limb_t l, slimb_t prec, int rnd_mode){ int add_one, inexact; limb_t bit1, bit0; if (rnd_mode == BF_RNDF) { bit0 = 1; /* faithful rounding does not honor the INEXACT flag */ } else { /* starting limb for bit 'prec + 1' */ bit0 = scan_bit_nz(r, l * LIMB_BITS - 1 - bf_max(0, prec + 1)); } /* get the bit at 'prec' */ bit1 = get_bit(r->tab, l, l * LIMB_BITS - 1 - prec); inexact = (bit1 | bit0) != 0; add_one = 0; switch(rnd_mode) { case BF_RNDZ: break; case BF_RNDN: if (bit1) { if (bit0) { add_one = 1; } else { /* round to even */ add_one = get_bit(r->tab, l, l * LIMB_BITS - 1 - (prec - 1)); } } break; case BF_RNDD: case BF_RNDU: if (r->sign == (rnd_mode == BF_RNDD)) add_one = inexact; break; case BF_RNDA: add_one = inexact; break; case BF_RNDNA: case BF_RNDF: add_one = bit1; break; default: abort(); } if (inexact) *pret |= BF_ST_INEXACT; return add_one;} static int bf_set_overflow(bf_t *r, int sign, limb_t prec, bf_flags_t flags){ slimb_t i, l, e_max; int rnd_mode; rnd_mode = flags & BF_RND_MASK; if (prec == BF_PREC_INF || rnd_mode == BF_RNDN || rnd_mode == BF_RNDNA || rnd_mode == BF_RNDA || (rnd_mode == BF_RNDD && sign == 1) || (rnd_mode == BF_RNDU && sign == 0)) { bf_set_inf(r, sign); } else { /* set to maximum finite number */ l = (prec + LIMB_BITS - 1) / LIMB_BITS; if (bf_resize(r, l)) { bf_set_nan(r); return BF_ST_MEM_ERROR; } r->tab[0] = limb_mask((-prec) & (LIMB_BITS - 1), LIMB_BITS - 1); for(i = 1; i < l; i++) r->tab[i] = (limb_t)-1; e_max = (limb_t)1 << (bf_get_exp_bits(flags) - 1); r->expn = e_max; r->sign = sign; } return BF_ST_OVERFLOW | BF_ST_INEXACT;} /* round to prec1 bits assuming 'r' is non zero and finite. 'r' is assumed to have length 'l' (1 <= l <= r->len). Note: 'prec1' can be infinite (BF_PREC_INF). 'ret' is 0 or BF_ST_INEXACT if the result is known to be inexact. Can fail with BF_ST_MEM_ERROR in case of overflow not returning infinity. */static int __bf_round(bf_t *r, limb_t prec1, bf_flags_t flags, limb_t l, int ret){ limb_t v, a; int shift, add_one, rnd_mode; slimb_t i, bit_pos, pos, e_min, e_max, e_range, prec; /* e_min and e_max are computed to match the IEEE 754 conventions */ e_range = (limb_t)1 << (bf_get_exp_bits(flags) - 1); e_min = -e_range + 3; e_max = e_range; if (flags & BF_FLAG_RADPNT_PREC) { /* 'prec' is the precision after the radix point */ if (prec1 != BF_PREC_INF) prec = r->expn + prec1; elseLines 1–500 of 8464. Load more
/* * Tiny arbitrary precision floating point library * * Copyright (c) 2017-2020 Fabrice Bellard * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. */ #include <stdlib.h> #include <stdio.h> #include <inttypes.h> #include <math.h> #include <string.h> #include <assert.h> #ifdef __AVX2__ #include <immintrin.h> #endif #include "cutils.h" #include "libbf.h" /* enable it to check the multiplication result */ //#define USE_MUL_CHECK /* enable it to use FFT/NTT multiplication */ #define USE_FFT_MUL /* enable decimal floating point support */ #define USE_BF_DEC //#define inline __attribute__((always_inline)) #ifdef __AVX2__ #define FFT_MUL_THRESHOLD 100 /* in limbs of the smallest factor */ #else #define FFT_MUL_THRESHOLD 100 /* in limbs of the smallest factor */ #endif /* XXX: adjust */ #define DIVNORM_LARGE_THRESHOLD 50 #define UDIV1NORM_THRESHOLD 3 #if LIMB_BITS == 64 #define FMT_LIMB1 "%" PRIx64 #define FMT_LIMB "%016" PRIx64 #define PRId_LIMB PRId64 #define PRIu_LIMB PRIu64 #else #define FMT_LIMB1 "%x" #define FMT_LIMB "%08x" #define PRId_LIMB "d" #define PRIu_LIMB "u" #endif typedef intptr_t mp_size_t; typedef int bf_op2_func_t(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags); #ifdef USE_FFT_MUL #define FFT_MUL_R_OVERLAP_A (1 << 0) #define FFT_MUL_R_OVERLAP_B (1 << 1) #define FFT_MUL_R_NORESIZE (1 << 2) static no_inline int fft_mul(bf_context_t *s, bf_t *res, limb_t *a_tab, limb_t a_len, limb_t *b_tab, limb_t b_len, int mul_flags); static void fft_clear_cache(bf_context_t *s); #endif #ifdef USE_BF_DEC static limb_t get_digit(const limb_t *tab, limb_t len, slimb_t pos); #endif /* could leading zeros */ static inline int clz(limb_t a) { if (a == 0) { return LIMB_BITS; } else { #if LIMB_BITS == 64 return clz64(a); #else return clz32(a); #endif } } static inline int ctz(limb_t a) { if (a == 0) { return LIMB_BITS; } else { #if LIMB_BITS == 64 return ctz64(a); #else return ctz32(a); #endif } } static inline int ceil_log2(limb_t a) { if (a <= 1) return 0; else return LIMB_BITS - clz(a - 1); } /* b must be >= 1 */ static inline slimb_t ceil_div(slimb_t a, slimb_t b) { if (a >= 0) return (a + b - 1) / b; else return a / b; } /* b must be >= 1 */ static inline slimb_t floor_div(slimb_t a, slimb_t b) { if (a >= 0) { return a / b; } else { return (a - b + 1) / b; } } /* return r = a modulo b (0 <= r <= b - 1. b must be >= 1 */ static inline limb_t smod(slimb_t a, slimb_t b) { a = a % (slimb_t)b; if (a < 0) a += b; return a; } /* signed addition with saturation */ static inline slimb_t sat_add(slimb_t a, slimb_t b) { slimb_t r; r = a + b; /* overflow ? */ if (((a ^ r) & (b ^ r)) < 0) r = (a >> (LIMB_BITS - 1)) ^ (((limb_t)1 << (LIMB_BITS - 1)) - 1); return r; } #define malloc(s) malloc_is_forbidden(s) #define free(p) free_is_forbidden(p) #define realloc(p, s) realloc_is_forbidden(p, s) void bf_context_init(bf_context_t *s, bf_realloc_func_t *realloc_func, void *realloc_opaque) { memset(s, 0, sizeof(*s)); s->realloc_func = realloc_func; s->realloc_opaque = realloc_opaque; } void bf_context_end(bf_context_t *s) { bf_clear_cache(s); } void bf_init(bf_context_t *s, bf_t *r) { r->ctx = s; r->sign = 0; r->expn = BF_EXP_ZERO; r->len = 0; r->tab = NULL; } /* return 0 if OK, -1 if alloc error */ int bf_resize(bf_t *r, limb_t len) { limb_t *tab; if (len != r->len) { tab = bf_realloc(r->ctx, r->tab, len * sizeof(limb_t)); if (!tab && len != 0) return -1; r->tab = tab; r->len = len; } return 0; } /* return 0 or BF_ST_MEM_ERROR */ int bf_set_ui(bf_t *r, uint64_t a) { r->sign = 0; if (a == 0) { r->expn = BF_EXP_ZERO; bf_resize(r, 0); /* cannot fail */ } #if LIMB_BITS == 32 else if (a <= 0xffffffff) #else else #endif { int shift; if (bf_resize(r, 1)) goto fail; shift = clz(a); r->tab[0] = a << shift; r->expn = LIMB_BITS - shift; } #if LIMB_BITS == 32 else { uint32_t a1, a0; int shift; if (bf_resize(r, 2)) goto fail; a0 = a; a1 = a >> 32; shift = clz(a1); r->tab[0] = a0 << shift; r->tab[1] = (a1 << shift) | (a0 >> (LIMB_BITS - shift)); r->expn = 2 * LIMB_BITS - shift; } #endif return 0; fail: bf_set_nan(r); return BF_ST_MEM_ERROR; } /* return 0 or BF_ST_MEM_ERROR */ int bf_set_si(bf_t *r, int64_t a) { int ret; if (a < 0) { ret = bf_set_ui(r, -a); r->sign = 1; } else { ret = bf_set_ui(r, a); } return ret; } void bf_set_nan(bf_t *r) { bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_NAN; r->sign = 0; } void bf_set_zero(bf_t *r, int is_neg) { bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_ZERO; r->sign = is_neg; } void bf_set_inf(bf_t *r, int is_neg) { bf_resize(r, 0); /* cannot fail */ r->expn = BF_EXP_INF; r->sign = is_neg; } /* return 0 or BF_ST_MEM_ERROR */ int bf_set(bf_t *r, const bf_t *a) { if (r == a) return 0; if (bf_resize(r, a->len)) { bf_set_nan(r); return BF_ST_MEM_ERROR; } r->sign = a->sign; r->expn = a->expn; memcpy(r->tab, a->tab, a->len * sizeof(limb_t)); return 0; } /* equivalent to bf_set(r, a); bf_delete(a) */ void bf_move(bf_t *r, bf_t *a) { bf_context_t *s = r->ctx; if (r == a) return; bf_free(s, r->tab); *r = *a; } static limb_t get_limbz(const bf_t *a, limb_t idx) { if (idx >= a->len) return 0; else return a->tab[idx]; } /* get LIMB_BITS at bit position 'pos' in tab */ static inline limb_t get_bits(const limb_t *tab, limb_t len, slimb_t pos) { limb_t i, a0, a1; int p; i = pos >> LIMB_LOG2_BITS; p = pos & (LIMB_BITS - 1); if (i < len) a0 = tab[i]; else a0 = 0; if (p == 0) { return a0; } else { i++; if (i < len) a1 = tab[i]; else a1 = 0; return (a0 >> p) | (a1 << (LIMB_BITS - p)); } } static inline limb_t get_bit(const limb_t *tab, limb_t len, slimb_t pos) { slimb_t i; i = pos >> LIMB_LOG2_BITS; if (i < 0 || i >= len) return 0; return (tab[i] >> (pos & (LIMB_BITS - 1))) & 1; } static inline limb_t limb_mask(int start, int last) { limb_t v; int n; n = last - start + 1; if (n == LIMB_BITS) v = -1; else v = (((limb_t)1 << n) - 1) << start; return v; } static limb_t mp_scan_nz(const limb_t *tab, mp_size_t n) { mp_size_t i; for(i = 0; i < n; i++) { if (tab[i] != 0) return 1; } return 0; } /* return != 0 if one bit between 0 and bit_pos inclusive is not zero. */ static inline limb_t scan_bit_nz(const bf_t *r, slimb_t bit_pos) { slimb_t pos; limb_t v; pos = bit_pos >> LIMB_LOG2_BITS; if (pos < 0) return 0; v = r->tab[pos] & limb_mask(0, bit_pos & (LIMB_BITS - 1)); if (v != 0) return 1; pos--; while (pos >= 0) { if (r->tab[pos] != 0) return 1; pos--; } return 0; } /* return the addend for rounding. Note that prec can be <= 0 (for BF_FLAG_RADPNT_PREC) */ static int bf_get_rnd_add(int *pret, const bf_t *r, limb_t l, slimb_t prec, int rnd_mode) { int add_one, inexact; limb_t bit1, bit0; if (rnd_mode == BF_RNDF) { bit0 = 1; /* faithful rounding does not honor the INEXACT flag */ } else { /* starting limb for bit 'prec + 1' */ bit0 = scan_bit_nz(r, l * LIMB_BITS - 1 - bf_max(0, prec + 1)); } /* get the bit at 'prec' */ bit1 = get_bit(r->tab, l, l * LIMB_BITS - 1 - prec); inexact = (bit1 | bit0) != 0; add_one = 0; switch(rnd_mode) { case BF_RNDZ: break; case BF_RNDN: if (bit1) { if (bit0) { add_one = 1; } else { /* round to even */ add_one = get_bit(r->tab, l, l * LIMB_BITS - 1 - (prec - 1)); } } break; case BF_RNDD: case BF_RNDU: if (r->sign == (rnd_mode == BF_RNDD)) add_one = inexact; break; case BF_RNDA: add_one = inexact; break; case BF_RNDNA: case BF_RNDF: add_one = bit1; break; default: abort(); } if (inexact) *pret |= BF_ST_INEXACT; return add_one; } static int bf_set_overflow(bf_t *r, int sign, limb_t prec, bf_flags_t flags) { slimb_t i, l, e_max; int rnd_mode; rnd_mode = flags & BF_RND_MASK; if (prec == BF_PREC_INF || rnd_mode == BF_RNDN || rnd_mode == BF_RNDNA || rnd_mode == BF_RNDA || (rnd_mode == BF_RNDD && sign == 1) || (rnd_mode == BF_RNDU && sign == 0)) { bf_set_inf(r, sign); } else { /* set to maximum finite number */ l = (prec + LIMB_BITS - 1) / LIMB_BITS; if (bf_resize(r, l)) { bf_set_nan(r); return BF_ST_MEM_ERROR; } r->tab[0] = limb_mask((-prec) & (LIMB_BITS - 1), LIMB_BITS - 1); for(i = 1; i < l; i++) r->tab[i] = (limb_t)-1; e_max = (limb_t)1 << (bf_get_exp_bits(flags) - 1); r->expn = e_max; r->sign = sign; } return BF_ST_OVERFLOW | BF_ST_INEXACT; } /* round to prec1 bits assuming 'r' is non zero and finite. 'r' is assumed to have length 'l' (1 <= l <= r->len). Note: 'prec1' can be infinite (BF_PREC_INF). 'ret' is 0 or BF_ST_INEXACT if the result is known to be inexact. Can fail with BF_ST_MEM_ERROR in case of overflow not returning infinity. */ static int __bf_round(bf_t *r, limb_t prec1, bf_flags_t flags, limb_t l, int ret) { limb_t v, a; int shift, add_one, rnd_mode; slimb_t i, bit_pos, pos, e_min, e_max, e_range, prec; /* e_min and e_max are computed to match the IEEE 754 conventions */ e_range = (limb_t)1 << (bf_get_exp_bits(flags) - 1); e_min = -e_range + 3; e_max = e_range; if (flags & BF_FLAG_RADPNT_PREC) { /* 'prec' is the precision after the radix point */ if (prec1 != BF_PREC_INF) prec = r->expn + prec1; else