1 /* mpfr_gamma_inc -- incomplete gamma function
2
3 Copyright 2016-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 /* The incomplete gamma function is defined for x >= 0 and a not a negative
27 integer by:
28
29 gamma_inc(a,x) := Gamma(a,x) = int(t^(a-1) * exp(-t), t=x..infinity)
30
31 = Gamma(a) - gamma(a,x) with:
32
33 gamma(a,x) = int(t^(a-1) * exp(-t), t=0..x).
34
35 The function gamma(a,x) satisfies the Taylor expansions (we use the second
36 one in the code below):
37
38 gamma(a,x) = x^a * sum((-x)^k/k!/(a+k), k=0..infinity)
39
40 gamma(a,x) = x^a * exp(-x) * sum(x^k/(a*(a+1)*...*(a+k)), k=0..infinity)
41 */
42
43 static int
44 mpfr_gamma_inc_negint (mpfr_ptr y, mpfr_srcptr a, mpfr_srcptr x, mpfr_rnd_t r);
45
46 int
47 mpfr_gamma_inc (mpfr_ptr y, mpfr_srcptr a, mpfr_srcptr x, mpfr_rnd_t rnd)
48 {
49 mpfr_prec_t w;
50 mpfr_t s, t, u;
51 int inex;
52 unsigned long k;
53 mpfr_exp_t e0, e1, e2, err;
54 MPFR_GROUP_DECL(group);
55 MPFR_ZIV_DECL(loop);
56 MPFR_SAVE_EXPO_DECL (expo);
57
58 if (MPFR_ARE_SINGULAR (a, x))
59 {
60 /* if a or x is NaN, return NaN */
61 if (MPFR_IS_NAN (a) || MPFR_IS_NAN (x))
62 {
63 MPFR_SET_NAN (y);
64 MPFR_RET_NAN;
65 }
66
67 /* Note: for x < 0, gamma_inc (a, x) is a complex number */
68
69 if (MPFR_IS_INF (a) || MPFR_IS_INF (x))
70 {
71 if (MPFR_IS_INF (a) && MPFR_IS_INF (x))
72 {
73 if ((MPFR_IS_POS (a) && MPFR_IS_POS (x)) || MPFR_IS_NEG (x))
74 {
75 /* (a) gamma_inc(+Inf,+Inf) = NaN because
76 gamma_inc(x,x) tends to +Inf but
77 gamma_inc(x,x^2) tends to +0.
78 (b) gamma_inc(+/-Inf,-Inf) = NaN, for example
79 gamma_inc (a, -a) is a complex number
80 for a not an integer */
81 MPFR_SET_NAN (y);
82 MPFR_RET_NAN;
83 }
84 else
85 {
86 /* gamma_inc(-Inf,+Inf) = +0 */
87 MPFR_SET_ZERO (y);
88 MPFR_SET_POS (y);
89 MPFR_RET (0); /* exact */
90 }
91 }
92 else /* only one of a, x is infinite */
93 {
94 if (MPFR_IS_INF (a))
95 {
96 MPFR_ASSERTD (MPFR_IS_INF (a) && MPFR_IS_FP (x));
97 if (MPFR_IS_POS (a))
98 {
99 /* gamma_inc(+Inf, x) = +Inf */
100 MPFR_SET_INF (y);
101 MPFR_SET_POS (y);
102 MPFR_RET (0); /* exact */
103 }
104 else /* a = -Inf */
105 {
106 /* gamma_inc(-Inf, x) = NaN for x <= 0
107 +Inf for 0 < x < 1
108 +0 for 1 <= x */
109 if (mpfr_cmp_ui (x, 0) <= 0)
110 {
111 MPFR_SET_NAN (y);
112 MPFR_RET_NAN;
113 }
114 else if (mpfr_cmp_ui (x, 1) < 0)
115 {
116 MPFR_SET_INF (y);
117 MPFR_SET_POS (y);
118 MPFR_RET (0); /* exact */
119 }
120 else
121 {
122 MPFR_SET_ZERO (y);
123 MPFR_SET_POS (y);
124 MPFR_RET (0); /* exact */
125 }
126 }
127 }
128 else
129 {
130 MPFR_ASSERTD (MPFR_IS_FP (a) && MPFR_IS_INF (x));
131 if (MPFR_IS_POS (x))
132 {
133 /* x is +Inf: integral tends to zero */
134 MPFR_SET_ZERO (y);
135 MPFR_SET_POS (y);
136 MPFR_RET (0); /* exact */
137 }
138 else /* NaN for x < 0 */
139 {
140 MPFR_SET_NAN (y);
141 MPFR_RET_NAN;
142 }
143 }
144 }
145 }
146
147 if (MPFR_IS_ZERO (a) || MPFR_IS_ZERO (x))
148 {
149 if (MPFR_IS_ZERO (a))
150 {
151 if (mpfr_cmp_ui (x, 0) < 0)
152 {
153 /* gamma_inc(a,x) = NaN for x < 0 */
154 MPFR_SET_NAN (y);
155 MPFR_RET_NAN;
156 }
157 else if (MPFR_IS_ZERO (x))
158 /* gamma_inc(a,0) = gamma(a) */
159 return mpfr_gamma (y, a, rnd); /* a=+0->+Inf, a=-0->-Inf */
160 else
161 {
162 /* gamma_inc (0, x) = int (exp(-t), t=x..infinity) = E1(x) */
163 mpfr_t minus_x;
164 MPFR_TMP_INIT_NEG(minus_x, x);
165 /* mpfr_eint(x) for x < 0 returns -E1(-x) */
166 inex = mpfr_eint (y, minus_x, MPFR_INVERT_RND(rnd));
167 MPFR_CHANGE_SIGN(y);
168 return -inex;
169 }
170 }
171 else /* x = 0: gamma_inc(a,0) = gamma(a) */
172 return mpfr_gamma (y, a, rnd);
173 }
174 }
175
176 /* for x < 0 return NaN */
177 if (MPFR_SIGN(x) < 0)
178 {
179 MPFR_SET_NAN (y);
180 MPFR_RET_NAN;
181 }
182
183 if (mpfr_integer_p (a) && MPFR_SIGN(a) < 0)
184 return mpfr_gamma_inc_negint (y, a, x, rnd);
185
186 MPFR_SAVE_EXPO_MARK (expo);
187
188 w = MPFR_PREC(y) + 13; /* working precision */
189
190 MPFR_GROUP_INIT_2(group, w, s, t);
191 mpfr_init2 (u, 2); /* u is special (see below) */
192 MPFR_ZIV_INIT (loop, w);
193 for (;;)
194 {
195 mpfr_exp_t expu, precu, exps;
196 mpfr_t s_abs;
197 mpfr_exp_t decay = 0;
198 MPFR_BLOCK_DECL (flags);
199
200 /* Note: in the error analysis below, theta represents any value of
201 absolute value less than 2^(-w) where w is the working precision (two
202 instances of theta may represent different values), cf Higham's book.
203 */
204
205 /* to ensure that u = a + k is exact, we have three cases:
206 (1) EXP(a) <= 0, then we need PREC(u) >= 1 - EXP(a) + PREC(a)
207 (2) EXP(a) - PREC(a) <= 0 < E(a), then PREC(u) >= PREC(a)
208 (3) 0 < EXP(a) - PREC(a), then PREC(u) >= EXP(a) */
209 precu = MPFR_GET_EXP(a) <= 0 ?
210 MPFR_ADD_PREC (MPFR_PREC(a), 1 - MPFR_EXP(a))
211 : (MPFR_EXP(a) <= MPFR_PREC(a)) ? MPFR_PREC(a) : MPFR_EXP(a);
212 MPFR_ASSERTD (precu + 1 <= MPFR_PREC_MAX);
213 mpfr_set_prec (u, precu + 1);
214 expu = (MPFR_EXP(a) > 0) ? MPFR_EXP(a) : 1;
215
216 /* estimate Taylor series */
217 mpfr_ui_div (t, 1, a, MPFR_RNDA); /* t = 1/a * (1 + theta) */
218 mpfr_set (s, t, MPFR_RNDA); /* s = 1/a * (1 + theta) */
219 if (MPFR_IS_NEG(a))
220 {
221 mpfr_init2 (s_abs, 32);
222 mpfr_abs (s_abs, s, MPFR_RNDU);
223 }
224 for (k = 1;; k++)
225 {
226 mpfr_mul (t, t, x, MPFR_RNDU); /* t = x^k/(a * ... * (a+k-1))
227 * (1 + theta)^(2k) */
228 inex = mpfr_add_ui (u, a, k, MPFR_RNDZ); /* u = a+k exactly */
229 MPFR_ASSERTD(inex == 0);
230 mpfr_div (t, t, u, MPFR_RNDA); /* t = x^k/(a * ... * (a+k))
231 * (1 + theta)^(2k+1) */
232 mpfr_add (s, s, t, MPFR_RNDZ);
233 /* when s is zero, we consider ulp(s) = ulp(t) */
234 exps = (MPFR_IS_ZERO(s)) ? MPFR_GET_EXP(t) : MPFR_GET_EXP(s);
235 if (MPFR_IS_NEG(a))
236 {
237 if (MPFR_IS_POS(t))
238 mpfr_add (s_abs, s_abs, t, MPFR_RNDU);
239 else
240 mpfr_sub (s_abs, s_abs, t, MPFR_RNDU);
241 }
242 /* we stop when |t| < ulp(s), u > 0 and |x/u| < 1/2, which ensures
243 that the tail is at most 2*ulp(s) */
244 MPFR_ASSERTD (MPFR_NOTZERO(t));
245 if (MPFR_GET_EXP(t) + w <= exps && MPFR_IS_POS(u) &&
246 MPFR_GET_EXP(x) + 1 < MPFR_GET_EXP(u))
247 break;
248
249 /* if there was an exponent shift in u, increase the precision of
250 u so that mpfr_add_ui (u, a, k) remains exact */
251 if (MPFR_EXP(u) > expu) /* exponent shift in u */
252 {
253 MPFR_ASSERTD(MPFR_EXP(u) == expu + 1);
254 expu = MPFR_EXP(u);
255 mpfr_set_prec (u, mpfr_get_prec (u) + 1);
256 }
257 }
258 if (MPFR_IS_NEG(a))
259 {
260 decay = MPFR_GET_EXP(s_abs) - MPFR_GET_EXP(s);
261 mpfr_clear (s_abs);
262 }
263 /* For a > 0, since all terms are positive, we have
264 s = S * (1 + theta)^(2k+3) with S being the infinite Taylor series.
265 For a < 0, the error is bounded by that on the sum s_abs of absolute
266 values of the terms, i.e., S_abs * [(1 + theta)^(2k+3) - 1]. Thus we
267 can simply use the same error analysis as for a > 0, adding an error
268 corresponding to the decay of exponent between s_abs and s. */
269
270 /* multiply by exp(-x) */
271 mpfr_exp (t, x, MPFR_RNDZ); /* t = exp(x) * (1+theta) */
272 mpfr_div (s, s, t, MPFR_RNDZ); /* s = <exact value> * (1+theta)^(2k+5) */
273
274 /* multiply by x^a */
275 mpfr_pow (t, x, a, MPFR_RNDZ); /* t = x^a * (1+theta) */
276 mpfr_mul (s, s, t, MPFR_RNDZ); /* s = Gamma(a,x) * (1+theta)^(2k+7) */
277
278 /* Since |theta| < 2^(-w) using the Taylor expansion of log(1+x)
279 we have log(1+theta) = theta1 with |theta1| < 1.16*2^(-w) for w >= 2,
280 thus (1+theta)^(2k+7) = exp((2k+7)*theta1).
281 Assuming 2k+7 = t*2^w for |t| < 0.5, we have
282 |(2k+7)*theta1| = |t*2^w*theta1| < 0.58.
283 For |u| < 0.58 we have |exp(u)-1| < 1.36*|u|
284 thus |(1+theta)^(2k+7) - 1| < 1.36*0.58*(2k+7)/2^w < 0.79*(2k+7)/2^w.
285 Since one ulp is at worst a relative error of 2^(1-w),
286 the error on s is at most 2^(decay+1)*(2k+7) ulps. */
287
288 /* subtract from gamma(a) */
289 MPFR_BLOCK (flags, mpfr_gamma (t, a, MPFR_RNDZ));
290 MPFR_ASSERTN (!MPFR_OVERFLOW (flags)); /* FIXME: support overflow */
291 /* t = gamma(a) * (1+theta) */
292 e0 = MPFR_GET_EXP (t);
293 e1 = (MPFR_IS_ZERO(s)) ? __gmpfr_emin : MPFR_GET_EXP (s);
294 mpfr_sub (s, t, s, MPFR_RNDZ);
295 /* if s is zero, we can assume ulp(s) = ulp(t), but anyway we won't
296 be able to round */
297 e2 = (MPFR_IS_ZERO(s)) ? e0 : MPFR_GET_EXP (s);
298 /* the final error is at most 1 ulp (for the final subtraction)
299 + 2^(e0-e2) ulps # for the error in t
300 + 2^(decay+1)*(2k+7) ulps * 2^(e1-e2) # for the error in gamma(a,x) */
301
302 e1 += decay + 1 + MPFR_INT_CEIL_LOG2 (2*k+7);
303 /* Now the error is <= 1 + 2^(e0-e2) + 2^(e1-e2).
304 Since the formula is symmetric in e0 and e1, we can assume without
305 loss of generality e0 >= e1, then:
306 if e0 = e1: err <= 1 + 2*2^(e0-e2) <= 2^(e0-e2+2)
307 if e0 > e1: err <= 1 + 1.5*2^(e0-e2)
308 <= 2^(e0-e2+1) if e0 > e2
309 <= 2^2 otherwise */
310 if (e0 == e1)
311 {
312 /* Check that e0 - e2 + 2 <= MPFR_EXP_MAX */
313 MPFR_ASSERTD (e2 >= 2 || e0 <= (MPFR_EXP_MAX - 2) + e2);
314 /* Check that e0 - e2 + 2 >= MPFR_EXP_MIN */
315 MPFR_ASSERTD (e2 <= 2 || e0 >= MPFR_EXP_MIN + (e2 - 2));
316 err = e0 - e2 + 2;
317 }
318 else
319 {
320 e0 = (e0 > e1) ? e0 : e1; /* max(e0,e1) */
321 MPFR_ASSERTD (e0 <= e2 || e2 >= 1 || e0 <= (MPFR_EXP_MAX - 1) + e2);
322 err = (e0 > e2) ? e0 - e2 + 1 : 2;
323 }
324
325 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, w - err, MPFR_PREC(y), rnd)))
326 break;
327
328 MPFR_ZIV_NEXT (loop, w);
329 MPFR_GROUP_REPREC_2(group, w, s, t);
330 }
331 MPFR_ZIV_FREE (loop);
332 mpfr_clear (u);
333
334 inex = mpfr_set (y, s, rnd);
335 MPFR_GROUP_CLEAR(group);
336
337 MPFR_SAVE_EXPO_FREE (expo);
338 return mpfr_check_range (y, inex, rnd);
339 }
340
341 /* For a negative integer, we have (formula 6.5.19):
342
343 gamma(-n,x) = (-1)^n/n! [E_1(x) - exp(-x) sum((-1)^j*j!/x^(j+1), j=0..n-1)]
344
345 See also https://arxiv.org/pdf/1407.0349v1.pdf.
346
347 Assumes 'a' is a negative integer.
348 */
349 static int
350 mpfr_gamma_inc_negint (mpfr_ptr y, mpfr_srcptr a, mpfr_srcptr x,
351 mpfr_rnd_t rnd)
352 {
353 mpfr_t s, t, abs_a, neg_x;
354 unsigned long j;
355 mpfr_prec_t w;
356 int inex;
357 mpfr_exp_t exp_s, new_exp_s, exp_t, err_s, logj;
358 MPFR_GROUP_DECL(group);
359 MPFR_ZIV_DECL(loop);
360 MPFR_SAVE_EXPO_DECL (expo);
361
362 MPFR_ASSERTD(mpfr_integer_p (a));
363 MPFR_ASSERTD(mpfr_cmp_ui (a, 0) < 0);
364
365 MPFR_TMP_INIT_ABS(abs_a, a);
366
367 /* below, theta represents any value such that |theta| <= 2^(-w) */
368
369 w = MPFR_PREC(y) + 10; /* initial working precision */
370
371 MPFR_SAVE_EXPO_MARK (expo);
372 MPFR_GROUP_INIT_2(group, w, s, t);
373 MPFR_ZIV_INIT (loop, w);
374 for (;;)
375 {
376 /* we require |a| <= 2^(w-3) for the error analysis below */
377 if (MPFR_GET_EXP(a) + 3 > w)
378 w = MPFR_GET_EXP(a) + 3;
379
380 mpfr_ui_div (t, 1, x, MPFR_RNDN); /* t = 1/x * (1 + theta) */
381 mpfr_set (s, t, MPFR_RNDN);
382 MPFR_ASSERTD (MPFR_NOTZERO(s));
383 exp_t = exp_s = MPFR_GET_EXP(s); /* max. exponent of s/t during loop */
384 new_exp_s = exp_s;
385
386 for (j = 1; mpfr_cmp_ui (abs_a, j) > 0; j++)
387 {
388 /* invariant: t = (-1)^(j-1)*(j-1)!/x^j * (1 + theta)^(2j-1) */
389 mpfr_mul_ui (t, t, j, MPFR_RNDN);
390 mpfr_neg (t, t, MPFR_RNDN); /* exact */
391 mpfr_div (t, t, x, MPFR_RNDN);
392 /* now t = (-1)^j*j!/x^(j+1) * (1 + theta)^(2j+1).
393 We have (1 + theta)^(2j+1) = exp((2j+1)*log(1+theta)).
394 For |u| <= 1/2, we have |log(1+u)| <= 1.4 |u| thus:
395 |(1+theta)^(2j+1)-1| <= max |exp(1.4*(2j+1)*u)-1| for |u|<=2^(-w).
396 Now for |v| <= 1/2 we have |exp(v)-1| <= 0.7*|v| thus:
397 |(1+theta)^(2j+1) - 1| <= 2*(2j+1)*2^(-w)
398 as long as 1.4*(2j+1)*2^(-w) <= 1/2, which is true when j<2^(w-3).
399 Since j < |a| it suffices that |a| <= 2^(w-3).
400 In that case the rel. error on t is bounded by 2*(2j+1)*2^(-w),
401 thus the error in ulps is bounded by 2*(2j+1) ulp(t). */
402 if (MPFR_IS_ZERO(t)) /* underflow on t */
403 break;
404 if (MPFR_GET_EXP(t) > exp_t)
405 exp_t = MPFR_GET_EXP(t);
406 mpfr_add (s, s, t, MPFR_RNDN);
407 /* if s is zero, we can assume its ulp is that of t */
408 new_exp_s = (MPFR_IS_ZERO(s)) ? MPFR_GET_EXP(t) : MPFR_GET_EXP(s);
409 if (new_exp_s > exp_s)
410 exp_s = new_exp_s;
411 }
412
413 /* the error on s is bounded by (j-1) * 2^(exp_s - EXP(s)) * 1/2
414 for the mpfr_add roundings, plus
415 sum(2*(2i+1), i=1..j-1) * 2^(exp_t - EXP(s)) for the error on t.
416 The latter sum is (2*j^2-2) * 2^(exp_t - EXP(s)). */
417
418 logj = MPFR_INT_CEIL_LOG2(j);
419 exp_s += logj - 1;
420 exp_t += 1 + 2 * logj;
421
422 /* now the error on s is bounded by 2^(exp_s-EXP(s))+2^(exp_t-EXP(s)) */
423
424 exp_s = (exp_s >= exp_t) ? exp_s + 1 : exp_t + 1;
425 err_s = exp_s - new_exp_s;
426
427 /* now the error on the sum S := sum((-1)^j*j!/x^(j+1), j=0..n-1)
428 is bounded by 2^err_s ulp(s) */
429
430 MPFR_TMP_INIT_NEG(neg_x, x);
431
432 mpfr_exp (t, neg_x, MPFR_RNDN); /* t = exp(-x) * (1 + theta) */
433 mpfr_mul (s, s, t, MPFR_RNDN);
434 if (MPFR_IS_ZERO(s))
435 {
436 MPFR_ASSERTD (MPFR_NOTZERO(t));
437 new_exp_s += MPFR_GET_EXP(t);
438 }
439 /* s = exp(-x) * (S +/- 2^err_s ulp(S)) * (1 + theta)^2.
440 = exp(-x) * (S +/- 2^err_s ulp(S)) * (1 +/- 3 ulp(1))
441 The error on s is bounded by:
442 exp(-x) * [2^err_s*ulp(S) + S*3*ulp(1) + 2^err_s*ulp(S)*3*ulp(1)]
443 <= ulp(s) * [2^(err_s+1) + 6 + 1]
444 <= ulp(s) * 2^(err_s+2) as long as err_s >= 2. */
445
446 err_s = (err_s >= 2) ? err_s + 2 : 4;
447 /* now the error on s is bounded by 2^err_s ulp(s) */
448
449 mpfr_eint (t, neg_x, MPFR_RNDN); /* t = -E1(-x) * (1 + theta) */
450 mpfr_neg (t, t, MPFR_RNDN); /* exact */
451
452 exp_s = (MPFR_IS_ZERO(s)) ? new_exp_s : MPFR_GET_EXP(s);
453 MPFR_ASSERTD (MPFR_NOTZERO(t));
454 exp_t = MPFR_GET_EXP(t);
455 mpfr_sub (s, t, s, MPFR_RNDN); /* E_1(x) - exp(-x) * S */
456 if (MPFR_IS_ZERO(s)) /* cancellation: increase working precision */
457 goto next_w;
458
459 /* err(s) <= 1/2 * ulp(s) [mpfr_sub]
460 + 2^err_s * 2^(exp_s-EXP(s)) * ulp(s) [previous error on s]
461 + 1/2 * 2^(exp_t-EXP(s)) * ulp(s) [error on t] */
462
463 exp_s += err_s;
464 exp_t -= 1;
465 exp_s = (exp_s >= exp_t) ? exp_s + 1 : exp_t + 1;
466 MPFR_ASSERTD (MPFR_NOTZERO(s));
467 err_s = exp_s - MPFR_GET_EXP(s);
468 /* err(s) <= 1/2 * ulp(s) + 2^err_s * ulp(s) */
469
470 /* divide by n! */
471 mpfr_gamma (t, abs_a, MPFR_RNDN); /* t = (n-1)! * (1 + theta) */
472 mpfr_mul (t, t, abs_a, MPFR_RNDN); /* t = n! * (1 + theta)^2 */
473 mpfr_div (s, s, t, MPFR_RNDN);
474 /* since (1 + theta)^2 converts to an error of at most 3 ulps
475 for w >= 2, the final error is at most:
476 2 * (1/2 + 2^err_s) * ulp(s) [error on previous s]
477 + 2 * 3 * ulp(s) [error on t]
478 + 1 * ulp(s) [product of errors]
479 = ulp(s) * (2^(err_s+1) + 8) */
480 err_s = (err_s >= 2) ? err_s + 1 : 4;
481
482 /* the final error is bounded by 2^err_s * ulp(s) */
483
484 /* Is there a better way to compute (-1)^n? */
485 mpfr_set_si (t, -1, MPFR_RNDN);
486 mpfr_pow (t, t, abs_a, MPFR_RNDN);
487 if (MPFR_IS_NEG(t))
488 mpfr_neg (s, s, MPFR_RNDN);
489
490 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, w - err_s, MPFR_PREC(y), rnd)))
491 break;
492
493 next_w:
494 MPFR_ZIV_NEXT (loop, w);
495 MPFR_GROUP_REPREC_2(group, w, s, t);
496 }
497 MPFR_ZIV_FREE (loop);
498
499 inex = mpfr_set (y, s, rnd);
500 MPFR_GROUP_CLEAR(group);
501
502 MPFR_SAVE_EXPO_FREE (expo);
503 return mpfr_check_range (y, inex, rnd);
504 }