Commit c1266acb authored by Julien Puydt's avatar Julien Puydt

New upstream version 2.13.0

parent 8c3c7c59
......@@ -18,9 +18,12 @@ compiler:
script:
- if [[ "${TRAVIS_OS_NAME}" == "osx" ]] && [[ "${CC}" == "gcc" ]]; then
brew update;
brew install gcc48;
export CC=gcc-4.8;
export CXX=g++-4.8;
brew install gcc;
brew link --overwrite gcc;
export CC=gcc;
export CXX=g++;
gcc --version;
g++ --version;
fi
- if [[ "${TRAVIS_OS_NAME}" == "osx" ]]; then
export ARB_TEST_MULTIPLIER=0.1;
......
......@@ -52,7 +52,7 @@ $(BUILD_DIR)/$(MOD_DIR)_%.o: %.c
$(QUIET_CC) $(CC) $(CFLAGS) $(INCS) -c $< -o $@ -MMD -MP -MF "$(BUILD_DIR)/$(MOD_DIR)_$*.d" -MT "$(BUILD_DIR)/$(MOD_DIR)_$*.d" -MT "$@"
$(MOD_LOBJ): $(LOBJS)
$(QUIET_CC) $(CC) $(ABI_FLAG) -Wl,-r $^ -o $@ -nostdlib
$(QUIET_CC) $(CC) $(ABI_FLAG) -r $^ -o $@ -nostdlib
-include $(LOBJS:.lo=.d)
......
......@@ -161,3 +161,6 @@ A separate wrapper of transcendental functions for use with the
C99 `complex double` type is available
(<https://github.com/fredrik-johansson/arbcmath>).
Other third-party wrappers include:
* Java wrapper using JNA: https://github.com/crowlogic/arb/
* Another Julia interface: https://github.com/JuliaArbTypes/ArbFloats.jl
......@@ -404,6 +404,20 @@ void acb_sgn(acb_t res, const acb_t z, slong prec);
void acb_csgn(arb_t res, const acb_t z);
void acb_real_abs(acb_t res, const acb_t z, int analytic, slong prec);
void acb_real_sgn(acb_t res, const acb_t z, int analytic, slong prec);
void acb_real_heaviside(acb_t res, const acb_t z, int analytic, slong prec);
void acb_real_floor(acb_t res, const acb_t z, int analytic, slong prec);
void acb_real_ceil(acb_t res, const acb_t z, int analytic, slong prec);
void acb_real_max(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec);
void acb_real_min(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec);
void acb_real_sqrtpos(acb_t res, const acb_t z, int analytic, slong prec);
void acb_sqrt_analytic(acb_t res, const acb_t z, int analytic, slong prec);
void acb_rsqrt_analytic(acb_t res, const acb_t z, int analytic, slong prec);
void acb_log_analytic(acb_t res, const acb_t z, int analytic, slong prec);
void acb_pow_analytic(acb_t res, const acb_t z, const acb_t w, int analytic, slong prec);
ACB_INLINE void
acb_add(acb_t z, const acb_t x, const acb_t y, slong prec)
{
......
......@@ -16,24 +16,38 @@ acb_cos(acb_t r, const acb_t z, slong prec)
{
#define a acb_realref(z)
#define b acb_imagref(z)
arb_t sa, ca, sb, cb;
arb_init(sa);
arb_init(ca);
arb_init(sb);
arb_init(cb);
arb_sin_cos(sa, ca, a, prec);
arb_sinh_cosh(sb, cb, b, prec);
arb_mul(acb_realref(r), ca, cb, prec);
arb_mul(acb_imagref(r), sa, sb, prec);
arb_neg(acb_imagref(r), acb_imagref(r));
arb_clear(sa);
arb_clear(ca);
arb_clear(sb);
arb_clear(cb);
if (arb_is_zero(b))
{
arb_cos(acb_realref(r), a, prec);
arb_zero(acb_imagref(r));
}
else if (arb_is_zero(a))
{
arb_cosh(acb_realref(r), b, prec);
arb_zero(acb_imagref(r));
}
else
{
arb_t sa, ca, sb, cb;
arb_init(sa);
arb_init(ca);
arb_init(sb);
arb_init(cb);
arb_sin_cos(sa, ca, a, prec);
arb_sinh_cosh(sb, cb, b, prec);
arb_mul(acb_realref(r), ca, cb, prec);
arb_mul(acb_imagref(r), sa, sb, prec);
arb_neg(acb_imagref(r), acb_imagref(r));
arb_clear(sa);
arb_clear(ca);
arb_clear(sb);
arb_clear(cb);
}
#undef a
#undef b
}
......@@ -16,28 +16,44 @@ acb_cos_pi(acb_t r, const acb_t z, slong prec)
{
#define a acb_realref(z)
#define b acb_imagref(z)
arb_t sa, ca, sb, cb;
arb_init(sa);
arb_init(ca);
arb_init(sb);
arb_init(cb);
arb_sin_cos_pi(sa, ca, a, prec);
arb_const_pi(cb, prec);
arb_mul(cb, cb, b, prec);
arb_sinh_cosh(sb, cb, cb, prec);
arb_mul(acb_realref(r), ca, cb, prec);
arb_mul(acb_imagref(r), sa, sb, prec);
arb_neg(acb_imagref(r), acb_imagref(r));
arb_clear(sa);
arb_clear(ca);
arb_clear(sb);
arb_clear(cb);
if (arb_is_zero(b))
{
arb_cos_pi(acb_realref(r), a, prec);
arb_zero(acb_imagref(r));
}
else if (arb_is_zero(a))
{
arb_t t;
arb_init(t);
arb_const_pi(t, prec);
arb_mul(t, t, b, prec);
arb_cosh(acb_realref(r), t, prec);
arb_zero(acb_imagref(r));
arb_clear(t);
}
else
{
arb_t sa, ca, sb, cb;
arb_init(sa);
arb_init(ca);
arb_init(sb);
arb_init(cb);
arb_sin_cos_pi(sa, ca, a, prec);
arb_const_pi(cb, prec);
arb_mul(cb, cb, b, prec);
arb_sinh_cosh(sb, cb, cb, prec);
arb_mul(acb_realref(r), ca, cb, prec);
arb_mul(acb_imagref(r), sa, sb, prec);
arb_neg(acb_imagref(r), acb_imagref(r));
arb_clear(sa);
arb_clear(ca);
arb_clear(sb);
arb_clear(cb);
}
#undef a
#undef b
}
......
......@@ -127,7 +127,30 @@ _acb_gamma(acb_t y, const acb_t x, slong prec, int inverse)
wp = FLINT_MAX(wp, 2);
wp = wp + FLINT_BIT_COUNT(wp);
acb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp);
if (acc < 3) /* try to avoid divisions blowing up */
{
if (arf_cmp_d(arb_midref(acb_realref(x)), -0.5) < 0)
{
reflect = 1;
r = 0;
}
else if (arf_cmp_si(arb_midref(acb_realref(x)), 1) < 0)
{
reflect = 0;
r = 1;
}
else
{
reflect = 0;
r = 0;
}
n = 1;
}
else
{
acb_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp);
}
acb_init(t);
acb_init(u);
......@@ -135,7 +158,6 @@ _acb_gamma(acb_t y, const acb_t x, slong prec, int inverse)
if (reflect)
{
/* gamma(x) = (rf(1-x, r) * pi) / (gamma(1-x+r) sin(pi x)) */
acb_sub_ui(t, x, 1, wp);
acb_neg(t, t);
acb_rising_ui_rec(u, t, r, wp);
......@@ -143,23 +165,47 @@ _acb_gamma(acb_t y, const acb_t x, slong prec, int inverse)
acb_mul_arb(u, u, acb_realref(v), wp);
acb_add_ui(t, t, r, wp);
acb_gamma_stirling_eval(v, t, n, 0, wp);
acb_exp(v, v, wp);
acb_sin_pi(t, x, wp);
acb_mul(v, v, t, wp);
if (inverse)
{
/* rgamma(x) = gamma(1-x+r) sin(pi x) / ((rf(1-x, r) * pi) */
acb_exp(v, v, wp);
acb_sin_pi(t, x, wp);
acb_mul(v, v, t, wp);
acb_mul(y, u, v, wp);
acb_div(y, v, u, prec);
}
else
{
/* gamma(x) = (rf(1-x, r) * pi) rgamma(1-x+r) csc(pi x) */
acb_neg(v, v);
acb_exp(v, v, wp);
acb_sin_pi(t, x, wp); /* todo: write a csc_pi function */
acb_div(v, v, t, wp);
acb_mul(y, v, u, prec);
}
}
else
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
acb_add_ui(t, x, r, wp);
acb_gamma_stirling_eval(u, t, n, 0, wp);
acb_exp(u, u, prec);
acb_rising_ui_rec(v, x, r, wp);
}
if (inverse)
acb_div(y, v, u, prec);
else
acb_div(y, u, v, prec);
if (inverse)
{
/* rgamma(x) = rf(x,r) rgamma(x+r) */
acb_neg(u, u);
acb_exp(u, u, prec);
acb_rising_ui_rec(v, x, r, wp);
acb_mul(y, v, u, prec);
}
else
{
/* gamma(x) = gamma(x+r) / rf(x,r) */
acb_exp(u, u, prec);
acb_rising_ui_rec(v, x, r, wp);
acb_div(y, u, v, prec);
}
}
acb_clear(t);
acb_clear(u);
......
......@@ -14,10 +14,28 @@
void
acb_get_mag_lower(mag_t z, const acb_t x)
{
arf_t t;
arf_init(t);
acb_get_abs_lbound_arf(t, x, MAG_BITS);
arf_get_mag_lower(z, t);
arf_clear(t);
if (arb_is_zero(acb_imagref(x)))
{
arb_get_mag_lower(z, acb_realref(x));
}
else if (arb_is_zero(acb_realref(x)))
{
arb_get_mag_lower(z, acb_imagref(x));
}
else
{
mag_t t;
mag_init(t);
arb_get_mag_lower(t, acb_realref(x));
arb_get_mag_lower(z, acb_imagref(x));
mag_mul_lower(t, t, t);
mag_mul_lower(z, z, z);
mag_add_lower(z, z, t);
mag_sqrt_lower(z, z);
mag_clear(t);
}
}
......@@ -11,30 +11,6 @@
#include "acb.h"
static int
close_to_one(const acb_t z)
{
mp_limb_t top;
if (arf_abs_bound_lt_2exp_si(arb_midref(acb_imagref(z))) > -3)
return 0;
if (ARF_EXP(arb_midref(acb_realref(z))) == 0)
{
ARF_GET_TOP_LIMB(top, arb_midref(acb_realref(z)));
return (top >> (FLINT_BITS - 4)) == 15;
}
else if (ARF_EXP(arb_midref(acb_realref(z))) == 1)
{
ARF_GET_TOP_LIMB(top, arb_midref(acb_realref(z)));
return (top >> (FLINT_BITS - 4)) == 8;
}
return 0;
}
void
acb_log(acb_t r, const acb_t z, slong prec)
{
......@@ -82,48 +58,42 @@ acb_log(acb_t r, const acb_t z, slong prec)
}
else
{
arb_t t, u;
arb_init(t);
arb_init(u);
if (close_to_one(z))
if (r != z)
{
arb_sub_ui(u, a, 1, prec + 8);
arb_mul(t, u, u, prec + 8);
arb_addmul(t, b, b, prec + 8);
arb_mul_2exp_si(u, u, 1);
arb_add(t, t, u, prec + 8);
arb_log1p(t, t, prec);
arb_mul_2exp_si(t, t, -1);
arb_log_hypot(acb_realref(r), a, b, prec);
if (arb_is_finite(acb_realref(r)))
arb_atan2(acb_imagref(r), b, a, prec);
else
arb_indeterminate(acb_imagref(r));
}
else
{
arb_mul(t, a, a, prec + 8);
arb_addmul(t, b, b, prec + 8);
if (arb_contains_zero(t) || arf_sgn(arb_midref(t)) < 0)
arb_zero_pm_inf(t);
arb_t t;
arb_init(t);
arb_log_hypot(t, a, b, prec);
if (arb_is_finite(t))
arb_atan2(acb_imagref(r), b, a, prec);
else
arb_log(t, t, prec);
arb_mul_2exp_si(t, t, -1);
arb_indeterminate(acb_imagref(r));
arb_swap(acb_realref(r), t);
arb_clear(t);
}
acb_arg(u, z, prec);
arb_swap(acb_realref(r), t);
arb_swap(acb_imagref(r), u);
arb_clear(t);
arb_clear(u);
}
if (!acb_is_finite(r))
acb_indeterminate(r);
#undef a
#undef b
}
void
acb_log_analytic(acb_ptr res, const acb_t z, int analytic, slong prec)
{
if (analytic && arb_contains_zero(acb_imagref(z)) &&
!arb_is_positive(acb_realref(z)))
{
acb_indeterminate(res);
}
else
{
acb_log(res, z, prec);
}
}
......@@ -233,3 +233,17 @@ acb_pow(acb_t z, const acb_t x, const acb_t y, slong prec)
}
}
void
acb_pow_analytic(acb_ptr res, const acb_t z, const acb_t w, int analytic, slong prec)
{
if (analytic && !acb_is_int(w) && arb_contains_zero(acb_imagref(z)) &&
!arb_is_positive(acb_realref(z)))
{
acb_indeterminate(res);
}
else
{
acb_pow(res, z, w, prec);
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_abs(acb_t res, const acb_t z, int analytic, slong prec)
{
if (!acb_is_finite(z) || (analytic && arb_contains_zero(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
if (arb_is_nonnegative(acb_realref(z)))
{
acb_set_round(res, z, prec);
}
else if (arb_is_negative(acb_realref(z)))
{
acb_neg_round(res, z, prec);
}
else
{
acb_t t;
acb_init(t);
acb_neg(t, res);
acb_union(res, z, t, prec);
acb_clear(t);
}
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_ceil(acb_t res, const acb_t z, int analytic, slong prec)
{
if (!acb_is_finite(z) || (analytic && arb_contains_int(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
arb_ceil(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_floor(acb_t res, const acb_t z, int analytic, slong prec)
{
if (!acb_is_finite(z) || (analytic && arb_contains_int(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
arb_floor(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_heaviside(acb_t res, const acb_t z, int analytic, slong prec)
{
acb_real_sgn(res, z, analytic, prec);
if (acb_is_finite(res))
{
acb_add_ui(res, res, 1, prec);
acb_mul_2exp_si(res, res, -1);
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_max(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec)
{
arb_t t;
if (!acb_is_finite(x) || !acb_is_finite(y))
{
acb_indeterminate(res);
return;
}
arb_init(t);
arb_sub(t, acb_realref(x), acb_realref(y), prec);
if (arb_is_positive(t))
acb_set_round(res, x, prec);
else if (arb_is_negative(t))
acb_set_round(res, y, prec);
else if (!analytic)
acb_union(res, x, y, prec);
else
acb_indeterminate(res);
arb_clear(t);
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_min(acb_t res, const acb_t x, const acb_t y, int analytic, slong prec)
{
arb_t t;
if (!acb_is_finite(x) || !acb_is_finite(y))
{
acb_indeterminate(res);
return;
}
arb_init(t);
arb_sub(t, acb_realref(x), acb_realref(y), prec);
if (arb_is_positive(t))
acb_set_round(res, y, prec);
else if (arb_is_negative(t))
acb_set_round(res, x, prec);
else if (!analytic)
acb_union(res, x, y, prec);
else
acb_indeterminate(res);
arb_clear(t);
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_sgn(acb_t res, const acb_t z, int analytic, slong prec)
{
if (!acb_is_finite(z) || (analytic && arb_contains_zero(acb_realref(z))))
{
acb_indeterminate(res);
}
else
{
acb_csgn(acb_realref(res), z);
arb_zero(acb_imagref(res));
}
}
/*
Copyright (C) 2018 Fredrik Johansson
This file is part of Arb.
Arb is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See <http://www.gnu.org/licenses/>.
*/
#include "acb.h"
void
acb_real_sqrtpos(acb_t res, const acb_t z, int analytic, slong prec)
{
if (arb_is_zero(acb_imagref(z)) && !analytic)
{
arb_sqrtpos(acb_realref(res), acb_realref(z), prec);
arb_zero(acb_imagref(res));
}
else if (arb_is_positive(acb_realref(z)) || !arb_contains_zero(acb_imagref(z)))
{
acb_sqrt(res, z, prec);
}
else
{
acb_indeterminate(res);
}
}
......@@ -11,12 +11,98 @@
#include "acb.h"