Skip to content

Commit

Permalink
use sieved power sum in zeta for short power series; some refactoring
Browse files Browse the repository at this point in the history
  • Loading branch information
@fredrik-johansson
fredrik-johansson committed Aug 2, 2013
1 parent 0659224 commit 6e8e2d2
Show file tree
Hide file tree
Showing 5 changed files with 324 additions and 36 deletions.
3 changes: 3 additions & 0 deletions zeta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,5 +50,8 @@ void zeta_series_em_bound(fmpr_t bound, const fmpcb_t s, const fmpcb_t a, long N
void zeta_series_em_vec_bound(fmprb_ptr vec, const fmpcb_t s, const fmpcb_t a, ulong N, ulong M, long d, long wp);
void zeta_series(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec);

void zeta_powsum_series_naive(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, long n, long len, long prec);
void zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, long prec);

#endif

164 changes: 164 additions & 0 deletions zeta/powsum_one_series_sieved.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,164 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/

#include "zeta.h"
#include "fmpcb.h"
#include "fmpcb_poly.h"

void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_ptr res,
const fmpcb_t N, const fmpcb_t c, long trunc, long prec);


#define POWER(_k) (powers + (((_k)-1)/2) * (len))
#define DIVISOR(_k) (divisors[((_k)-1)/2])

#define COMPUTE_POWER(t, k) \
do { \
if (integer) \
{ \
fmprb_neg(w, fmpcb_realref(s)); \
fmprb_set_ui(v, k); \
fmprb_pow(fmpcb_realref(t), v, w, prec); \
fmprb_zero(fmpcb_imagref(t)); \
if (len != 1) \
{ \
fmprb_log_ui(logk, k, prec); \
fmprb_neg(logk, logk); \
} \
} \
else \
{ \
fmprb_log_ui(logk, k, prec); \
fmprb_neg(logk, logk); \
fmprb_mul(w, logk, fmpcb_imagref(s), prec); \
fmprb_sin_cos(fmpcb_imagref(t), fmpcb_realref(t), w, prec); \
if (critical_line) \
{ \
fmprb_rsqrt_ui(w, k, prec); \
fmpcb_mul_fmprb(t, t, w, prec); \
} \
else \
{ \
fmprb_mul(w, fmpcb_realref(s), logk, prec); \
fmprb_exp(w, w, prec); \
fmpcb_mul_fmprb(t, t, w, prec); \
} \
} \
for (i = 1; i < len; i++) \
{ \
fmpcb_mul_fmprb(t + i, t + i - 1, logk, prec); \
fmpcb_div_ui(t + i, t + i, i, prec); \
} \
} while (0); \

void
zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, long prec)
{
long * divisors;
long powers_alloc;
long i, j, k, power_of_two, horner_point;
int critical_line, integer;

fmpcb_ptr powers;
fmpcb_ptr t, u, x;
fmpcb_ptr p1, p2;
fmprb_t logk, v, w;

critical_line = fmprb_is_exact(fmpcb_realref(s)) &&
(fmpr_cmp_2exp_si(fmprb_midref(fmpcb_realref(s)), -1) == 0);

integer = fmprb_is_zero(fmpcb_imagref(s)) && fmprb_is_int(fmpcb_realref(s));

divisors = flint_calloc(n / 2 + 1, sizeof(long));
powers_alloc = (n / 4 + 1) * len;
powers = _fmpcb_vec_init(powers_alloc);

for (i = 3; i <= n; i += 2)
for (j = 3 * i; j <= n; j += 2 * i)
DIVISOR(j) = i;

t = _fmpcb_vec_init(len);
u = _fmpcb_vec_init(len);
x = _fmpcb_vec_init(len);
fmprb_init(logk);
fmprb_init(v);
fmprb_init(w);

power_of_two = 1;
while (power_of_two * 2 <= n)
power_of_two *= 2;
horner_point = n / power_of_two;

_fmpcb_vec_zero(z, len);

COMPUTE_POWER(x, 2);

for (k = 1; k <= n; k += 2)
{
/* t = k^(-s) */
if (DIVISOR(k) == 0)
{
COMPUTE_POWER(t, k);
}
else
{
p1 = POWER(DIVISOR(k));
p2 = POWER(k / DIVISOR(k));

if (len == 1)
fmpcb_mul(t, p1, p2, prec);
else
_fmpcb_poly_mullow(t, p1, len, p2, len, len, prec);
}

if (k * 2 < n)
_fmpcb_vec_set(POWER(k), t, len);

_fmpcb_vec_add(u, u, t, len, prec);

while (k == horner_point && power_of_two != 1)
{
_fmpcb_poly_mullow(t, z, len, x, len, len, prec);
_fmpcb_vec_add(z, t, u, len, prec);

power_of_two /= 2;
horner_point = n / power_of_two;
horner_point -= (horner_point % 2 == 0);
}
}

_fmpcb_poly_mullow(t, z, len, x, len, len, prec);
_fmpcb_vec_add(z, t, u, len, prec);

flint_free(divisors);
_fmpcb_vec_clear(powers, powers_alloc);
_fmpcb_vec_clear(t, len);
_fmpcb_vec_clear(u, len);
_fmpcb_vec_clear(x, len);
fmprb_clear(logk);
fmprb_clear(v);
fmprb_clear(w);
}

56 changes: 56 additions & 0 deletions zeta/powsum_series_naive.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/

#include "zeta.h"
#include "fmpcb.h"
#include "fmpcb_poly.h"

void _fmpcb_poly_fmpcb_invpow_cpx(fmpcb_ptr res,
const fmpcb_t N, const fmpcb_t c, long trunc, long prec);

void
zeta_powsum_series_naive(fmpcb_ptr z,
const fmpcb_t s, const fmpcb_t a, long n, long len, long prec)
{
long k;
fmpcb_ptr t;
fmpcb_t c;

t = _fmpcb_vec_init(len);
fmpcb_init(c);

_fmpcb_vec_zero(z, len);

for (k = 0; k < n; k++)
{
fmpcb_add_ui(c, a, k, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, c, s, len, prec);
_fmpcb_vec_add(z, z, t, len, prec);
}

_fmpcb_vec_clear(t, len);
fmpcb_clear(c);
}

42 changes: 6 additions & 36 deletions zeta/series_em_sum.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,7 @@ zeta_series_em_sum(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, u
fmprb_t x;
fmpz_t c;
long i;
ulong r, n;
ulong r;
int aint;

BERNOULLI_ENSURE_CACHED(2 * M);
Expand All @@ -129,47 +129,17 @@ zeta_series_em_sum(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, u
fmprb_init(x);
fmpz_init(c);

/* sum 1/(n+a)^(s+x) */
if (fmpcb_is_one(a) && (d == 1))
/* sum 1/(k+a)^(s+x) */
if (fmpcb_is_one(a) && d <= 3)
{
fmpcb_ptr pows;
long j;

pows = _fmpcb_vec_init(N + 1);
fmpcb_one(pows + 1);

for (i = 2; i <= N; i++)
{
if (fmpcb_is_zero(pows + i))
{
fmprb_log_ui(fmpcb_realref(pows + i), i, prec);
fmprb_zero(fmpcb_imagref(pows + i));
fmpcb_mul(pows + i, pows + i, s, prec);
fmpcb_neg(pows + i, pows + i);
fmpcb_exp(pows + i, pows + i, prec);
}

for (j = 2; j <= i && i * j <= N; j++)
if (fmpcb_is_zero(pows + i * j))
fmpcb_mul(pows + i * j, pows + i, pows + j, prec);
}

for (i = 1; i <= N; i++)
fmpcb_add(sum, sum, pows + i, prec);

_fmpcb_vec_clear(pows, N + 1);
zeta_powsum_one_series_sieved(sum, s, N, d, prec);
}
else
{
for (n = 0; n < N; n++)
{
/* printf("sum 1: %ld %ld\n", n, N); */
fmpcb_add_ui(Na, a, n, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, Na, s, d, prec);
_fmpcb_vec_add(sum, sum, t, d, prec);
}
zeta_powsum_series_naive(sum, s, a, N, d, prec);
}


/* t = 1/(N+a)^(s+x); we might need one extra term for deflation */
fmpcb_add_ui(Na, a, N, prec);
_fmpcb_poly_fmpcb_invpow_cpx(t, Na, s, d + 1, prec);
Expand Down
95 changes: 95 additions & 0 deletions zeta/test/t-powsum_one_series_sieved.c
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB 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 General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2013 Fredrik Johansson
******************************************************************************/

#include "zeta.h"

int main()
{
long iter;
flint_rand_t state;

printf("powsum_one_series_sieved....");
fflush(stdout);

flint_randinit(state);

for (iter = 0; iter < 2000; iter++)
{
fmpcb_t s, a;
fmpcb_ptr z1, z2;
long i, n, len, prec;

fmpcb_init(s);
fmpcb_init(a);

if (n_randint(state, 2))
{
fmpcb_randtest(s, state, 1 + n_randint(state, 200), 3);
}
else
{
fmprb_set_ui(fmpcb_realref(s), 1);
fmprb_mul_2exp_si(fmpcb_realref(s), fmpcb_realref(s), -1);
fmprb_randtest(fmpcb_imagref(s), state, 1 + n_randint(state, 200), 4);
}

fmpcb_one(a);

prec = 2 + n_randint(state, 200);
n = n_randint(state, 100);
len = 1 + n_randint(state, 10);

z1 = _fmpcb_vec_init(len);
z2 = _fmpcb_vec_init(len);

zeta_powsum_series_naive(z1, s, a, n, len, prec);
zeta_powsum_one_series_sieved(z2, s, n, len, prec);

for (i = 0; i < len; i++)
{
if (!fmpcb_overlaps(z1 + i, z2 + i))
{
printf("FAIL: overlap\n\n");
printf("iter = %ld\n", iter);
printf("n = %ld, prec = %ld, len = %ld, i = %ld\n\n", n, prec, len, i);
printf("s = "); fmpcb_printd(s, prec / 3.33); printf("\n\n");
printf("a = "); fmpcb_printd(a, prec / 3.33); printf("\n\n");
printf("z1 = "); fmpcb_printd(z1 + i, prec / 3.33); printf("\n\n");
printf("z2 = "); fmpcb_printd(z2 + i, prec / 3.33); printf("\n\n");
abort();
}
}

fmpcb_clear(a);
fmpcb_clear(s);
_fmpcb_vec_clear(z1, len);
_fmpcb_vec_clear(z2, len);
}

flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

0 comments on commit 6e8e2d2

Please sign in to comment.