/* SIPE - Small Integer Plus Exponent * * Copyright (c) 2008-2013 Vincent Lefevre . * Arenaire/AriC projects, INRIA / LIP / ENS-Lyon, France. * * Note: In non-SIPE_FLOAT mode, GCC (or compatible) is required due to * the use of GCC extensions and implementation-defined behaviors. * Exponents are assumed to be small enough so that no integer overflow * or floating-point underflow/overflow can occur. Code is optimized for * precisions (prec) known at compile time and for non-zero numbers. * In non-SIPE_FLOAT mode, non-zero numbers must be normalized. * Normalization can be done by using SIPE_ROUND. For sipe_nextabove, * sipe_nextbelow, and sipe_nexttozero, the argument must not be 0. * Signed zeros are not supported. * * This code is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2.1 of * the License, or (at your option) any later version. * * This code is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this file; if not, write to the Free Software Foundation, * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ #define SIPE_SVNID "$Id: sipe.h 60133 2013-05-03 12:54:10Z vinc17/ypig $" #ifndef SIPE_MINSIZE #define SIPE_MINSIZE 32 #endif #if SIPE_MINSIZE != 32 && SIPE_MINSIZE != 64 #error "SIPE_MINSIZE must be 32 or 64" #endif #include #include #include #include #if SIPE_CHECK #include #define SIPE_ASSERT(EXPR) assert(EXPR) #else #define SIPE_ASSERT(EXPR) ((void) 0) #endif #define SIPE_INT_TYPE_AUX(S) int_fast ## S ## _t #define SIPE_INT_TYPE(S) SIPE_INT_TYPE_AUX(S) typedef SIPE_INT_TYPE(SIPE_MINSIZE) sipe_int_t; /* WARNING! SIPE_FLOAT mode is new (written on 2013-04-07 in a plane) and * it isn't well tested yet. It also has currently no fma/fms support. */ #ifdef SIPE_FLOAT typedef int sipe_exp_t; #define SIPE_CONCAT2_AUX(A,B) A##B #define SIPE_CONCAT2(A,B) SIPE_CONCAT2_AUX(A,B) #define SIPE_CONCAT3_AUX(A,B,C) A##B##C #define SIPE_CONCAT3(A,B,C) SIPE_CONCAT3_AUX(A,B,C) #include #include #if SIPE_FLOAT == 1 typedef float sipe_t; #define SIPE_NATIVE_PREC FLT_MANT_DIG #define SIPE_SFX f #elif SIPE_FLOAT == 2 typedef double sipe_t; #define SIPE_NATIVE_PREC DBL_MANT_DIG #define SIPE_SFX #elif SIPE_FLOAT == 3 typedef long double sipe_t; #define SIPE_NATIVE_PREC LDBL_MANT_DIG #define SIPE_SFX l #endif #define SIPE_PREC_MAX ((SIPE_NATIVE_PREC-1)/2) #define SIPE_POWD(PREC) (1UL << (SIPE_NATIVE_PREC - PREC)) static inline sipe_int_t sipe_intm (sipe_t x, int prec) { int e; return SIPE_CONCAT2(frexp,SIPE_SFX) (x, &e) * (1UL << prec); } static inline sipe_exp_t sipe_qexp (sipe_t x, int prec) { int e; SIPE_CONCAT2(frexp,SIPE_SFX) (x, &e); return e - prec; } #define SIPE_ROUND_ERR(X,C,ERR,PREC) \ do \ { \ sipe_t _y, _z; \ _y = (X) * (SIPE_POWD(PREC) + 1); \ _z = (X) - _y; \ if (C) \ { \ sipe_t _r = _z + _y; \ ERR ((X) - _r); \ (X) = _r; \ } \ else \ (X) = _z + _y; \ } \ while (0) #define SIPE_ROUND(X,PREC) SIPE_ROUND_ERR(X,0,(void),PREC) static inline sipe_int_t sipe_to_int (sipe_t x) { return (sipe_int_t) x; } #define SIPE_DEFOP(OP,OPS) \ static inline sipe_t sipe_##OP (sipe_t x, sipe_t y, int prec) \ { \ sipe_t r = x OPS y; \ SIPE_ROUND (r, prec); \ return r; \ } SIPE_DEFOP(add,+) SIPE_DEFOP(sub,-) SIPE_DEFOP(mul,*) static inline sipe_t sipe_neg (sipe_t x, int prec) { return -x; } static inline sipe_t sipe_set_si (sipe_int_t x, int prec) { return (sipe_t) x; } #define SIPE_DEFNEXT(DIR,VAL) \ static inline sipe_t sipe_next##DIR (sipe_t x, int prec) \ { \ return x + \ (SIPE_CONCAT2(nextafter,SIPE_SFX) (x, (sipe_t) VAL) - x) \ * SIPE_POWD(prec); \ } SIPE_DEFNEXT(above,+INFINITY) SIPE_DEFNEXT(below,-INFINITY) SIPE_DEFNEXT(tozero,0) static inline int sipe_zero_p (sipe_t x) { return x == 0; } #define SIPE_DEFCMP(OP,OPS) \ static inline int sipe_##OP (sipe_t x, sipe_t y) \ { \ return x OPS y; \ } SIPE_DEFCMP(eq,==) SIPE_DEFCMP(ne,!=) SIPE_DEFCMP(le,<=) SIPE_DEFCMP(lt,<) SIPE_DEFCMP(ge,>=) SIPE_DEFCMP(gt,>) #define SIPE_DEFMINMAX(FCT) \ static inline sipe_t sipe_##FCT (sipe_t x, sipe_t y) \ { \ return SIPE_CONCAT3(f,FCT,SIPE_SFX) (x, y); \ } SIPE_DEFMINMAX(min) SIPE_DEFMINMAX(max) #define SIPE_DEFCMPMAG(OP,TYPE,X,Y,R) \ static inline TYPE sipe_##OP##mag (sipe_t x, sipe_t y) \ { \ sipe_t absx, absy; \ absx = SIPE_CONCAT2(fabs,SIPE_SFX) (x); \ absy = SIPE_CONCAT2(fabs,SIPE_SFX) (y); \ if (absx < absy) return X; \ if (absx > absy) return Y; \ return R; \ } SIPE_DEFCMPMAG(min,sipe_t,x,y,x<0?x:y) SIPE_DEFCMPMAG(max,sipe_t,y,x,x<0?y:x) SIPE_DEFCMPMAG(cmp,int,-1,1,0) #define SIPE_DEFOPSI(OP) \ static inline sipe_t sipe_##OP##_si (sipe_t x, sipe_int_t y, \ int prec) \ { \ return sipe_##OP (x, sipe_set_si (y, prec), prec); \ } SIPE_DEFOPSI(add) SIPE_DEFOPSI(sub) SIPE_DEFOPSI(mul) static inline void sipe_outbin (FILE *stream, sipe_t x, int prec) { int e; if (x == 0) { putc ('0', stream); return; } else if (x < 0) { putc ('-', stream); x = - x; } fputs ("1.", stream); x = SIPE_CONCAT2(frexp,SIPE_SFX) (x, &e) * 2 - 1; while (--prec) { x *= 2; putc ((int) x ? '1' : '0', stream); x -= (int) x; } fprintf (stream, "e%d", e - 1); SIPE_ASSERT (x == 0); } /* Multiplication: rounded result and error. */ #define SIPE_2MUL(R,S,X,Y,PREC) \ do \ { \ (R) = (X) * (Y); \ SIPE_ROUND_ERR (R, 1, (S) =, PREC); \ } \ while (0) #else /* SIPE_FLOAT */ #ifndef __GNUC__ #error "GCC is required" #endif typedef int_fast8_t sipe_exp_t; #define SIPE_TYPESIZE(T) (sizeof(T) * CHAR_BIT) #define SIPE_SIZE SIPE_TYPESIZE(sipe_int_t) #define SIPE_ONE ((sipe_int_t) 1) #define SIPE_TWO_TO(N) (SIPE_ONE << (N)) /* To avoid integer overflows in the code below, one must have: max_delta + prec <= SIPE_SIZE - 1 (sign bit), where max_delta is the maximum value of delta that does not lead to early rounding. */ #if SIPE_USE_FMA #define SIPE_PREC_MAX ((SIPE_SIZE - 2) / 3) #else #define SIPE_PREC_MAX ((SIPE_SIZE - 2) / 2) #endif #define SIPE_LIKELY(C) (__builtin_expect (!!(C), 1)) #define SIPE_UNLIKELY(C) (__builtin_expect (!!(C), 0)) #define SIPE_ABSINT(X) ((X) >= 0 ? (X) : - (X)) /* Count the leading zeros of i > 0 */ static inline int sipe_clz (sipe_int_t i) { if (SIPE_TYPESIZE (int) >= SIPE_SIZE) return __builtin_clz ((unsigned int) i) - (SIPE_TYPESIZE (int) - SIPE_SIZE); if (SIPE_TYPESIZE (long) >= SIPE_SIZE) return __builtin_clzl ((unsigned long) i) - (SIPE_TYPESIZE (long) - SIPE_SIZE); if (SIPE_TYPESIZE (long long) >= SIPE_SIZE) return __builtin_clzll ((unsigned long long) i) - (SIPE_TYPESIZE (long long) - SIPE_SIZE); fprintf (stderr, "sipe: unsupported sipe_int_t size for sipe_clz"); exit (119); } /* Number X is X.i * 2^(X.e), where X.i = 0 or 2^(p-1) <= |X.i| < 2^p; if X.i = 0, then X.e should be set to 0 in order to avoid undefined behavior in some cases, such as squaring 0 multiple times (this has no impact on the performance in practice). */ typedef struct { sipe_int_t i; sipe_exp_t e; } sipe_t; static inline sipe_int_t sipe_intm (sipe_t x, int prec) { return x.i; } static inline sipe_exp_t sipe_qexp (sipe_t x, int prec) { return x.e; } /* From some tests done with GCC on Intel Xeon and AMD Opteron, SIPE_ROUND_ZOPT = 0 yields faster code if profiling is used. Otherwise SIPE_ROUND_ZOPT = 1 is faster on AMD Opteron. */ #ifndef SIPE_ROUND_ZOPT #define SIPE_ROUND_ZOPT 0 #endif /* Auxiliary function for SIPE_ROUND and SIPE_2MUL */ #define SIPE_ROUND_ERR(X,C,ERR,PREC) \ do \ if (SIPE_LIKELY ((X).i != 0)) \ { \ sipe_int_t _i, _j; \ int _s, _ns; \ \ _i = SIPE_ABSINT ((X).i); \ _s = (PREC) - SIPE_SIZE + sipe_clz (_i); \ if ((SIPE_ROUND_ZOPT) ? _s >= 0 : _s > 0) \ { \ (X).i <<= _s; \ (X).e -= _s; \ ERR (sipe_t) { 0, 0 }; \ } \ else if (!(SIPE_ROUND_ZOPT) && _s == 0) \ { \ ERR (sipe_t) { 0, 0 }; \ } \ else \ { \ _ns = - 1 - _s; \ _j = _i >> _ns; \ if ((_j & 2) | (_i - (_j << _ns))) \ _j++; \ _j >>= 1; \ if (SIPE_UNLIKELY (_j == SIPE_TWO_TO (PREC))) \ { \ _j >>= 1; \ _ns++; \ } \ if (C) \ { \ _i -= _j << (_ns + 1); \ if (_i == 0) \ ERR (sipe_t) { 0, 0 }; \ else \ { \ _s = (PREC) - SIPE_SIZE + \ sipe_clz (SIPE_ABSINT (_i)); \ SIPE_ASSERT (_s >= 0); \ ERR (sipe_t) \ { ((X).i >= 0 ? _i : - _i) << _s, (X).e - _s }; \ } \ } \ (X).i = (X).i >= 0 ? _j : - _j; \ (X).e += _ns + 1; \ } \ } \ else \ { \ (X).e = 0; \ ERR (sipe_t) { 0, 0 }; \ } \ while (0) /* Round and normalize a variable X. */ #define SIPE_ROUND(X,PREC) SIPE_ROUND_ERR(X,0,(void),PREC) /* Convert a SIPE value x, which must be an integer representable in a sipe_int_t, into an integer, i.e. the conversion must be exact. Note: >> on a negative integer is implementation-defined, but we depend on GCC, and GCC documents >> as acting on negative numbers by sign extension (it behaves as wanted because x.i is assumed to be a multiple of -x.e). We also need to make sure that the shift count is less than the width of sipe_int_t (required by ISO C99 6.5.7#3). If x.i == 0, then x.e == 0 and x.i is returned directly. If x.i != 0, then: * if x.e < 0, then |x.i| / 2^|x.e| must be an integer different from 0, thus >= 1, so that |x.e| is necessarily less than the width of sipe_int_t; * if x.e > 0, then |x.i| * 2^|x.e| >= 2^|x.e|, which must be representable in a sipe_int_t, so that |x.e| is necessarily less than the width of sipe_int_t. */ static inline sipe_int_t sipe_to_int (sipe_t x) { return x.e < 0 ? x.i >> - x.e : x.e > 0 ? x.i << x.e : x.i; } #ifdef __GNU_MP_VERSION /* Convert a SIPE value x (not necessarily an integer) into an mpz_t integer, rounding the value toward zero (truncation). |x.e| must fit in an mp_bitcnt_t (unsigned long at least up to GMP 5), which is not a problem in practice. We could consider the case x.i == 0 separately, but don't need to. */ static inline void sipe_to_mpz (mpz_t z, sipe_t x) { mpz_set_si (z, x.i); if (x.e < 0) mpz_tdiv_q_2exp (z, z, - x.e); else mpz_mul_2exp (z, z, x.e); } #endif #define SIPE_DEFADDSUB(OP,ADD,OPS) \ static inline sipe_t sipe_##OP (sipe_t x, sipe_t y, int prec) \ { \ sipe_exp_t delta = x.e - y.e; \ sipe_t r; \ \ if (SIPE_UNLIKELY (x.i == 0)) \ return (ADD) ? y : (sipe_t) { - y.i, y.e }; \ if (SIPE_UNLIKELY (y.i == 0) || delta > prec + 1) \ return x; \ if (delta < - (prec + 1)) \ return (ADD) ? y : (sipe_t) { - y.i, y.e }; \ r = delta < 0 ? \ ((sipe_t) { (x.i) OPS (y.i << - delta), x.e }) : \ ((sipe_t) { (x.i << delta) OPS (y.i), y.e }); \ SIPE_ROUND (r, prec); \ return r; \ } SIPE_DEFADDSUB(add,1,+) SIPE_DEFADDSUB(sub,0,-) static inline sipe_t sipe_neg (sipe_t x, int prec) { return (sipe_t) { - x.i, x.e }; } static inline sipe_t sipe_set_si (sipe_int_t x, int prec) { sipe_t r = { x, 0 }; SIPE_ROUND (r, prec); return r; } #define SIPE_DEFNEXT(DIR,OPS,OPN) \ static inline sipe_t sipe_next##DIR (sipe_t x, int prec) \ { \ sipe_t r; \ if (SIPE_UNLIKELY (x.i == OPN (SIPE_TWO_TO (prec - 1)))) \ { \ r.i = OPN (SIPE_TWO_TO (prec) - 1); \ r.e = x.e - 1; \ } \ else \ { \ r.i = x.i OPS 1; \ r.e = x.e; \ SIPE_ROUND (r, prec); \ } \ return r; \ } SIPE_DEFNEXT(above,+,-) SIPE_DEFNEXT(below,-,+) static inline sipe_t sipe_nexttozero (sipe_t x, int prec) { if (x.i < 0) return sipe_nextabove (x, prec); else return sipe_nextbelow (x, prec); } static inline sipe_t sipe_mul (sipe_t x, sipe_t y, int prec) { sipe_t r; r.i = x.i * y.i; r.e = x.e + y.e; SIPE_ROUND (r, prec); return r; } #if SIPE_USE_FMA #define SIPE_DEFFMAFMS(OP,FMA,OPS) \ static inline sipe_t sipe_##OP (sipe_t x, sipe_t y, sipe_t z, \ int prec) \ { \ sipe_t r; \ sipe_exp_t delta; \ \ r.i = x.i * y.i; \ if (SIPE_UNLIKELY (r.i == 0)) \ return (FMA) ? z : (sipe_t) { - z.i, z.e }; \ r.e = x.e + y.e; \ if (SIPE_UNLIKELY (z.i == 0)) \ { \ SIPE_ROUND (r, prec); \ return r; \ } \ delta = r.e - z.e; \ if (delta > prec) \ { \ /* Warning! The sign of z.i is important if r is the \ middle of two consecutive machine numbers. */ \ r.i = 2 * r.i OPS (z.i < 0 ? -1 : 1); \ r.e--; \ SIPE_ROUND (r, prec); \ return r; \ } \ if (delta < - (2 * prec + 1)) \ return (FMA) ? z : (sipe_t) { - z.i, z.e }; \ r = delta < 0 ? \ ((sipe_t) { (r.i) OPS (z.i << - delta), r.e }) : \ ((sipe_t) { (r.i << delta) OPS (z.i), z.e }); \ SIPE_ROUND (r, prec); \ return r; \ } SIPE_DEFFMAFMS(fma,1,+) SIPE_DEFFMAFMS(fms,0,-) #endif static inline int sipe_zero_p (sipe_t x) { return x.i == 0; } static inline int sipe_eq (sipe_t x, sipe_t y) { return x.i == y.i && x.e == y.e; } static inline int sipe_ne (sipe_t x, sipe_t y) { return x.i != y.i || x.e != y.e; } #define SIPE_DEFCMP(OP,TYPE,E,L,G) \ static inline TYPE sipe_##OP##pos (sipe_t x, sipe_t y) \ { \ if (x.e < y.e) return (L); \ if (x.e > y.e) return (G); \ return ((E) ? x.i < y.i : x.i <= y.i) ? (L) : (G); \ } \ static inline TYPE sipe_##OP##neg (sipe_t x, sipe_t y) \ { \ if (x.e < y.e) return (G); \ if (x.e > y.e) return (L); \ return ((E) ? x.i < y.i : x.i <= y.i) ? (L) : (G); \ } \ static inline TYPE sipe_##OP (sipe_t x, sipe_t y) \ { \ if ((E) && x.i == 0 && y.i == 0) return (G); \ if (x.i <= 0 && y.i >= 0) return (L); \ if (x.i >= 0 && y.i <= 0) return (G); \ return x.i < 0 ? \ sipe_##OP##neg (x, y) : \ sipe_##OP##pos (x, y); \ } SIPE_DEFCMP(le,int,0,1,0) SIPE_DEFCMP(lt,int,1,1,0) SIPE_DEFCMP(ge,int,1,0,1) SIPE_DEFCMP(gt,int,0,0,1) SIPE_DEFCMP(min,sipe_t,0,x,y) SIPE_DEFCMP(max,sipe_t,0,y,x) #define SIPE_DEFCMPMAG(OP,X,Y) \ static inline sipe_t sipe_##OP##mag (sipe_t x, sipe_t y) \ { \ sipe_int_t absxi, absyi; \ if (x.i == 0) return X; \ if (y.i == 0) return Y; \ if (x.e < y.e) return X; \ if (x.e > y.e) return Y; \ absxi = SIPE_ABSINT (x.i); \ absyi = SIPE_ABSINT (y.i); \ if (absxi < absyi) return X; \ if (absxi > absyi) return Y; \ return x.i < 0 ? X : Y; \ } SIPE_DEFCMPMAG(min,x,y) SIPE_DEFCMPMAG(max,y,x) static inline int sipe_cmpmag (sipe_t x, sipe_t y) { sipe_int_t absxi, absyi; if (x.i == 0 && y.i == 0) return 0; if (x.i == 0) return -1; if (y.i == 0) return 1; if (x.e < y.e) return -1; if (x.e > y.e) return 1; absxi = SIPE_ABSINT (x.i); absyi = SIPE_ABSINT (y.i); if (absxi < absyi) return -1; if (absxi > absyi) return 1; return 0; } #define SIPE_DEFOPSI(OP) \ static inline sipe_t sipe_##OP##_si (sipe_t x, sipe_int_t y, \ int prec) \ { \ return sipe_##OP (x, sipe_set_si (y, prec), prec); \ } SIPE_DEFOPSI(add) SIPE_DEFOPSI(sub) SIPE_DEFOPSI(mul) static inline void sipe_outbin (FILE *stream, sipe_t x, int prec) { sipe_int_t bit; if (x.i == 0) { putc ('0', stream); return; } else if (x.i < 0) { putc ('-', stream); x.i = - x.i; } fputs ("1.", stream); for (bit = SIPE_TWO_TO (prec - 2); bit != 0; bit >>= 1) putc (x.i & bit ? '1' : '0', stream); fprintf (stream, "e%" PRIdFAST8, x.e + prec - 1); } /* Multiplication: rounded result and error. */ #define SIPE_2MUL(R,S,X,Y,PREC) \ do \ { \ (R) = (sipe_t) { (X).i * (Y).i, (X).e + (Y).e }; \ SIPE_ROUND_ERR (R, 1, (S) =, PREC); \ } \ while (0) #endif /* ========================== End of sipe.h ========================== */