1 /* Compute a product of X, X+1, ..., with an error estimate. 2 Copyright (C) 2013-2018 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 5 The GNU C Library is free software; you can redistribute it and/or 6 modify it under the terms of the GNU Lesser General Public 7 License as published by the Free Software Foundation; either 8 version 2.1 of the License, or (at your option) any later version. 9 10 The GNU C Library is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13 Lesser General Public License for more details. 14 15 You should have received a copy of the GNU Lesser General Public 16 License along with the GNU C Library; if not, see 17 <http://www.gnu.org/licenses/>. */ 18 19 #include "quadmath-imp.h" 20 21 /* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N 22 - 1, in the form R * (1 + *EPS) where the return value R is an 23 approximation to the product and *EPS is set to indicate the 24 approximate error in the return value. X is such that all the 25 values X + 1, ..., X + N - 1 are exactly representable, and X_EPS / 26 X is small enough that factors quadratic in it can be 27 neglected. */ 28 29 __float128 30 __quadmath_gamma_productq (__float128 x, __float128 x_eps, int n, __float128 *eps) 31 { 32 SET_RESTORE_ROUNDF128 (FE_TONEAREST); 33 __float128 ret = x; 34 *eps = x_eps / x; 35 for (int i = 1; i < n; i++) 36 { 37 *eps += x_eps / (x + i); 38 __float128 lo; 39 mul_splitq (&ret, &lo, ret, x + i); 40 *eps += lo / ret; 41 } 42 return ret; 43 } 44