mirror of
https://review.haiku-os.org/buildtools
synced 2025-01-19 12:51:22 +01:00
b58ddff026
* merged mpfr 3.0.0 and gmp 5.0.1 in buildtools trunk git-svn-id: file:///srv/svn/repos/haiku/buildtools/trunk@37378 a95241bf-73f2-0310-859d-f6bbb57e9c96
421 lines
14 KiB
C
421 lines
14 KiB
C
/* mpfr_y0, mpfr_y1, mpfr_yn -- Bessel functions of 2nd kind, integer order.
|
|
http://www.opengroup.org/onlinepubs/009695399/functions/y0.html
|
|
|
|
Copyright 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
|
|
Contributed by the Arenaire and Cacao projects, INRIA.
|
|
|
|
This file is part of the GNU MPFR Library.
|
|
|
|
The GNU MPFR Library 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 3 of the License, or (at your
|
|
option) any later version.
|
|
|
|
The GNU MPFR Library 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 Lesser General Public License
|
|
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
|
|
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
|
|
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
|
|
|
|
#define MPFR_NEED_LONGLONG_H
|
|
#include "mpfr-impl.h"
|
|
|
|
static int mpfr_yn_asympt (mpfr_ptr, long, mpfr_srcptr, mpfr_rnd_t);
|
|
|
|
int
|
|
mpfr_y0 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
|
|
{
|
|
return mpfr_yn (res, 0, z, r);
|
|
}
|
|
|
|
int
|
|
mpfr_y1 (mpfr_ptr res, mpfr_srcptr z, mpfr_rnd_t r)
|
|
{
|
|
return mpfr_yn (res, 1, z, r);
|
|
}
|
|
|
|
/* compute in s an approximation of S1 = sum((n-k)!/k!*y^k,k=0..n)
|
|
return e >= 0 the exponent difference between the maximal value of |s|
|
|
during the for loop and the final value of |s|.
|
|
*/
|
|
static mpfr_exp_t
|
|
mpfr_yn_s1 (mpfr_ptr s, mpfr_srcptr y, unsigned long n)
|
|
{
|
|
unsigned long k;
|
|
mpz_t f;
|
|
mpfr_exp_t e, emax;
|
|
|
|
mpz_init_set_ui (f, 1);
|
|
/* we compute n!*S1 = sum(a[k]*y^k,k=0..n) where a[k] = n!*(n-k)!/k!,
|
|
a[0] = (n!)^2, a[1] = n!*(n-1)!, ..., a[n-1] = n, a[n] = 1 */
|
|
mpfr_set_ui (s, 1, MPFR_RNDN); /* a[n] */
|
|
emax = MPFR_EXP(s);
|
|
for (k = n; k-- > 0;)
|
|
{
|
|
/* a[k]/a[k+1] = (n-k)!/k!/(n-(k+1))!*(k+1)! = (k+1)*(n-k) */
|
|
mpfr_mul (s, s, y, MPFR_RNDN);
|
|
mpz_mul_ui (f, f, n - k);
|
|
mpz_mul_ui (f, f, k + 1);
|
|
/* invariant: f = a[k] */
|
|
mpfr_add_z (s, s, f, MPFR_RNDN);
|
|
e = MPFR_EXP(s);
|
|
if (e > emax)
|
|
emax = e;
|
|
}
|
|
/* now we have f = (n!)^2 */
|
|
mpz_sqrt (f, f);
|
|
mpfr_div_z (s, s, f, MPFR_RNDN);
|
|
mpz_clear (f);
|
|
return emax - MPFR_EXP(s);
|
|
}
|
|
|
|
/* compute in s an approximation of
|
|
S3 = c*sum((h(k)+h(n+k))*y^k/k!/(n+k)!,k=0..infinity)
|
|
where h(k) = 1 + 1/2 + ... + 1/k
|
|
k=0: h(n)
|
|
k=1: 1+h(n+1)
|
|
k=2: 3/2+h(n+2)
|
|
Returns e such that the error is bounded by 2^e ulp(s).
|
|
*/
|
|
static mpfr_exp_t
|
|
mpfr_yn_s3 (mpfr_ptr s, mpfr_srcptr y, mpfr_srcptr c, unsigned long n)
|
|
{
|
|
unsigned long k, zz;
|
|
mpfr_t t, u;
|
|
mpz_t p, q; /* p/q will store h(k)+h(n+k) */
|
|
mpfr_exp_t exps, expU;
|
|
|
|
zz = mpfr_get_ui (y, MPFR_RNDU); /* y = z^2/4 */
|
|
MPFR_ASSERTN (zz < ULONG_MAX - 2);
|
|
zz += 2; /* z^2 <= 2^zz */
|
|
mpz_init_set_ui (p, 0);
|
|
mpz_init_set_ui (q, 1);
|
|
/* initialize p/q to h(n) */
|
|
for (k = 1; k <= n; k++)
|
|
{
|
|
/* p/q + 1/k = (k*p+q)/(q*k) */
|
|
mpz_mul_ui (p, p, k);
|
|
mpz_add (p, p, q);
|
|
mpz_mul_ui (q, q, k);
|
|
}
|
|
mpfr_init2 (t, MPFR_PREC(s));
|
|
mpfr_init2 (u, MPFR_PREC(s));
|
|
mpfr_fac_ui (t, n, MPFR_RNDN);
|
|
mpfr_div (t, c, t, MPFR_RNDN); /* c/n! */
|
|
mpfr_mul_z (u, t, p, MPFR_RNDN);
|
|
mpfr_div_z (s, u, q, MPFR_RNDN);
|
|
exps = MPFR_EXP (s);
|
|
expU = exps;
|
|
for (k = 1; ;k ++)
|
|
{
|
|
/* update t */
|
|
mpfr_mul (t, t, y, MPFR_RNDN);
|
|
mpfr_div_ui (t, t, k, MPFR_RNDN);
|
|
mpfr_div_ui (t, t, n + k, MPFR_RNDN);
|
|
/* update p/q:
|
|
p/q + 1/k + 1/(n+k) = [p*k*(n+k) + q*(n+k) + q*k]/(q*k*(n+k)) */
|
|
mpz_mul_ui (p, p, k);
|
|
mpz_mul_ui (p, p, n + k);
|
|
mpz_addmul_ui (p, q, n + 2 * k);
|
|
mpz_mul_ui (q, q, k);
|
|
mpz_mul_ui (q, q, n + k);
|
|
mpfr_mul_z (u, t, p, MPFR_RNDN);
|
|
mpfr_div_z (u, u, q, MPFR_RNDN);
|
|
exps = MPFR_EXP (u);
|
|
if (exps > expU)
|
|
expU = exps;
|
|
mpfr_add (s, s, u, MPFR_RNDN);
|
|
exps = MPFR_EXP (s);
|
|
if (exps > expU)
|
|
expU = exps;
|
|
if (MPFR_EXP (u) + (mpfr_exp_t) MPFR_PREC (u) < MPFR_EXP (s) &&
|
|
zz / (2 * k) < k + n)
|
|
break;
|
|
}
|
|
mpfr_clear (t);
|
|
mpfr_clear (u);
|
|
mpz_clear (p);
|
|
mpz_clear (q);
|
|
exps = expU - MPFR_EXP (s);
|
|
/* the error is bounded by (6k^2+33/2k+11) 2^exps ulps
|
|
<= 8*(k+2)^2 2^exps ulps */
|
|
return 3 + 2 * MPFR_INT_CEIL_LOG2(k + 2) + exps;
|
|
}
|
|
|
|
int
|
|
mpfr_yn (mpfr_ptr res, long n, mpfr_srcptr z, mpfr_rnd_t r)
|
|
{
|
|
int inex;
|
|
unsigned long absn;
|
|
|
|
MPFR_LOG_FUNC (("x[%#R]=%R n=%d rnd=%d", z, z, n, r),
|
|
("y[%#R]=%R", res, res));
|
|
|
|
absn = SAFE_ABS (unsigned long, n);
|
|
|
|
if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (z)))
|
|
{
|
|
if (MPFR_IS_NAN (z))
|
|
{
|
|
MPFR_SET_NAN (res); /* y(n,NaN) = NaN */
|
|
MPFR_RET_NAN;
|
|
}
|
|
/* y(n,z) tends to zero when z goes to +Inf, oscillating around
|
|
0. We choose to return +0 in that case. */
|
|
else if (MPFR_IS_INF (z))
|
|
{
|
|
if (MPFR_SIGN(z) > 0)
|
|
return mpfr_set_ui (res, 0, r);
|
|
else /* y(n,-Inf) = NaN */
|
|
{
|
|
MPFR_SET_NAN (res);
|
|
MPFR_RET_NAN;
|
|
}
|
|
}
|
|
else /* y(n,z) tends to -Inf for n >= 0 or n even, to +Inf otherwise,
|
|
when z goes to zero */
|
|
{
|
|
MPFR_SET_INF(res);
|
|
if (n >= 0 || (n & 1) == 0)
|
|
MPFR_SET_NEG(res);
|
|
else
|
|
MPFR_SET_POS(res);
|
|
MPFR_RET(0);
|
|
}
|
|
}
|
|
|
|
/* for z < 0, y(n,z) is imaginary except when j(n,|z|) = 0, which we
|
|
assume does not happen for a rational z. */
|
|
if (MPFR_SIGN(z) < 0)
|
|
{
|
|
MPFR_SET_NAN (res);
|
|
MPFR_RET_NAN;
|
|
}
|
|
|
|
/* now z is not singular, and z > 0 */
|
|
|
|
/* Deal with tiny arguments. We have:
|
|
y0(z) = 2 log(z)/Pi + 2 (euler - log(2))/Pi + O(log(z)*z^2), more
|
|
precisely for 0 <= z <= 1/2, with g(z) = 2/Pi + 2(euler-log(2))/Pi/log(z),
|
|
g(z) - 0.41*z^2 < y0(z)/log(z) < g(z)
|
|
thus since log(z) is negative:
|
|
g(z)*log(z) < y0(z) < (g(z) - z^2/2)*log(z)
|
|
and since |g(z)| >= 0.63 for 0 <= z <= 1/2, the relative error on
|
|
y0(z)/log(z) is bounded by 0.41*z^2/0.63 <= 0.66*z^2.
|
|
Note: we use both the main term in log(z) and the constant term, because
|
|
otherwise the relative error would be only in 1/log(|log(z)|).
|
|
*/
|
|
if (n == 0 && MPFR_EXP(z) < - (mpfr_exp_t) (MPFR_PREC(res) / 2))
|
|
{
|
|
mpfr_t l, h, t, logz;
|
|
mpfr_prec_t prec;
|
|
int ok, inex2;
|
|
|
|
prec = MPFR_PREC(res) + 10;
|
|
mpfr_init2 (l, prec);
|
|
mpfr_init2 (h, prec);
|
|
mpfr_init2 (t, prec);
|
|
mpfr_init2 (logz, prec);
|
|
/* first enclose log(z) + euler - log(2) = log(z/2) + euler */
|
|
mpfr_log (logz, z, MPFR_RNDD); /* lower bound of log(z) */
|
|
mpfr_set (h, logz, MPFR_RNDU); /* exact */
|
|
mpfr_nextabove (h); /* upper bound of log(z) */
|
|
mpfr_const_euler (t, MPFR_RNDD); /* lower bound of euler */
|
|
mpfr_add (l, logz, t, MPFR_RNDD); /* lower bound of log(z) + euler */
|
|
mpfr_nextabove (t); /* upper bound of euler */
|
|
mpfr_add (h, h, t, MPFR_RNDU); /* upper bound of log(z) + euler */
|
|
mpfr_const_log2 (t, MPFR_RNDU); /* upper bound of log(2) */
|
|
mpfr_sub (l, l, t, MPFR_RNDD); /* lower bound of log(z/2) + euler */
|
|
mpfr_nextbelow (t); /* lower bound of log(2) */
|
|
mpfr_sub (h, h, t, MPFR_RNDU); /* upper bound of log(z/2) + euler */
|
|
mpfr_const_pi (t, MPFR_RNDU); /* upper bound of Pi */
|
|
mpfr_div (l, l, t, MPFR_RNDD); /* lower bound of (log(z/2)+euler)/Pi */
|
|
mpfr_nextbelow (t); /* lower bound of Pi */
|
|
mpfr_div (h, h, t, MPFR_RNDD); /* upper bound of (log(z/2)+euler)/Pi */
|
|
mpfr_mul_2ui (l, l, 1, MPFR_RNDD); /* lower bound on g(z)*log(z) */
|
|
mpfr_mul_2ui (h, h, 1, MPFR_RNDU); /* upper bound on g(z)*log(z) */
|
|
/* we now have l <= g(z)*log(z) <= h, and we need to add -z^2/2*log(z)
|
|
to h */
|
|
mpfr_mul (t, z, z, MPFR_RNDU); /* upper bound on z^2 */
|
|
/* since logz is negative, a lower bound corresponds to an upper bound
|
|
for its absolute value */
|
|
mpfr_neg (t, t, MPFR_RNDD);
|
|
mpfr_div_2ui (t, t, 1, MPFR_RNDD);
|
|
mpfr_mul (t, t, logz, MPFR_RNDU); /* upper bound on z^2/2*log(z) */
|
|
/* an underflow may happen in the above instructions, clear flag */
|
|
mpfr_clear_underflow ();
|
|
mpfr_add (h, h, t, MPFR_RNDU);
|
|
inex = mpfr_prec_round (l, MPFR_PREC(res), r);
|
|
inex2 = mpfr_prec_round (h, MPFR_PREC(res), r);
|
|
/* we need h=l and inex=inex2 */
|
|
ok = (inex == inex2) && (mpfr_cmp (l, h) == 0);
|
|
if (ok)
|
|
mpfr_set (res, h, r); /* exact */
|
|
mpfr_clear (l);
|
|
mpfr_clear (h);
|
|
mpfr_clear (t);
|
|
mpfr_clear (logz);
|
|
if (ok)
|
|
return inex;
|
|
}
|
|
|
|
/* small argument check for y1(z) = -2/Pi/z + O(log(z)):
|
|
for 0 <= z <= 1, |y1(z) + 2/Pi/z| <= 0.25 */
|
|
if (n == 1 && MPFR_EXP(z) + 1 < - (mpfr_exp_t) MPFR_PREC(res))
|
|
{
|
|
mpfr_t y;
|
|
mpfr_prec_t prec;
|
|
mpfr_exp_t err1;
|
|
int ok;
|
|
MPFR_BLOCK_DECL (flags);
|
|
|
|
/* since 2/Pi > 0.5, and |y1(z)| >= |2/Pi/z|, if z <= 2^(-emax-1),
|
|
then |y1(z)| > 2^emax */
|
|
prec = MPFR_PREC(res) + 10;
|
|
mpfr_init2 (y, prec);
|
|
mpfr_const_pi (y, MPFR_RNDU); /* Pi*(1+u)^2, where here and below u
|
|
represents a quantity <= 1/2^prec */
|
|
mpfr_mul (y, y, z, MPFR_RNDU); /* Pi*z * (1+u)^4, upper bound */
|
|
MPFR_BLOCK (flags, mpfr_ui_div (y, 2, y, MPFR_RNDZ));
|
|
/* 2/Pi/z * (1+u)^6, lower bound, with possible overflow */
|
|
if (MPFR_OVERFLOW (flags))
|
|
{
|
|
mpfr_clear (y);
|
|
return mpfr_overflow (res, r, -1);
|
|
}
|
|
mpfr_neg (y, y, MPFR_RNDN);
|
|
/* (1+u)^6 can be written 1+7u [for another value of u], thus the
|
|
error on 2/Pi/z is less than 7ulp(y). The truncation error is less
|
|
than 1/4, thus if ulp(y)>=1/4, the total error is less than 8ulp(y),
|
|
otherwise it is less than 1/4+7/8 <= 2. */
|
|
if (MPFR_EXP(y) + 2 >= MPFR_PREC(y)) /* ulp(y) >= 1/4 */
|
|
err1 = 3;
|
|
else /* ulp(y) <= 1/8 */
|
|
err1 = (mpfr_exp_t) MPFR_PREC(y) - MPFR_EXP(y) + 1;
|
|
ok = MPFR_CAN_ROUND (y, prec - err1, MPFR_PREC(res), r);
|
|
if (ok)
|
|
inex = mpfr_set (res, y, r);
|
|
mpfr_clear (y);
|
|
if (ok)
|
|
return inex;
|
|
}
|
|
|
|
/* we can use the asymptotic expansion as soon as z > p log(2)/2,
|
|
but to get some margin we use it for z > p/2 */
|
|
if (mpfr_cmp_ui (z, MPFR_PREC(res) / 2 + 3) > 0)
|
|
{
|
|
inex = mpfr_yn_asympt (res, n, z, r);
|
|
if (inex != 0)
|
|
return inex;
|
|
}
|
|
|
|
/* General case */
|
|
{
|
|
mpfr_prec_t prec;
|
|
mpfr_exp_t err1, err2, err3;
|
|
mpfr_t y, s1, s2, s3;
|
|
MPFR_ZIV_DECL (loop);
|
|
|
|
mpfr_init (y);
|
|
mpfr_init (s1);
|
|
mpfr_init (s2);
|
|
mpfr_init (s3);
|
|
|
|
prec = MPFR_PREC(res) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (res)) + 13;
|
|
MPFR_ZIV_INIT (loop, prec);
|
|
for (;;)
|
|
{
|
|
mpfr_set_prec (y, prec);
|
|
mpfr_set_prec (s1, prec);
|
|
mpfr_set_prec (s2, prec);
|
|
mpfr_set_prec (s3, prec);
|
|
|
|
mpfr_mul (y, z, z, MPFR_RNDN);
|
|
mpfr_div_2ui (y, y, 2, MPFR_RNDN); /* z^2/4 */
|
|
|
|
/* store (z/2)^n temporarily in s2 */
|
|
mpfr_pow_ui (s2, z, absn, MPFR_RNDN);
|
|
mpfr_div_2si (s2, s2, absn, MPFR_RNDN);
|
|
|
|
/* compute S1 * (z/2)^(-n) */
|
|
if (n == 0)
|
|
{
|
|
mpfr_set_ui (s1, 0, MPFR_RNDN);
|
|
err1 = 0;
|
|
}
|
|
else
|
|
err1 = mpfr_yn_s1 (s1, y, absn - 1);
|
|
mpfr_div (s1, s1, s2, MPFR_RNDN); /* (z/2)^(-n) * S1 */
|
|
/* See algorithms.tex: the relative error on s1 is bounded by
|
|
(3n+3)*2^(e+1-prec). */
|
|
err1 = MPFR_INT_CEIL_LOG2 (3 * absn + 3) + err1 + 1;
|
|
/* rel_err(s1) <= 2^(err1-prec), thus err(s1) <= 2^err1 ulps */
|
|
|
|
/* compute (z/2)^n * S3 */
|
|
mpfr_neg (y, y, MPFR_RNDN); /* -z^2/4 */
|
|
err3 = mpfr_yn_s3 (s3, y, s2, absn); /* (z/2)^n * S3 */
|
|
/* the error on s3 is bounded by 2^err3 ulps */
|
|
|
|
/* add s1+s3 */
|
|
err1 += MPFR_EXP(s1);
|
|
mpfr_add (s1, s1, s3, MPFR_RNDN);
|
|
/* the error is bounded by 1/2 + 2^err1*2^(- EXP(s1))
|
|
+ 2^err3*2^(EXP(s3) - EXP(s1)) */
|
|
err3 += MPFR_EXP(s3);
|
|
err1 = (err3 > err1) ? err3 + 1 : err1 + 1;
|
|
err1 -= MPFR_EXP(s1);
|
|
err1 = (err1 >= 0) ? err1 + 1 : 1;
|
|
/* now the error on s1 is bounded by 2^err1*ulp(s1) */
|
|
|
|
/* compute S2 */
|
|
mpfr_div_2ui (s2, z, 1, MPFR_RNDN); /* z/2 */
|
|
mpfr_log (s2, s2, MPFR_RNDN); /* log(z/2) */
|
|
mpfr_const_euler (s3, MPFR_RNDN);
|
|
err2 = MPFR_EXP(s2) > MPFR_EXP(s3) ? MPFR_EXP(s2) : MPFR_EXP(s3);
|
|
mpfr_add (s2, s2, s3, MPFR_RNDN); /* log(z/2) + gamma */
|
|
err2 -= MPFR_EXP(s2);
|
|
mpfr_mul_2ui (s2, s2, 1, MPFR_RNDN); /* 2*(log(z/2) + gamma) */
|
|
mpfr_jn (s3, absn, z, MPFR_RNDN); /* Jn(z) */
|
|
mpfr_mul (s2, s2, s3, MPFR_RNDN); /* 2*(log(z/2) + gamma)*Jn(z) */
|
|
err2 += 4; /* the error on s2 is bounded by 2^err2 ulps, see
|
|
algorithms.tex */
|
|
|
|
/* add all three sums */
|
|
err1 += MPFR_EXP(s1); /* the error on s1 is bounded by 2^err1 */
|
|
err2 += MPFR_EXP(s2); /* the error on s2 is bounded by 2^err2 */
|
|
mpfr_sub (s2, s2, s1, MPFR_RNDN); /* s2 - (s1+s3) */
|
|
err2 = (err1 > err2) ? err1 + 1 : err2 + 1;
|
|
err2 -= MPFR_EXP(s2);
|
|
err2 = (err2 >= 0) ? err2 + 1 : 1;
|
|
/* now the error on s2 is bounded by 2^err2*ulp(s2) */
|
|
mpfr_const_pi (y, MPFR_RNDN); /* error bounded by 1 ulp */
|
|
mpfr_div (s2, s2, y, MPFR_RNDN); /* error bounded by
|
|
2^(err2+1)*ulp(s2) */
|
|
err2 ++;
|
|
|
|
if (MPFR_LIKELY (MPFR_CAN_ROUND (s2, prec - err2, MPFR_PREC(res), r)))
|
|
break;
|
|
MPFR_ZIV_NEXT (loop, prec);
|
|
}
|
|
MPFR_ZIV_FREE (loop);
|
|
|
|
inex = (n >= 0 || (n & 1) == 0)
|
|
? mpfr_set (res, s2, r)
|
|
: mpfr_neg (res, s2, r);
|
|
|
|
mpfr_clear (y);
|
|
mpfr_clear (s1);
|
|
mpfr_clear (s2);
|
|
mpfr_clear (s3);
|
|
}
|
|
|
|
return inex;
|
|
}
|
|
|
|
#define MPFR_YN
|
|
#include "jyn_asympt.c"
|