Import libbf https://bellard.org/libbf/ as third_party(2020-01-19 version)

Signed-off-by: Seonghyun Kim <sh8281.kim@samsung.com>
This commit is contained in:
Seonghyun Kim 2020-10-12 16:44:55 +09:00 committed by Hyukwoo Park
commit b03e0ab040
8 changed files with 9461 additions and 1 deletions

View file

@ -89,6 +89,31 @@ ADD_SUBDIRECTORY (third_party/GCutil)
SET (ESCARGOT_LIBRARIES ${ESCARGOT_LIBRARIES} gc-lib)
# LIBBF
ADD_LIBRARY (libbf STATIC
${ESCARGOT_THIRD_PARTY_ROOT}/libbf/libbf.c
${ESCARGOT_THIRD_PARTY_ROOT}/libbf/cutils.c)
TARGET_INCLUDE_DIRECTORIES (libbf PUBLIC ${ESCARGOT_THIRD_PARTY_ROOT}/libbf)
SET (LIBBF_CFLAGS
${ESCARGOT_GCUTIL_CFLAGS} # we can share arch flags with gcutil
${CFLAGS_FROM_ENV}
-w
-g3
-fdata-sections
-ffunction-sections
-fno-omit-frame-pointer
-fvisibility=hidden)
IF (${ESCARGOT_MODE} STREQUAL "debug")
SET (LIBBF_CFLAGS ${ESCARGOT_CXXFLAGS_DEBUG} ${LIBBF_CFLAGS})
ELSEIF (${ESCARGOT_MODE} STREQUAL "release")
SET (LIBBF_CFLAGS ${ESCARGOT_CXXFLAGS_RELEASE} ${LIBBF_CFLAGS})
ENDIF()
TARGET_COMPILE_OPTIONS (libbf PRIVATE ${LIBBF_CFLAGS})
SET (ESCARGOT_LIBRARIES ${ESCARGOT_LIBRARIES} libbf)
# RUNTIME ICU BINDER
SET (RIB_CFLAGS ${ESCARGOT_CXXFLAGS})

1
third_party/libbf/VERSION vendored Normal file
View file

@ -0,0 +1 @@
2020-01-19

178
third_party/libbf/cutils.c vendored Normal file
View file

@ -0,0 +1,178 @@
/*
* C utilities
*
* Copyright (c) 2017 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 <stdarg.h>
#include <string.h>
#include "cutils.h"
void pstrcpy(char *buf, int buf_size, const char *str)
{
int c;
char *q = buf;
if (buf_size <= 0)
return;
for(;;) {
c = *str++;
if (c == 0 || q >= buf + buf_size - 1)
break;
*q++ = c;
}
*q = '\0';
}
/* strcat and truncate. */
char *pstrcat(char *buf, int buf_size, const char *s)
{
int len;
len = strlen(buf);
if (len < buf_size)
pstrcpy(buf + len, buf_size - len, s);
return buf;
}
int strstart(const char *str, const char *val, const char **ptr)
{
const char *p, *q;
p = str;
q = val;
while (*q != '\0') {
if (*p != *q)
return 0;
p++;
q++;
}
if (ptr)
*ptr = p;
return 1;
}
void dbuf_init2(DynBuf *s, void *opaque, DynBufReallocFunc *realloc_func)
{
memset(s, 0, sizeof(*s));
s->opaque = opaque;
s->realloc_func = realloc_func;
}
static void *dbuf_default_realloc(void *opaque, void *ptr, size_t size)
{
return realloc(ptr, size);
}
void dbuf_init(DynBuf *s)
{
dbuf_init2(s, NULL, dbuf_default_realloc);
}
/* return < 0 if error */
int dbuf_realloc(DynBuf *s, size_t new_size)
{
size_t size;
uint8_t *new_buf;
if (new_size > s->allocated_size) {
if (s->error)
return -1;
size = s->allocated_size * 3 / 2;
if (size > new_size)
new_size = size;
new_buf = s->realloc_func(s->opaque, s->buf, new_size);
if (!new_buf) {
s->error = TRUE;
return -1;
}
s->buf = new_buf;
s->allocated_size = new_size;
}
return 0;
}
int dbuf_write(DynBuf *s, size_t offset, const uint8_t *data, size_t len)
{
size_t end;
end = offset + len;
if (dbuf_realloc(s, end))
return -1;
memcpy(s->buf + offset, data, len);
if (end > s->size)
s->size = end;
return 0;
}
int dbuf_put(DynBuf *s, const uint8_t *data, size_t len)
{
if (unlikely((s->size + len) > s->allocated_size)) {
if (dbuf_realloc(s, s->size + len))
return -1;
}
memcpy(s->buf + s->size, data, len);
s->size += len;
return 0;
}
int dbuf_putc(DynBuf *s, uint8_t c)
{
return dbuf_put(s, &c, 1);
}
int dbuf_putstr(DynBuf *s, const char *str)
{
return dbuf_put(s, (const uint8_t *)str, strlen(str));
}
int __attribute__((format(printf, 2, 3))) dbuf_printf(DynBuf *s,
const char *fmt, ...)
{
va_list ap;
char buf[128];
int len;
va_start(ap, fmt);
len = vsnprintf(buf, sizeof(buf), fmt, ap);
va_end(ap);
if (len < sizeof(buf)) {
/* fast case */
return dbuf_put(s, (uint8_t *)buf, len);
} else {
if (dbuf_realloc(s, s->size + len + 1))
return -1;
va_start(ap, fmt);
vsnprintf((char *)(s->buf + s->size), s->allocated_size - s->size,
fmt, ap);
va_end(ap);
s->size += len;
}
return 0;
}
void dbuf_free(DynBuf *s)
{
/* we test s->buf as a fail safe to avoid crashing if dbuf_free()
is called twice */
if (s->buf) {
s->realloc_func(s->opaque, s->buf, 0);
}
memset(s, 0, sizeof(*s));
}

157
third_party/libbf/cutils.h vendored Normal file
View file

@ -0,0 +1,157 @@
/*
* C utilities
*
* Copyright (c) 2017 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.
*/
#ifndef _CUTILS_H
#define _CUTILS_H
#include <inttypes.h>
#define likely(x) __builtin_expect(!!(x), 1)
#define unlikely(x) __builtin_expect(!!(x), 0)
#define force_inline inline __attribute__((always_inline))
#define no_inline __attribute__((noinline))
#define __maybe_unused __attribute__((unused))
#define xglue(x, y) x ## y
#define glue(x, y) xglue(x, y)
#define stringify(s) tostring(s)
#define tostring(s) #s
#ifndef offsetof
#define offsetof(type, field) ((size_t) &((type *)0)->field)
#endif
#define countof(x) (sizeof(x) / sizeof(x[0]))
typedef int BOOL;
enum {
FALSE = 0,
TRUE = 1,
};
void pstrcpy(char *buf, int buf_size, const char *str);
char *pstrcat(char *buf, int buf_size, const char *s);
int strstart(const char *str, const char *val, const char **ptr);
static inline int max_int(int a, int b)
{
if (a > b)
return a;
else
return b;
}
static inline int min_int(int a, int b)
{
if (a < b)
return a;
else
return b;
}
/* WARNING: undefined if a = 0 */
static inline int clz32(unsigned int a)
{
return __builtin_clz(a);
}
/* WARNING: undefined if a = 0 */
static inline int clz64(uint64_t a)
{
return __builtin_clzll(a);
}
/* WARNING: undefined if a = 0 */
static inline int ctz32(unsigned int a)
{
return __builtin_ctz(a);
}
/* WARNING: undefined if a = 0 */
static inline int ctz64(uint64_t a)
{
return __builtin_ctzll(a);
}
struct __attribute__((packed)) packed_u32 {
uint32_t v;
};
struct __attribute__((packed)) packed_u16 {
uint16_t v;
};
static inline uint32_t get_u32(const uint8_t *tab)
{
return ((struct packed_u32 *)tab)->v;
}
static inline void put_u32(uint8_t *tab, uint32_t val)
{
((struct packed_u32 *)tab)->v = val;
}
static inline uint32_t get_u16(const uint8_t *tab)
{
return ((struct packed_u16 *)tab)->v;
}
static inline void put_u16(uint8_t *tab, uint16_t val)
{
((struct packed_u16 *)tab)->v = val;
}
typedef void *DynBufReallocFunc(void *opaque, void *ptr, size_t size);
typedef struct {
uint8_t *buf;
size_t size;
size_t allocated_size;
BOOL error; /* true if a memory allocation error occurred */
DynBufReallocFunc *realloc_func;
void *opaque; /* for realloc_func */
} DynBuf;
void dbuf_init(DynBuf *s);
void dbuf_init2(DynBuf *s, void *opaque, DynBufReallocFunc *realloc_func);
int dbuf_realloc(DynBuf *s, size_t new_size);
int dbuf_write(DynBuf *s, size_t offset, const uint8_t *data, size_t len);
int dbuf_put(DynBuf *s, const uint8_t *data, size_t len);
int dbuf_putc(DynBuf *s, uint8_t c);
int dbuf_putstr(DynBuf *s, const char *str);
static inline int dbuf_put_u32(DynBuf *s, uint32_t val)
{
return dbuf_put(s, (uint8_t *)&val, 4);
}
static inline int dbuf_put_u16(DynBuf *s, uint16_t val)
{
return dbuf_put(s, (uint8_t *)&val, 2);
}
int __attribute__((format(printf, 2, 3))) dbuf_printf(DynBuf *s,
const char *fmt, ...);
void dbuf_free(DynBuf *s);
static inline BOOL dbuf_error(DynBuf *s) {
return s->error;
}
#endif

8404
third_party/libbf/libbf.c vendored Normal file

File diff suppressed because it is too large Load diff

534
third_party/libbf/libbf.h vendored Normal file
View file

@ -0,0 +1,534 @@
/*
* 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.
*/
#ifndef LIBBF_H
#define LIBBF_H
#include <stddef.h>
#include <stdint.h>
#if defined(__x86_64__)
#define LIMB_LOG2_BITS 6
#else
#define LIMB_LOG2_BITS 5
#endif
#define LIMB_BITS (1 << LIMB_LOG2_BITS)
#if LIMB_BITS == 64
typedef __int128 int128_t;
typedef unsigned __int128 uint128_t;
typedef int64_t slimb_t;
typedef uint64_t limb_t;
typedef uint128_t dlimb_t;
#define BF_RAW_EXP_MIN INT64_MIN
#define BF_RAW_EXP_MAX INT64_MAX
#define LIMB_DIGITS 19
#define BF_DEC_BASE UINT64_C(10000000000000000000)
#else
typedef int32_t slimb_t;
typedef uint32_t limb_t;
typedef uint64_t dlimb_t;
#define BF_RAW_EXP_MIN INT32_MIN
#define BF_RAW_EXP_MAX INT32_MAX
#define LIMB_DIGITS 9
#define BF_DEC_BASE 1000000000U
#endif
/* in bits */
/* minimum number of bits for the exponent */
#define BF_EXP_BITS_MIN 3
/* maximum number of bits for the exponent */
#define BF_EXP_BITS_MAX (LIMB_BITS - 3)
/* extended range for exponent, used internally */
#define BF_EXT_EXP_BITS_MAX (BF_EXP_BITS_MAX + 1)
/* minimum possible precision */
#define BF_PREC_MIN 2
/* minimum possible precision */
#define BF_PREC_MAX (((limb_t)1 << (LIMB_BITS - 2)) - 2)
/* some operations support infinite precision */
#define BF_PREC_INF (BF_PREC_MAX + 1) /* infinite precision */
#if LIMB_BITS == 64
#define BF_CHKSUM_MOD (UINT64_C(975620677) * UINT64_C(9795002197))
#else
#define BF_CHKSUM_MOD 975620677U
#endif
#define BF_EXP_ZERO BF_RAW_EXP_MIN
#define BF_EXP_INF (BF_RAW_EXP_MAX - 1)
#define BF_EXP_NAN BF_RAW_EXP_MAX
/* +/-zero is represented with expn = BF_EXP_ZERO and len = 0,
+/-infinity is represented with expn = BF_EXP_INF and len = 0,
NaN is represented with expn = BF_EXP_NAN and len = 0 (sign is ignored)
*/
typedef struct {
struct bf_context_t *ctx;
int sign;
slimb_t expn;
limb_t len;
limb_t *tab;
} bf_t;
typedef struct {
/* must be kept identical to bf_t */
struct bf_context_t *ctx;
int sign;
slimb_t expn;
limb_t len;
limb_t *tab;
} bfdec_t;
typedef enum {
BF_RNDN, /* round to nearest, ties to even */
BF_RNDZ, /* round to zero */
BF_RNDD, /* round to -inf (the code relies on (BF_RNDD xor BF_RNDU) = 1) */
BF_RNDU, /* round to +inf */
BF_RNDNA, /* round to nearest, ties away from zero */
BF_RNDA, /* round away from zero */
BF_RNDF, /* faithful rounding (nondeterministic, either RNDD or RNDU,
inexact flag is always set) */
} bf_rnd_t;
/* allow subnormal numbers. Only available if the number of exponent
bits is <= BF_EXP_BITS_USER_MAX and prec != BF_PREC_INF. */
#define BF_FLAG_SUBNORMAL (1 << 3)
/* 'prec' is the precision after the radix point instead of the whole
mantissa. Can only be used with bf_round() and
bfdec_[add|sub|mul|div|sqrt|round](). */
#define BF_FLAG_RADPNT_PREC (1 << 4)
#define BF_RND_MASK 0x7
#define BF_EXP_BITS_SHIFT 5
#define BF_EXP_BITS_MASK 0x3f
/* shortcut for bf_set_exp_bits(BF_EXT_EXP_BITS_MAX) */
#define BF_FLAG_EXT_EXP (BF_EXP_BITS_MASK << BF_EXP_BITS_SHIFT)
/* contains the rounding mode and number of exponents bits */
typedef uint32_t bf_flags_t;
typedef void *bf_realloc_func_t(void *opaque, void *ptr, size_t size);
typedef struct {
bf_t val;
limb_t prec;
} BFConstCache;
typedef struct bf_context_t {
void *realloc_opaque;
bf_realloc_func_t *realloc_func;
BFConstCache log2_cache;
BFConstCache pi_cache;
struct BFNTTState *ntt_state;
} bf_context_t;
static inline int bf_get_exp_bits(bf_flags_t flags)
{
int e;
e = (flags >> BF_EXP_BITS_SHIFT) & BF_EXP_BITS_MASK;
if (e == BF_EXP_BITS_MASK)
return BF_EXP_BITS_MAX + 1;
else
return BF_EXP_BITS_MAX - e;
}
static inline bf_flags_t bf_set_exp_bits(int n)
{
return ((BF_EXP_BITS_MAX - n) & BF_EXP_BITS_MASK) << BF_EXP_BITS_SHIFT;
}
/* returned status */
#define BF_ST_INVALID_OP (1 << 0)
#define BF_ST_DIVIDE_ZERO (1 << 1)
#define BF_ST_OVERFLOW (1 << 2)
#define BF_ST_UNDERFLOW (1 << 3)
#define BF_ST_INEXACT (1 << 4)
/* indicate that a memory allocation error occured. NaN is returned */
#define BF_ST_MEM_ERROR (1 << 5)
#define BF_RADIX_MAX 36 /* maximum radix for bf_atof() and bf_ftoa() */
static inline slimb_t bf_max(slimb_t a, slimb_t b)
{
if (a > b)
return a;
else
return b;
}
static inline slimb_t bf_min(slimb_t a, slimb_t b)
{
if (a < b)
return a;
else
return b;
}
void bf_context_init(bf_context_t *s, bf_realloc_func_t *realloc_func,
void *realloc_opaque);
void bf_context_end(bf_context_t *s);
/* free memory allocated for the bf cache data */
void bf_clear_cache(bf_context_t *s);
static inline void *bf_realloc(bf_context_t *s, void *ptr, size_t size)
{
return s->realloc_func(s->realloc_opaque, ptr, size);
}
/* 'size' must be != 0 */
static inline void *bf_malloc(bf_context_t *s, size_t size)
{
return bf_realloc(s, NULL, size);
}
static inline void bf_free(bf_context_t *s, void *ptr)
{
/* must test ptr otherwise equivalent to malloc(0) */
if (ptr)
bf_realloc(s, ptr, 0);
}
void bf_init(bf_context_t *s, bf_t *r);
static inline void bf_delete(bf_t *r)
{
bf_context_t *s = r->ctx;
/* we accept to delete a zeroed bf_t structure */
if (s && r->tab) {
bf_realloc(s, r->tab, 0);
}
}
static inline void bf_neg(bf_t *r)
{
r->sign ^= 1;
}
static inline int bf_is_finite(const bf_t *a)
{
return (a->expn < BF_EXP_INF);
}
static inline int bf_is_nan(const bf_t *a)
{
return (a->expn == BF_EXP_NAN);
}
static inline int bf_is_zero(const bf_t *a)
{
return (a->expn == BF_EXP_ZERO);
}
static inline void bf_memcpy(bf_t *r, const bf_t *a)
{
*r = *a;
}
int bf_set_ui(bf_t *r, uint64_t a);
int bf_set_si(bf_t *r, int64_t a);
void bf_set_nan(bf_t *r);
void bf_set_zero(bf_t *r, int is_neg);
void bf_set_inf(bf_t *r, int is_neg);
int bf_set(bf_t *r, const bf_t *a);
void bf_move(bf_t *r, bf_t *a);
int bf_get_float64(const bf_t *a, double *pres, bf_rnd_t rnd_mode);
int bf_set_float64(bf_t *a, double d);
int bf_cmpu(const bf_t *a, const bf_t *b);
int bf_cmp_full(const bf_t *a, const bf_t *b);
int bf_cmp(const bf_t *a, const bf_t *b);
static inline int bf_cmp_eq(const bf_t *a, const bf_t *b)
{
return bf_cmp(a, b) == 0;
}
static inline int bf_cmp_le(const bf_t *a, const bf_t *b)
{
return bf_cmp(a, b) <= 0;
}
static inline int bf_cmp_lt(const bf_t *a, const bf_t *b)
{
return bf_cmp(a, b) < 0;
}
int bf_add(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags);
int bf_sub(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags);
int bf_add_si(bf_t *r, const bf_t *a, int64_t b1, limb_t prec, bf_flags_t flags);
int bf_mul(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags);
int bf_mul_ui(bf_t *r, const bf_t *a, uint64_t b1, limb_t prec, bf_flags_t flags);
int bf_mul_si(bf_t *r, const bf_t *a, int64_t b1, limb_t prec,
bf_flags_t flags);
int bf_mul_2exp(bf_t *r, slimb_t e, limb_t prec, bf_flags_t flags);
int bf_div(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec, bf_flags_t flags);
#define BF_DIVREM_EUCLIDIAN BF_RNDF
int bf_divrem(bf_t *q, bf_t *r, const bf_t *a, const bf_t *b,
limb_t prec, bf_flags_t flags, int rnd_mode);
int bf_rem(bf_t *r, const bf_t *a, const bf_t *b, limb_t prec,
bf_flags_t flags, int rnd_mode);
int bf_remquo(slimb_t *pq, bf_t *r, const bf_t *a, const bf_t *b, limb_t prec,
bf_flags_t flags, int rnd_mode);
/* round to integer with infinite precision */
int bf_rint(bf_t *r, int rnd_mode);
int bf_round(bf_t *r, limb_t prec, bf_flags_t flags);
int bf_sqrtrem(bf_t *r, bf_t *rem1, const bf_t *a);
int bf_sqrt(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
slimb_t bf_get_exp_min(const bf_t *a);
int bf_logic_or(bf_t *r, const bf_t *a, const bf_t *b);
int bf_logic_xor(bf_t *r, const bf_t *a, const bf_t *b);
int bf_logic_and(bf_t *r, const bf_t *a, const bf_t *b);
/* additional flags for bf_atof */
/* do not accept hex radix prefix (0x or 0X) if radix = 0 or radix = 16 */
#define BF_ATOF_NO_HEX (1 << 16)
/* accept binary (0b or 0B) or octal (0o or 0O) radix prefix if radix = 0 */
#define BF_ATOF_BIN_OCT (1 << 17)
/* Do not parse NaN or Inf */
#define BF_ATOF_NO_NAN_INF (1 << 18)
/* return the exponent separately */
#define BF_ATOF_EXPONENT (1 << 19)
int bf_atof(bf_t *a, const char *str, const char **pnext, int radix,
limb_t prec, bf_flags_t flags);
/* this version accepts prec = BF_PREC_INF and returns the radix
exponent */
int bf_atof2(bf_t *r, slimb_t *pexponent,
const char *str, const char **pnext, int radix,
limb_t prec, bf_flags_t flags);
int bf_mul_pow_radix(bf_t *r, const bf_t *T, limb_t radix,
slimb_t expn, limb_t prec, bf_flags_t flags);
/* Conversion of floating point number to string. Return a null
terminated string or NULL if memory error. *plen contains its
length if plen != NULL. The exponent letter is "e" for base 10,
"p" for bases 2, 8, 16 with a binary exponent and "@" for the other
bases. */
#define BF_FTOA_FORMAT_MASK (3 << 16)
/* fixed format: prec significant digits rounded with (flags &
BF_RND_MASK). Exponential notation is used if too many zeros are
needed.*/
#define BF_FTOA_FORMAT_FIXED (0 << 16)
/* fractional format: prec digits after the decimal point rounded with
(flags & BF_RND_MASK) */
#define BF_FTOA_FORMAT_FRAC (1 << 16)
/* free format:
For binary radices with bf_ftoa() and for bfdec_ftoa(): use the minimum
number of digits to represent 'a'. The precision and the rounding
mode are ignored.
For the non binary radices with bf_ftoa(): use as many digits as
necessary so that bf_atof() return the same number when using
precision 'prec', rounding to nearest and the subnormal
configuration of 'flags'. The result is meaningful only if 'a' is
already rounded to 'prec' bits. If the subnormal flag is set, the
exponent in 'flags' must also be set to the desired exponent range.
*/
#define BF_FTOA_FORMAT_FREE (2 << 16)
/* same as BF_FTOA_FORMAT_FREE but uses the minimum number of digits
(takes more computation time). Identical to BF_FTOA_FORMAT_FREE for
binary radices with bf_ftoa() and for bfdec_ftoa(). */
#define BF_FTOA_FORMAT_FREE_MIN (3 << 16)
/* force exponential notation for fixed or free format */
#define BF_FTOA_FORCE_EXP (1 << 20)
/* add 0x prefix for base 16, 0o prefix for base 8 or 0b prefix for
base 2 if non zero value */
#define BF_FTOA_ADD_PREFIX (1 << 21)
/* return "Infinity" instead of "Inf" and add a "+" for positive
exponents */
#define BF_FTOA_JS_QUIRKS (1 << 22)
char *bf_ftoa(size_t *plen, const bf_t *a, int radix, limb_t prec,
bf_flags_t flags);
/* modulo 2^n instead of saturation. NaN and infinity return 0 */
#define BF_GET_INT_MOD (1 << 0)
int bf_get_int32(int *pres, const bf_t *a, int flags);
int bf_get_int64(int64_t *pres, const bf_t *a, int flags);
/* the following functions are exported for testing only. */
void mp_print_str(const char *str, const limb_t *tab, limb_t n);
void bf_print_str(const char *str, const bf_t *a);
int bf_resize(bf_t *r, limb_t len);
int bf_get_fft_size(int *pdpl, int *pnb_mods, limb_t len);
int bf_normalize_and_round(bf_t *r, limb_t prec1, bf_flags_t flags);
int bf_can_round(const bf_t *a, slimb_t prec, bf_rnd_t rnd_mode, slimb_t k);
slimb_t bf_mul_log2_radix(slimb_t a1, unsigned int radix, int is_inv,
int is_ceil1);
int mp_mul(bf_context_t *s, limb_t *result,
const limb_t *op1, limb_t op1_size,
const limb_t *op2, limb_t op2_size);
limb_t mp_add(limb_t *res, const limb_t *op1, const limb_t *op2,
limb_t n, limb_t carry);
limb_t mp_add_ui(limb_t *tab, limb_t b, size_t n);
int mp_sqrtrem(bf_context_t *s, limb_t *tabs, limb_t *taba, limb_t n);
int mp_recip(bf_context_t *s, limb_t *tabr, const limb_t *taba, limb_t n);
limb_t bf_isqrt(limb_t a);
/* transcendental functions */
int bf_const_log2(bf_t *T, limb_t prec, bf_flags_t flags);
int bf_const_pi(bf_t *T, limb_t prec, bf_flags_t flags);
int bf_exp(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_log(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
#define BF_POW_JS_QUIRKS (1 << 16) /* (+/-1)^(+/-Inf) = NaN, 1^NaN = NaN */
int bf_pow(bf_t *r, const bf_t *x, const bf_t *y, limb_t prec, bf_flags_t flags);
int bf_cos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_sin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_tan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_atan(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_atan2(bf_t *r, const bf_t *y, const bf_t *x,
limb_t prec, bf_flags_t flags);
int bf_asin(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
int bf_acos(bf_t *r, const bf_t *a, limb_t prec, bf_flags_t flags);
/* decimal floating point */
static inline void bfdec_init(bf_context_t *s, bfdec_t *r)
{
bf_init(s, (bf_t *)r);
}
static inline void bfdec_delete(bfdec_t *r)
{
bf_delete((bf_t *)r);
}
static inline void bfdec_neg(bfdec_t *r)
{
r->sign ^= 1;
}
static inline int bfdec_is_finite(const bfdec_t *a)
{
return (a->expn < BF_EXP_INF);
}
static inline int bfdec_is_nan(const bfdec_t *a)
{
return (a->expn == BF_EXP_NAN);
}
static inline int bfdec_is_zero(const bfdec_t *a)
{
return (a->expn == BF_EXP_ZERO);
}
static inline void bfdec_memcpy(bfdec_t *r, const bfdec_t *a)
{
bf_memcpy((bf_t *)r, (const bf_t *)a);
}
int bfdec_set_ui(bfdec_t *r, uint64_t a);
int bfdec_set_si(bfdec_t *r, int64_t a);
static inline void bfdec_set_nan(bfdec_t *r)
{
bf_set_nan((bf_t *)r);
}
static inline void bfdec_set_zero(bfdec_t *r, int is_neg)
{
bf_set_zero((bf_t *)r, is_neg);
}
static inline void bfdec_set_inf(bfdec_t *r, int is_neg)
{
bf_set_inf((bf_t *)r, is_neg);
}
static inline int bfdec_set(bfdec_t *r, const bfdec_t *a)
{
return bf_set((bf_t *)r, (bf_t *)a);
}
static inline void bfdec_move(bfdec_t *r, bfdec_t *a)
{
bf_move((bf_t *)r, (bf_t *)a);
}
static inline int bfdec_cmpu(const bfdec_t *a, const bfdec_t *b)
{
return bf_cmpu((const bf_t *)a, (const bf_t *)b);
}
static inline int bfdec_cmp_full(const bfdec_t *a, const bfdec_t *b)
{
return bf_cmp_full((const bf_t *)a, (const bf_t *)b);
}
static inline int bfdec_cmp(const bfdec_t *a, const bfdec_t *b)
{
return bf_cmp((const bf_t *)a, (const bf_t *)b);
}
static inline int bfdec_cmp_eq(const bfdec_t *a, const bfdec_t *b)
{
return bfdec_cmp(a, b) == 0;
}
static inline int bfdec_cmp_le(const bfdec_t *a, const bfdec_t *b)
{
return bfdec_cmp(a, b) <= 0;
}
static inline int bfdec_cmp_lt(const bfdec_t *a, const bfdec_t *b)
{
return bfdec_cmp(a, b) < 0;
}
int bfdec_add(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec,
bf_flags_t flags);
int bfdec_sub(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec,
bf_flags_t flags);
int bfdec_add_si(bfdec_t *r, const bfdec_t *a, int64_t b1, limb_t prec,
bf_flags_t flags);
int bfdec_mul(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec,
bf_flags_t flags);
int bfdec_mul_si(bfdec_t *r, const bfdec_t *a, int64_t b1, limb_t prec,
bf_flags_t flags);
int bfdec_div(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec,
bf_flags_t flags);
int bfdec_divrem(bfdec_t *q, bfdec_t *r, const bfdec_t *a, const bfdec_t *b,
limb_t prec, bf_flags_t flags, int rnd_mode);
int bfdec_rem(bfdec_t *r, const bfdec_t *a, const bfdec_t *b, limb_t prec,
bf_flags_t flags, int rnd_mode);
int bfdec_rint(bfdec_t *r, int rnd_mode);
int bfdec_sqrt(bfdec_t *r, const bfdec_t *a, limb_t prec, bf_flags_t flags);
int bfdec_round(bfdec_t *r, limb_t prec, bf_flags_t flags);
int bfdec_get_int32(int *pres, const bfdec_t *a);
int bfdec_pow_ui(bfdec_t *r, const bfdec_t *a, limb_t b);
char *bfdec_ftoa(size_t *plen, const bfdec_t *a, limb_t prec, bf_flags_t flags);
int bfdec_atof(bfdec_t *r, const char *str, const char **pnext,
limb_t prec, bf_flags_t flags);
/* the following functions are exported for testing only. */
extern const limb_t mp_pow_dec[LIMB_DIGITS + 1];
void bfdec_print_str(const char *str, const bfdec_t *a);
static inline int bfdec_resize(bfdec_t *r, limb_t len)
{
return bf_resize((bf_t *)r, len);
}
int bfdec_normalize_and_round(bfdec_t *r, limb_t prec1, bf_flags_t flags);
#endif /* LIBBF_H */

161
third_party/libbf/readme.txt vendored Normal file
View file

@ -0,0 +1,161 @@
Tiny Big Float library
----------------------
Copyright (c) 2017-2020 Fabrice Bellard
LibBF is a small library to handle arbitrary precision binary or
decimal floating point numbers. Its compiled size is about 90 KB of
x86 code and has no dependency on other libraries. It is not the
fastest library nor the smallest but it tries to be simple while using
asymptotically optimal algorithms. The basic arithmetic operations
have a near linear running time.
The TinyPI example computes billions of digits of Pi using the
Chudnovsky formula.
1) Features
-----------
- Arbitrary precision floating point numbers in base 2 using the IEEE
754 semantics (including subnormal numbers, infinities and
NaN).
- All operations are exactly rounded using the 5 IEEE 754 rounding
modes (round to nearest with ties to even or away from zero, round
to zero, -/+ infinity). The additional non-deterministic faithful
rounding mode is supported when a lower or deterministic running
time is necessary.
- Stateless API (each function takes as input the rounding mode,
mantissa and exponent precisions in bits and return the IEEE status
flags).
- The basic arithmetic operations (addition, subtraction,
multiplication, division, square root) have a near linear running
time.
- Multiplication using a SIMD optimized Number Theoretic Transform.
- Exactly rounded floating point input and output in any base between
2 and 36 with near linear runnning time. Floating point output can
select the smallest amount of digits to get the required precision.
- Transcendental functions are supported (exp, log, pow, sin, cos, tan,
asin, acos, atan, atan2).
- Operations on arbitrarily large integers are supported by using a
special "infinite" precision. Integer division with remainder and
logical operations (assuming two complement binary representation)
are implemented.
- Arbitrary precision floating point numbers in base 10 corresponding
to the IEEE 754 2008 semantics with the limitation that the mantissa
is always normalized. The basic arithmetic operations, output and
input are supported with a quadratic running time.
- Easy to embed: a few C files need to be copied, the memory allocator
can be redefined, the memory allocation failures are tested.
- MIT license.
2) Compilation
--------------
Edit the top of the Makefile to select the build options. By default,
the MPFR library is used to compile the test tools (bftest and
bfbench) but it is not needed to build libbf. The included SoftFP code
(softfp* files) is only used by the bftest test tool.
TinyPI example: the "tinypi" executable uses the portable code. The
"tinypi-avx2" executable uses the AVX2 implementation. An x86 CPU of
at least the Intel Haswell generation is necessary for AVX2.
3) Design principles
--------------------
- Base 2 and IEEE 754 semantics were chosen so that it is possible to
get good performance and to compare the results with other libraries
or hardware implementations. Moreover, base 2 arbitrary precision is
easier to analyse and implement.
- The support of subnormal numbers and of a configurable number of
bits for the exponent allows the exact emulation of IEEE 754
floating hardware.
- The stateless API ensures that there is no global state to save and
restore between operations. The rounding mode, subnormal flag and
number of exponent bits are ored to a single "flags" parameter to
limit the verbosity of the API. The number of exponent bits 'n' is
specified as '(M-n)' where M is the maximum number of exponent bits
so that '0' always indicates the maximum number of exponent bits.
- All the IEEE 754 status flags are returned by each operation. The
user can easily or them when necessary.
- Unlike other libraries (such as MPFR [2]), the numbers have no
attached precision. The general rule is that each operation is
internally computed with infinite precision and then rounded with
the precision and rounding mode specified for the operation.
- In many computations it is necessary to use arbitrarily large
integers. LibBF support them without adding another number type by
providing a special "infinite" precision. There is a small overhead
of course because they are manipulated as floating point numbers but
there is no cost to convert between floating point numbers and
integers.
- The faithful rounding mode (i.e. the result is rounded to - or
+infinity non deterministically) is supported for all operations. It
usually gives a faster and deterministic running time. The
transcendental functions, inverse or inverse square root are
internally implemented to give a faithful rounding. When a
non-faithful rounding is requested by the user, the Ziv rounding
algorithm is invoked.
4) Implementation notes
-----------------------
- The code was tested on a 64 bit x86 CPU. It should be portable to
other CPUs. The portable version handles numbers with up to 4*10^16
digits. The AVX2 version handles numbers with up to 8*10^12 digits.
- 32 bits: the code compiles on 32 bit architectures but it is not
designed to be efficient nor scalable in this case. The size of the
numbers is limited to about 10 million digits.
- The Number Theoretic Transform is not the fastest algorithm for
small to medium numbers (i.e. a few million digits), but it gets
better when the size of the numbers grows. There is no round-off
errors as with Fast Fourier Transform, the memory usage is much
smaller and it is potentially easier to parallelize. This code
contains an original SIMD (AVX2 on x86) implementation using 64 bit
floating point numbers. It relies on the fact that the fused
multiply accumulate (FMA) operation gives access to the full
precision of the product of two 64 bit floating point numbers. The
portable code relies on the fact that the C compiler supports a
double word integer type (i.e. 128 bit integers on 64 bit). The
modulo operations were replaced with multiplications which are
usually faster.
- Base conversion: the algorithm is not the fastest one but it is
simple and still gives a near linear running time.
- This library reuses some ideas from TachusPI (
http://bellard.org/pi/pi2700e9/tpi.html ) . It is about 4 times
slower to compute Pi but is much smaller and simpler.
5) Known limitations
--------------------
- In some operations (such as the transcendental ones), there is no
rigourous proof of the rounding error. We expect to improve it by
reusing ideas from the MPFR algorithms. Some unlikely
overflow/underflow cases are also not handled in exp or pow.
- The transcendental operations are not speed optimized and do not use
an asymptotically optimal algorithm (the running time is in
O(n^(1/2)*M(n)) where M(n) is the time to multiply two n bit
numbers). A possible solution would be to implement a binary
splitting algorithm for exp and sin/cos (see [1]) and to use a
Newton based inversion to get log and atan.
- Memory allocation errors are not always correctly reported for the
transcendental operations.
6) References
-------------
[1] Modern Computer Arithmetic, Richard Brent and Paul Zimmermann,
Cambridge University Press, 2010
(https://members.loria.fr/PZimmermann/mca/pub226.html).
[2] The GNU MPFR Library (http://www.mpfr.org/)

View file

@ -35,7 +35,7 @@ TERM_EMPTY = '\033[0m'
clang_format_exts = ['.cpp', '.h']
skip_dirs = ['build', 'CMakeFiles', 'docs', 'out', 'test', 'tools', 'double_conversion', 'GCutil', 'lz4', 'rapidjson', 'windows', '.git']
skip_dirs = ['build', 'CMakeFiles', 'docs', 'out', 'test', 'tools', 'double_conversion', 'GCutil', 'lz4', 'rapidjson', 'windows', '.git', 'libbf']
skip_files = []