multithreaded refinement of zeta zeros; tweak example program

This commit is contained in:
Fredrik Johansson 2022-05-24 10:26:39 +02:00
parent 8e2436c7a1
commit 31578cdefe
3 changed files with 63 additions and 37 deletions

View file

@ -10,6 +10,7 @@
(at your option) any later version. See <http://www.gnu.org/licenses/>. (at your option) any later version. See <http://www.gnu.org/licenses/>.
*/ */
#include "flint/thread_support.h"
#include "acb_dirichlet.h" #include "acb_dirichlet.h"
#include "arb_calc.h" #include "arb_calc.h"
@ -1376,6 +1377,20 @@ _acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t)
exact_zeta_multi_nzeros(res, t, 1); exact_zeta_multi_nzeros(res, t, 1);
} }
typedef struct
{
arb_ptr res;
arf_interval_ptr p;
slong prec;
}
work_t;
static void
refinement_worker(slong i, work_t * work)
{
_acb_dirichlet_refine_hardy_z_zero(work->res + i, &(work->p[i].a), &(work->p[i].b), work->prec);
}
void void
acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec) acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
{ {
@ -1390,13 +1405,15 @@ acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
} }
else else
{ {
slong i; work_t work;
arf_interval_ptr p = _arf_interval_vec_init(len); arf_interval_ptr p = _arf_interval_vec_init(len);
acb_dirichlet_isolate_hardy_z_zeros(p, n, len); acb_dirichlet_isolate_hardy_z_zeros(p, n, len);
for (i = 0; i < len; i++)
{ work.res = res;
_acb_dirichlet_refine_hardy_z_zero(res + i, &p[i].a, &p[i].b, prec); work.p = p;
} work.prec = prec;
flint_parallel_do((do_func_t) refinement_worker, &work, len, -1, FLINT_PARALLEL_STRIDED);
_arf_interval_vec_clear(p, len); _arf_interval_vec_clear(p, len);
} }
} }

View file

@ -53,41 +53,37 @@ zeta_zeros.c
This program computes one or several consecutive zeros of the This program computes one or several consecutive zeros of the
Riemann zeta function on the critical line:: Riemann zeta function on the critical line::
> build/examples/zeta_zeros -num 10 -digits 30 > build/examples/zeta_zeros -n 1 -count 10 -digits 30
Computing 10 zeros starting with zero #1, with prec = 102 bits... 1 14.1347251417346937904572519836
cpu/wall(s): 0.01 0.009 2 21.0220396387715549926284795939
virt/peak/res/peak(MB): 21.30 21.30 7.32 7.32 3 25.0108575801456887632137909926
14.1347251417346937904572519836 4 30.4248761258595132103118975306
21.0220396387715549926284795939 5 32.9350615877391896906623689641
25.0108575801456887632137909926 6 37.5861781588256712572177634807
30.4248761258595132103118975306 7 40.9187190121474951873981269146
32.9350615877391896906623689641 8 43.3270732809149995194961221654
37.5861781588256712572177634807 9 48.0051508811671597279424727494
40.9187190121474951873981269146 10 49.7738324776723021819167846786
43.3270732809149995194961221654 cpu/wall(s): 0.01 0.01
48.0051508811671597279424727494 virt/peak/res/peak(MB): 21.28 21.28 7.29 7.29
49.7738324776723021819167846786
Five zeros starting with the millionth:: Five zeros starting with the millionth::
> build/examples/zeta_zeros -start 1000000 -num 5 -digits 20 > build/examples/zeta_zeros -n 1000000 -count 5 -digits 20
Computing 5 zeros starting with zero #1000000, with prec = 69 bits... 1000000 600269.67701244495552
cpu/wall(s): 0.029 0.029 1000001 600270.30109071169866
virt/peak/res/peak(MB): 21.16 21.16 5.82 5.82 1000002 600270.74787059436613
600269.67701244495552 1000003 600271.48637367364820
600270.30109071169866 1000004 600271.76148042593778
600270.74787059436613 cpu/wall(s): 0.03 0.03
600271.48637367364820 virt/peak/res/peak(MB): 21.41 21.41 7.41 7.41
600271.76148042593778
The program supports the following options:: The program supports the following options::
zeta_zeros [-start n] [-num n] [-threads t] [-digits d] [-suppress] [-gradual] [-stream] zeta_zeros [-n n] [-count n] [-prec n] [-digits n] [-threads n] [-platt] [-noplatt] [-v] [-verbose] [-h] [-help]
With ``-platt``, Platt's algorithm is used, which may be faster when With ``-platt``, Platt's algorithm is used, which may be faster when
computing many zeros of large index simultaneously. computing many zeros of large index simultaneously.
With ``-stream``, the zeros are computed and output one
by one instead of simultaneously.
bernoulli.c bernoulli.c
------------------------------------------------------------------------------- -------------------------------------------------------------------------------

View file

@ -13,7 +13,7 @@ void print_zeros(arb_srcptr p, const fmpz_t n, slong len, slong digits)
{ {
fmpz_print(k); fmpz_print(k);
flint_printf("\t"); flint_printf("\t");
arb_printn(p+i, digits, 0); arb_printn(p+i, digits, ARB_STR_NO_RADIUS);
flint_printf("\n"); flint_printf("\n");
fmpz_add_ui(k, k, 1); fmpz_add_ui(k, k, 1);
} }
@ -22,7 +22,7 @@ void print_zeros(arb_srcptr p, const fmpz_t n, slong len, slong digits)
void print_help() void print_help()
{ {
flint_printf("zeta_zeros [-n n] [-count n] [-prec n] [-threads n] " flint_printf("zeta_zeros [-n n] [-count n] [-prec n] [-digits n] [-threads n] "
"[-platt] [-noplatt] [-v] [-verbose] [-h] [-help]\n\n"); "[-platt] [-noplatt] [-v] [-verbose] [-h] [-help]\n\n");
flint_printf("Reports the imaginary parts of consecutive nontrivial zeros " flint_printf("Reports the imaginary parts of consecutive nontrivial zeros "
"of the Riemann zeta function.\n"); "of the Riemann zeta function.\n");
@ -49,7 +49,7 @@ int main(int argc, char *argv[])
int verbose = 0; int verbose = 0;
int platt = 0; int platt = 0;
int noplatt = 0; int noplatt = 0;
slong i, buffersize, prec; slong i, buffersize, prec, digits;
fmpz_t requested, count, nstart, n; fmpz_t requested, count, nstart, n;
arb_ptr p; arb_ptr p;
@ -71,6 +71,7 @@ int main(int argc, char *argv[])
fmpz_set_si(requested, -1); fmpz_set_si(requested, -1);
buffersize = max_buffersize; buffersize = max_buffersize;
prec = -1; prec = -1;
digits = 2;
for (i = 1; i < argc; i++) for (i = 1; i < argc; i++)
{ {
@ -122,6 +123,17 @@ int main(int argc, char *argv[])
{ {
requires_value(argc, argv, i); requires_value(argc, argv, i);
prec = atol(argv[i+1]); prec = atol(argv[i+1]);
digits = prec / 3.32192809488736 + 1;
if (prec < 2)
{
invalid_value(argv, i);
}
}
else if (!strcmp(argv[i], "-digits"))
{
requires_value(argc, argv, i);
digits = atol(argv[i+1]);
prec = digits * 3.32192809488736 + 3;
if (prec < 2) if (prec < 2)
{ {
invalid_value(argv, i); invalid_value(argv, i);
@ -173,6 +185,7 @@ int main(int argc, char *argv[])
{ {
prec = 64 + fmpz_clog_ui(nstart, 2); prec = 64 + fmpz_clog_ui(nstart, 2);
if (platt) prec *= 2; if (platt) prec *= 2;
digits = prec / 3.32192809488736 + 1;
} }
if (verbose) if (verbose)
@ -225,7 +238,7 @@ int main(int argc, char *argv[])
num = FLINT_MIN(1 << fmpz_get_si(iter), num); num = FLINT_MIN(1 << fmpz_get_si(iter), num);
} }
acb_dirichlet_hardy_z_zeros(p, n, num, prec); acb_dirichlet_hardy_z_zeros(p, n, num, prec);
print_zeros(p, n, num, prec*0.3 + 2); print_zeros(p, n, num, digits);
fmpz_add_si(n, n, num); fmpz_add_si(n, n, num);
fmpz_add_si(count, count, num); fmpz_add_si(count, count, num);
fmpz_add_si(iter, iter, 1); fmpz_add_si(iter, iter, 1);
@ -259,7 +272,7 @@ int main(int argc, char *argv[])
"is wrong with the internal tuning parameters."); "is wrong with the internal tuning parameters.");
flint_abort(); flint_abort();
} }
print_zeros(p, n, found, prec*0.3 + 2); print_zeros(p, n, found, digits);
fmpz_add_si(n, n, found); fmpz_add_si(n, n, found);
fmpz_add_si(count, count, found); fmpz_add_si(count, count, found);
} }
@ -274,6 +287,6 @@ int main(int argc, char *argv[])
fmpz_clear(requested); fmpz_clear(requested);
fmpz_clear(count); fmpz_clear(count);
flint_cleanup(); flint_cleanup_master();
return 0; return 0;
} }