]>
Commit | Line | Data |
---|---|---|
050f29c1 | 1 | /* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... |
04277e02 | 2 | Copyright (C) 2015-2019 Free Software Foundation, Inc. |
050f29c1 JM |
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 | |
5a82c748 | 17 | <https://www.gnu.org/licenses/>. */ |
050f29c1 JM |
18 | |
19 | #include <math.h> | |
20 | #include <math_private.h> | |
4482ff22 | 21 | #include <mul_splitl.h> |
050f29c1 JM |
22 | |
23 | /* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + | |
24 | 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that | |
25 | all the values X + 1, ..., X + N - 1 are exactly representable, and | |
26 | X_EPS / X is small enough that factors quadratic in it can be | |
27 | neglected. */ | |
28 | ||
29 | long double | |
30 | __lgamma_productl (long double t, long double x, long double x_eps, int n) | |
31 | { | |
32 | long double ret = 0, ret_eps = 0; | |
33 | for (int i = 0; i < n; i++) | |
34 | { | |
35 | long double xi = x + i; | |
36 | long double quot = t / xi; | |
37 | long double mhi, mlo; | |
4482ff22 | 38 | mul_splitl (&mhi, &mlo, quot, xi); |
050f29c1 JM |
39 | long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); |
40 | /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ | |
41 | long double rhi, rlo; | |
4482ff22 | 42 | mul_splitl (&rhi, &rlo, ret, quot); |
050f29c1 JM |
43 | long double rpq = ret + quot; |
44 | long double rpq_eps = (ret - rpq) + quot; | |
45 | long double nret = rpq + rhi; | |
46 | long double nret_eps = (rpq - nret) + rhi; | |
47 | ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot | |
48 | + quot_lo + quot_lo * (ret + ret_eps)); | |
49 | ret = nret; | |
50 | } | |
51 | return ret + ret_eps; | |
52 | } |