1 /* mpfr_y0, mpfr_y1, mpfr_yn -- Bessel functions of 2nd kind, integer order.
2 http://www.opengroup.org/onlinepubs/009695399/functions/y0.html
4 Copyright 2007-2015 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
27 static int mpfr_yn_asympt (mpfr_ptr
, long, mpfr_srcptr
, mpfr_rnd_t
);
30 mpfr_y0 (mpfr_ptr res
, mpfr_srcptr z
, mpfr_rnd_t r
)
32 return mpfr_yn (res
, 0, z
, r
);
36 mpfr_y1 (mpfr_ptr res
, mpfr_srcptr z
, mpfr_rnd_t r
)
38 return mpfr_yn (res
, 1, z
, r
);
41 /* compute in s an approximation of S1 = sum((n-k)!/k!*y^k,k=0..n)
42 return e >= 0 the exponent difference between the maximal value of |s|
43 during the for loop and the final value of |s|.
46 mpfr_yn_s1 (mpfr_ptr s
, mpfr_srcptr y
, unsigned long n
)
52 mpz_init_set_ui (f
, 1);
53 /* we compute n!*S1 = sum(a[k]*y^k,k=0..n) where a[k] = n!*(n-k)!/k!,
54 a[0] = (n!)^2, a[1] = n!*(n-1)!, ..., a[n-1] = n, a[n] = 1 */
55 mpfr_set_ui (s
, 1, MPFR_RNDN
); /* a[n] */
59 /* a[k]/a[k+1] = (n-k)!/k!/(n-(k+1))!*(k+1)! = (k+1)*(n-k) */
60 mpfr_mul (s
, s
, y
, MPFR_RNDN
);
61 mpz_mul_ui (f
, f
, n
- k
);
62 mpz_mul_ui (f
, f
, k
+ 1);
63 /* invariant: f = a[k] */
64 mpfr_add_z (s
, s
, f
, MPFR_RNDN
);
69 /* now we have f = (n!)^2 */
71 mpfr_div_z (s
, s
, f
, MPFR_RNDN
);
73 return emax
- MPFR_EXP(s
);
76 /* compute in s an approximation of
77 S3 = c*sum((h(k)+h(n+k))*y^k/k!/(n+k)!,k=0..infinity)
78 where h(k) = 1 + 1/2 + ... + 1/k
82 Returns e such that the error is bounded by 2^e ulp(s).
85 mpfr_yn_s3 (mpfr_ptr s
, mpfr_srcptr y
, mpfr_srcptr c
, unsigned long n
)
89 mpz_t p
, q
; /* p/q will store h(k)+h(n+k) */
90 mpfr_exp_t exps
, expU
;
92 zz
= mpfr_get_ui (y
, MPFR_RNDU
); /* y = z^2/4 */
93 MPFR_ASSERTN (zz
< ULONG_MAX
- 2);
94 zz
+= 2; /* z^2 <= 2^zz */
95 mpz_init_set_ui (p
, 0);
96 mpz_init_set_ui (q
, 1);
97 /* initialize p/q to h(n) */
98 for (k
= 1; k
<= n
; k
++)
100 /* p/q + 1/k = (k*p+q)/(q*k) */
101 mpz_mul_ui (p
, p
, k
);
103 mpz_mul_ui (q
, q
, k
);
105 mpfr_init2 (t
, MPFR_PREC(s
));
106 mpfr_init2 (u
, MPFR_PREC(s
));
107 mpfr_fac_ui (t
, n
, MPFR_RNDN
);
108 mpfr_div (t
, c
, t
, MPFR_RNDN
); /* c/n! */
109 mpfr_mul_z (u
, t
, p
, MPFR_RNDN
);
110 mpfr_div_z (s
, u
, q
, MPFR_RNDN
);
116 mpfr_mul (t
, t
, y
, MPFR_RNDN
);
117 mpfr_div_ui (t
, t
, k
, MPFR_RNDN
);
118 mpfr_div_ui (t
, t
, n
+ k
, MPFR_RNDN
);
120 p/q + 1/k + 1/(n+k) = [p*k*(n+k) + q*(n+k) + q*k]/(q*k*(n+k)) */
121 mpz_mul_ui (p
, p
, k
);
122 mpz_mul_ui (p
, p
, n
+ k
);
123 mpz_addmul_ui (p
, q
, n
+ 2 * k
);
124 mpz_mul_ui (q
, q
, k
);
125 mpz_mul_ui (q
, q
, n
+ k
);
126 mpfr_mul_z (u
, t
, p
, MPFR_RNDN
);
127 mpfr_div_z (u
, u
, q
, MPFR_RNDN
);
131 mpfr_add (s
, s
, u
, MPFR_RNDN
);
135 if (MPFR_EXP (u
) + (mpfr_exp_t
) MPFR_PREC (u
) < MPFR_EXP (s
) &&
136 zz
/ (2 * k
) < k
+ n
)
143 exps
= expU
- MPFR_EXP (s
);
144 /* the error is bounded by (6k^2+33/2k+11) 2^exps ulps
145 <= 8*(k+2)^2 2^exps ulps */
146 return 3 + 2 * MPFR_INT_CEIL_LOG2(k
+ 2) + exps
;
150 mpfr_yn (mpfr_ptr res
, long n
, mpfr_srcptr z
, mpfr_rnd_t r
)
154 MPFR_SAVE_EXPO_DECL (expo
);
157 (("n=%ld x[%Pu]=%.*Rg rnd=%d", n
, mpfr_get_prec (z
), mpfr_log_prec
, z
, r
),
158 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (res
), mpfr_log_prec
, res
, inex
));
160 absn
= SAFE_ABS (unsigned long, n
);
162 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z
)))
166 MPFR_SET_NAN (res
); /* y(n,NaN) = NaN */
169 /* y(n,z) tends to zero when z goes to +Inf, oscillating around
170 0. We choose to return +0 in that case. */
171 else if (MPFR_IS_INF (z
))
173 if (MPFR_SIGN(z
) > 0)
174 return mpfr_set_ui (res
, 0, r
);
175 else /* y(n,-Inf) = NaN */
181 else /* y(n,z) tends to -Inf for n >= 0 or n even, to +Inf otherwise,
182 when z goes to zero */
185 if (n
>= 0 || ((unsigned long) n
& 1) == 0)
194 /* for z < 0, y(n,z) is imaginary except when j(n,|z|) = 0, which we
195 assume does not happen for a rational z. */
196 if (MPFR_SIGN(z
) < 0)
202 /* now z is not singular, and z > 0 */
204 MPFR_SAVE_EXPO_MARK (expo
);
206 /* Deal with tiny arguments. We have:
207 y0(z) = 2 log(z)/Pi + 2 (euler - log(2))/Pi + O(log(z)*z^2), more
208 precisely for 0 <= z <= 1/2, with g(z) = 2/Pi + 2(euler-log(2))/Pi/log(z),
209 g(z) - 0.41*z^2 < y0(z)/log(z) < g(z)
210 thus since log(z) is negative:
211 g(z)*log(z) < y0(z) < (g(z) - z^2/2)*log(z)
212 and since |g(z)| >= 0.63 for 0 <= z <= 1/2, the relative error on
213 y0(z)/log(z) is bounded by 0.41*z^2/0.63 <= 0.66*z^2.
214 Note: we use both the main term in log(z) and the constant term, because
215 otherwise the relative error would be only in 1/log(|log(z)|).
217 if (n
== 0 && MPFR_EXP(z
) < - (mpfr_exp_t
) (MPFR_PREC(res
) / 2))
219 mpfr_t l
, h
, t
, logz
;
223 prec
= MPFR_PREC(res
) + 10;
224 mpfr_init2 (l
, prec
);
225 mpfr_init2 (h
, prec
);
226 mpfr_init2 (t
, prec
);
227 mpfr_init2 (logz
, prec
);
228 /* first enclose log(z) + euler - log(2) = log(z/2) + euler */
229 mpfr_log (logz
, z
, MPFR_RNDD
); /* lower bound of log(z) */
230 mpfr_set (h
, logz
, MPFR_RNDU
); /* exact */
231 mpfr_nextabove (h
); /* upper bound of log(z) */
232 mpfr_const_euler (t
, MPFR_RNDD
); /* lower bound of euler */
233 mpfr_add (l
, logz
, t
, MPFR_RNDD
); /* lower bound of log(z) + euler */
234 mpfr_nextabove (t
); /* upper bound of euler */
235 mpfr_add (h
, h
, t
, MPFR_RNDU
); /* upper bound of log(z) + euler */
236 mpfr_const_log2 (t
, MPFR_RNDU
); /* upper bound of log(2) */
237 mpfr_sub (l
, l
, t
, MPFR_RNDD
); /* lower bound of log(z/2) + euler */
238 mpfr_nextbelow (t
); /* lower bound of log(2) */
239 mpfr_sub (h
, h
, t
, MPFR_RNDU
); /* upper bound of log(z/2) + euler */
240 mpfr_const_pi (t
, MPFR_RNDU
); /* upper bound of Pi */
241 mpfr_div (l
, l
, t
, MPFR_RNDD
); /* lower bound of (log(z/2)+euler)/Pi */
242 mpfr_nextbelow (t
); /* lower bound of Pi */
243 mpfr_div (h
, h
, t
, MPFR_RNDD
); /* upper bound of (log(z/2)+euler)/Pi */
244 mpfr_mul_2ui (l
, l
, 1, MPFR_RNDD
); /* lower bound on g(z)*log(z) */
245 mpfr_mul_2ui (h
, h
, 1, MPFR_RNDU
); /* upper bound on g(z)*log(z) */
246 /* we now have l <= g(z)*log(z) <= h, and we need to add -z^2/2*log(z)
248 mpfr_mul (t
, z
, z
, MPFR_RNDU
); /* upper bound on z^2 */
249 /* since logz is negative, a lower bound corresponds to an upper bound
250 for its absolute value */
251 mpfr_neg (t
, t
, MPFR_RNDD
);
252 mpfr_div_2ui (t
, t
, 1, MPFR_RNDD
);
253 mpfr_mul (t
, t
, logz
, MPFR_RNDU
); /* upper bound on z^2/2*log(z) */
254 mpfr_add (h
, h
, t
, MPFR_RNDU
);
255 inex
= mpfr_prec_round (l
, MPFR_PREC(res
), r
);
256 inex2
= mpfr_prec_round (h
, MPFR_PREC(res
), r
);
257 /* we need h=l and inex=inex2 */
258 ok
= (inex
== inex2
) && mpfr_equal_p (l
, h
);
260 mpfr_set (res
, h
, r
); /* exact */
269 /* small argument check for y1(z) = -2/Pi/z + O(log(z)):
270 for 0 <= z <= 1, |y1(z) + 2/Pi/z| <= 0.25 */
271 if (n
== 1 && MPFR_EXP(z
) + 1 < - (mpfr_exp_t
) MPFR_PREC(res
))
277 MPFR_BLOCK_DECL (flags
);
279 /* since 2/Pi > 0.5, and |y1(z)| >= |2/Pi/z|, if z <= 2^(-emax-1),
280 then |y1(z)| > 2^emax */
281 prec
= MPFR_PREC(res
) + 10;
282 mpfr_init2 (y
, prec
);
283 mpfr_const_pi (y
, MPFR_RNDU
); /* Pi*(1+u)^2, where here and below u
284 represents a quantity <= 1/2^prec */
285 mpfr_mul (y
, y
, z
, MPFR_RNDU
); /* Pi*z * (1+u)^4, upper bound */
286 MPFR_BLOCK (flags
, mpfr_ui_div (y
, 2, y
, MPFR_RNDZ
));
287 /* 2/Pi/z * (1+u)^6, lower bound, with possible overflow */
288 if (MPFR_OVERFLOW (flags
))
291 MPFR_SAVE_EXPO_FREE (expo
);
292 return mpfr_overflow (res
, r
, -1);
294 mpfr_neg (y
, y
, MPFR_RNDN
);
295 /* (1+u)^6 can be written 1+7u [for another value of u], thus the
296 error on 2/Pi/z is less than 7ulp(y). The truncation error is less
297 than 1/4, thus if ulp(y)>=1/4, the total error is less than 8ulp(y),
298 otherwise it is less than 1/4+7/8 <= 2. */
299 if (MPFR_EXP(y
) + 2 >= MPFR_PREC(y
)) /* ulp(y) >= 1/4 */
301 else /* ulp(y) <= 1/8 */
302 err1
= (mpfr_exp_t
) MPFR_PREC(y
) - MPFR_EXP(y
) + 1;
303 ok
= MPFR_CAN_ROUND (y
, prec
- err1
, MPFR_PREC(res
), r
);
305 inex
= mpfr_set (res
, y
, r
);
311 /* we can use the asymptotic expansion as soon as z > p log(2)/2,
312 but to get some margin we use it for z > p/2 */
313 if (mpfr_cmp_ui (z
, MPFR_PREC(res
) / 2 + 3) > 0)
315 inex
= mpfr_yn_asympt (res
, n
, z
, r
);
323 mpfr_exp_t err1
, err2
, err3
;
324 mpfr_t y
, s1
, s2
, s3
;
325 MPFR_ZIV_DECL (loop
);
332 prec
= MPFR_PREC(res
) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res
)) + 13;
333 MPFR_ZIV_INIT (loop
, prec
);
336 mpfr_set_prec (y
, prec
);
337 mpfr_set_prec (s1
, prec
);
338 mpfr_set_prec (s2
, prec
);
339 mpfr_set_prec (s3
, prec
);
341 mpfr_mul (y
, z
, z
, MPFR_RNDN
);
342 mpfr_div_2ui (y
, y
, 2, MPFR_RNDN
); /* z^2/4 */
344 /* store (z/2)^n temporarily in s2 */
345 mpfr_pow_ui (s2
, z
, absn
, MPFR_RNDN
);
346 mpfr_div_2si (s2
, s2
, absn
, MPFR_RNDN
);
348 /* compute S1 * (z/2)^(-n) */
351 mpfr_set_ui (s1
, 0, MPFR_RNDN
);
355 err1
= mpfr_yn_s1 (s1
, y
, absn
- 1);
356 mpfr_div (s1
, s1
, s2
, MPFR_RNDN
); /* (z/2)^(-n) * S1 */
357 /* See algorithms.tex: the relative error on s1 is bounded by
358 (3n+3)*2^(e+1-prec). */
359 err1
= MPFR_INT_CEIL_LOG2 (3 * absn
+ 3) + err1
+ 1;
360 /* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */
362 /* compute (z/2)^n * S3 */
363 mpfr_neg (y
, y
, MPFR_RNDN
); /* -z^2/4 */
364 err3
= mpfr_yn_s3 (s3
, y
, s2
, absn
); /* (z/2)^n * S3 */
365 /* the error on s3 is bounded by 2^err3 ulps */
368 err1
+= MPFR_EXP(s1
);
369 mpfr_add (s1
, s1
, s3
, MPFR_RNDN
);
370 /* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1))
371 + 2^err3*2^(EXP(s3) - EXP(s1)) */
372 err3
+= MPFR_EXP(s3
);
373 err1
= (err3
> err1
) ? err3
+ 1 : err1
+ 1;
374 err1
-= MPFR_EXP(s1
);
375 err1
= (err1
>= 0) ? err1
+ 1 : 1;
376 /* now the error on s1 is bounded by 2^err1*ulp(s1) */
379 mpfr_div_2ui (s2
, z
, 1, MPFR_RNDN
); /* z/2 */
380 mpfr_log (s2
, s2
, MPFR_RNDN
); /* log(z/2) */
381 mpfr_const_euler (s3
, MPFR_RNDN
);
382 err2
= MPFR_EXP(s2
) > MPFR_EXP(s3
) ? MPFR_EXP(s2
) : MPFR_EXP(s3
);
383 mpfr_add (s2
, s2
, s3
, MPFR_RNDN
); /* log(z/2) + gamma */
384 err2
-= MPFR_EXP(s2
);
385 mpfr_mul_2ui (s2
, s2
, 1, MPFR_RNDN
); /* 2*(log(z/2) + gamma) */
386 mpfr_jn (s3
, absn
, z
, MPFR_RNDN
); /* Jn(z) */
387 mpfr_mul (s2
, s2
, s3
, MPFR_RNDN
); /* 2*(log(z/2) + gamma)*Jn(z) */
388 err2
+= 4; /* the error on s2 is bounded by 2^err2 ulps, see
391 /* add all three sums */
392 err1
+= MPFR_EXP(s1
); /* the error on s1 is bounded by 2^err1 */
393 err2
+= MPFR_EXP(s2
); /* the error on s2 is bounded by 2^err2 */
394 mpfr_sub (s2
, s2
, s1
, MPFR_RNDN
); /* s2 - (s1+s3) */
395 err2
= (err1
> err2
) ? err1
+ 1 : err2
+ 1;
396 err2
-= MPFR_EXP(s2
);
397 err2
= (err2
>= 0) ? err2
+ 1 : 1;
398 /* now the error on s2 is bounded by 2^err2*ulp(s2) */
399 mpfr_const_pi (y
, MPFR_RNDN
); /* error bounded by 1 ulp */
400 mpfr_div (s2
, s2
, y
, MPFR_RNDN
); /* error bounded by
401 2^(err2+1)*ulp(s2) */
404 if (MPFR_LIKELY (MPFR_CAN_ROUND (s2
, prec
- err2
, MPFR_PREC(res
), r
)))
406 MPFR_ZIV_NEXT (loop
, prec
);
408 MPFR_ZIV_FREE (loop
);
410 /* Assume two's complement for the test n & 1 */
411 inex
= mpfr_set4 (res
, s2
, r
, n
>= 0 || (n
& 1) == 0 ?
412 MPFR_SIGN (s2
) : - MPFR_SIGN (s2
));
421 MPFR_SAVE_EXPO_FREE (expo
);
422 return mpfr_check_range (res
, inex
, r
);
426 #include "jyn_asympt.c"