/* SIPE - Small Integer Plus Exponent * * Copyright (c) 2008-2009 Vincent Lefevre . * Arenaire project, INRIA / LIP / ENS-Lyon, France. * * Note: GCC extensions and implementation-defined behaviors are used! * Exponents are assumed to be small enough so that no integer overflow * can occur. Code is optimized for precisions (prec) known at compile * time and for non-zero numbers. Non-zero numbers must be normalized. * Normalization can be done by using SIPE_ROUND. For sipe_nextabove * and sipe_nextbelow, 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 33176 2009-11-17 11:50:19Z vinc17/vin $" #ifndef __GNUC__ #error "GCC is required" #endif #include #include #include #include #ifndef SIPE_MINSIZE #define SIPE_MINSIZE 32 #endif #if SIPE_MINSIZE != 32 && SIPE_MINSIZE != 64 #error "SIPE_MINSIZE must be 32 or 64" #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; 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) /* 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 */ typedef struct { sipe_int_t i; sipe_exp_t e; } sipe_t; #define SIPE_ROUND(X,PREC) \ do \ if (SIPE_LIKELY (X.i != 0)) \ { \ sipe_int_t _i; \ int _s; \ \ _i = X.i < 0 ? - X.i : X.i; \ _s = (PREC) - SIPE_SIZE + sipe_clz (_i); \ if (_s > 0) \ { \ X.i <<= _s; \ X.e -= _s; \ } \ else if (_s < 0) \ { \ sipe_int_t _j; \ int _ns; \ \ _ns = - 1 - _s; \ _j = _i >> _ns; \ if ((_j & 2) | (_i - (_j << _ns))) \ _j++; \ _j >>= 1; \ if (SIPE_UNLIKELY (_j == SIPE_ONE << (PREC))) \ { \ _j >>= 1; \ _ns++; \ } \ X.i = X.i < 0 ? - _j : _j; \ X.e += _ns + 1; \ } \ } \ while (0); /* Convert a SIPE value x, which must be an integer representable in a sipe_int_t, into an integer. Note: >> on a negative integer is implementation-defined, but we depend on GCC, and GCC documents >> as acting on negative numbers by sign extension. */ 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 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; } static inline sipe_t sipe_add_si (sipe_t x, sipe_int_t y, int prec) { sipe_t r = { y, 0 }; SIPE_ROUND (r, prec); return sipe_add (x, r, prec); } static inline sipe_t sipe_sub_si (sipe_t x, sipe_int_t y, int prec) { sipe_t r = { y, 0 }; SIPE_ROUND (r, prec); return sipe_sub (x, r, prec); } #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_ONE << (prec - 1)))) \ { \ r.i = OPN ((SIPE_ONE << 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_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; } static inline sipe_t sipe_mul_si (sipe_t x, sipe_int_t y, int prec) { sipe_t r = { y, 0 }; SIPE_ROUND (r, prec); return sipe_mul (x, r, prec); } #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; \ delta = r.e - z.e; \ if (SIPE_UNLIKELY (z.i == 0)) \ { \ SIPE_ROUND (r, prec); \ return r; \ } \ if (delta > 2 * prec + 1) \ { \ /* Warning! The sign of z.i is important if r is the \ middle of two consecutive machine numbers. */ \ r.i = 2 * r.i + (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_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) ? (G) : (L); \ } \ 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; } 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_ONE << (prec - 2); bit != 0; bit >>= 1) putc (x.i & bit ? '1' : '0', stream); fprintf (stream, "e%" PRIdFAST8, x.e + prec - 1); }