hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix log2 estimates

Browse source
This commit is contained in:
p15-git-acc 2019-02-18 11:08:04 -06:00
parent 9d33e8305f
commit baa1514ab3
3 changed files with 63 additions and 43 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

96
acb_dirichlet/backlund_s.c
View file

@ -11,61 +11,77 @@
#include "acb_dirichlet.h"
/* S(t) = N(t) - theta(t)/pi - 1 */
static void
_backlund_s(arb_t res, const arb_t t, slong prec)
{
arb_t N;
acb_t z;
arb_init(N);
acb_init(z);
acb_dirichlet_zeta_nzeros(N, t, prec);
acb_set_arb(z, t);
acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
arb_const_pi(acb_imagref(z), prec);
arb_div(acb_realref(z), acb_realref(z), acb_imagref(z), prec);
arb_sub(res, N, acb_realref(z), prec);
arb_sub_ui(res, res, 1, prec);
arb_clear(N);
acb_clear(z);
}
void
acb_dirichlet_backlund_s(arb_t res, const arb_t t, slong prec)
{
mag_t m;
mag_t m, b;
mag_init(m);
mag_init(b);
arb_get_mag(m, t);
if (!arb_is_nonnegative(t))
{
arb_indeterminate(res);
return;
}
mag_init(m);
/*
* If the error radius of t is greater than 8/log(t) then use a fast bound.
* The threshold doesn't matter for correctness, but note that the average
* spacing of zeros is 2*pi/log(t) asymptotically.
* Otherwise if the radius is smaller use S(t) = N(t) - theta(t)/pi - 1.
*/
arb_get_mag(m, t);
mag_log(m, m);
mag_mul_2exp_si(m, m, -3);
mag_inv(m, m);
if (mag_cmp(arb_radref(t), m) > 0)
else if (mag_cmp_2exp_si(m, 6) < 0) /* If t is small then N(t) is fast. */
{
arb_zero(res);
acb_dirichlet_backlund_s_bound(arb_radref(res), t);
_backlund_s(res, t, prec);
}
else
{
arb_t N;
acb_t z;
arb_init(N);
acb_init(z);
/*
* The difference (N(t) - theta(t)/pi) loses precision from
* cancellation. Add extra bits to mitigate the loss.
* If the error radius of t is greater than 8/log(t) then use a fast
* bound. The threshold doesn't matter for correctness, but note that
* the average spacing of zeros is 2*pi/log(t) asymptotically.
* Otherwise if the radius is smaller use S(t) = N(t) - theta(t)/pi - 1.
*/
acb_set_arb(z, t);
acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
prec += acb_bits(z);
acb_dirichlet_zeta_nzeros(N, t, prec);
acb_set_arb(z, t);
acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
arb_const_pi(acb_imagref(z), prec);
arb_div(acb_realref(z), acb_realref(z), acb_imagref(z), prec);
arb_sub(res, N, acb_realref(z), prec);
arb_sub_ui(res, res, 1, prec);
arb_clear(N);
acb_clear(z);
mag_log(b, m);
mag_mul_2exp_si(b, b, -3);
mag_inv(b, b);
if (mag_cmp(arb_radref(t), b) > 0)
{
arb_zero(res);
acb_dirichlet_backlund_s_bound(arb_radref(res), t);
}
else
{
acb_t z;
acb_init(z);
/*
* The difference (N(t) - theta(t)/pi) loses precision from
* cancellation. Add extra bits to mitigate the loss.
*/
acb_set_arb(z, t);
acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
_backlund_s(res, t, prec + mag_get_d_log2_approx(m));
acb_clear(z);
}
}
mag_clear(m);
mag_clear(b);
}

8
acb_dirichlet/hardy_z_zero.c
View file

@ -113,6 +113,7 @@ void
acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)
{
arf_interval_t r, s;
mag_t m;
slong bits;
if (fmpz_cmp_si(n, 1) < 0)
@ -123,11 +124,13 @@ acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)
arf_interval_init(r);
arf_interval_init(s);
mag_init(m);
acb_dirichlet_isolate_hardy_z_zero(&r->a, &r->b, n);
bits = arf_bits(&r->b);
arb_set_interval_arf(res, &r->a, &r->b, bits + 8);
arf_get_mag(m, &r->b);
bits = FLINT_MAX(0, mag_get_d_log2_approx(m)) + 8;
arb_set_interval_arf(res, &r->a, &r->b, bits);
bits = arb_rel_accuracy_bits(res);
if (bits < prec)
@ -140,4 +143,5 @@ acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)
arf_interval_clear(r);
arf_interval_clear(s);
mag_clear(m);
}

2
acb_dirichlet/test/t-backlund_s.c
View file

@ -32,7 +32,7 @@ int main()
arb_init(b);
arb_init(c);
arb_randtest(a, state, 1 + n_randint(state, 1000), 3);
arb_randtest(a, state, 1 + n_randint(state, 1000), 4);
acb_dirichlet_backlund_s(b, a, prec1);