Evaluate the sum

Question

I want to evaluate the sum $$\sum_{k=1}^\infty \left(\frac{i+1}{\sqrt{2}}\right)^k\cdot k^{-\alpha}$$ where $i=\sqrt{-1}$ and $\alpha\in[\frac{3}{4},1]$ with 8 digits accuracy.

If I am willing to expend up to a second of CPU time per calculation, how should I proceed?

Suppose I want a really fast scheme for the evaluation and I am willing to expend a lot of work. How should I proceed?

Any leads, clues or suggestions would be greatly appreciated.

Answer 1

What have you tried so far? This completely naive implementation manages to compute 7 (maybe 7.5) digits in 2.5 seconds on my laptop:

#include <iostream>
#include <complex>
#include <cmath>
#include <iomanip>

int main ()
{
  const double alpha = 1;
  std::cout.precision(16);

  std::complex<double> sum = 0;
  for (unsigned int k=1; k<10000000; ++k)
    {
      sum += std::pow(std::complex<double>(1,1)/std::sqrt(2.), k)
         *
         std::pow(k, -alpha);

      if (k % 1000000 == 0)
    std::cout << k << ' ' << sum << std::endl;
    }
}

Results are:

g/x> c++ -O3 x.cc
g/x> time ./a.out 
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)

real    0m2.553s
user    0m2.548s
sys     0m0.001s

I'm certain this can be vastly improved, for example by using the fact that $(1+i)^k$ can only attain 8 different values that could be tabulated, that $\sqrt{2}^k$ can be computed very easily, etc. I have no doubt that with 20 minutes of work I could make this run in under one second for 8 digits. It might be harder to do for smaller values of $\alpha$, however, where convergence is slower.

Answer 2

Since the question has a Fortran tag, here is Wolfgang's solution in Fortran:

implicit none
integer, parameter :: dp = kind(0.d0)
complex(dp), parameter :: i_ = (0, 1)

real(dp) :: alpha = 1
complex(dp) :: s = 0
integer :: k
do k = 1, 10000000
    s = s + ((i_+1)/sqrt(2._dp))**k * k**(-alpha)
    if (modulo(k, 1000000) == 0) print *, k, s
end do
end

Result are:

$ gfortran -O3 x.f90 
$ time ./a.out 
     1000000 ( 0.26740049836809593     ,  1.1780960379893362     )
     2000000 ( 0.26740024836936921     ,  1.1780966415427796     )
     3000000 ( 0.26740016503629038     ,  1.1780968427272547     )
     4000000 ( 0.26740012336968905     ,  1.1780969433194985     )
     5000000 ( 0.26740009836973161     ,  1.1780970036749265     )
     6000000 ( 0.26740008170307927     ,  1.1780970439117682     )
     7000000 ( 0.26740006979833991     ,  1.1780970726523770     )
     8000000 ( 0.26740006086978074     ,  1.1780970942078899     )
     9000000 ( 0.26740005392532235     ,  1.1780971109732488     )
    10000000 ( 0.26740004836978204     ,  1.1780971243856078     )

real    0m1.477s
user    0m1.472s
sys 0m0.000s

For comparison, on my computer the C++ code takes:

$ g++ -O3 x.cc
$ time ./a.out 
1000000 (0.2674004983680959,1.178096037989336)
2000000 (0.2674002483693692,1.17809664154278)
3000000 (0.2674001650362904,1.178096842727255)
4000000 (0.267400123369689,1.178096943319499)
5000000 (0.2674000983697316,1.178097003674927)
6000000 (0.2674000817030793,1.178097043911768)
7000000 (0.2674000697983399,1.178097072652377)
8000000 (0.2674000608697807,1.17809709420789)
9000000 (0.2674000539253223,1.178097110973249)

real    0m2.123s
user    0m2.116s
sys 0m0.004s

I use gcc 4.6.3.

Answer 3

mpmath has this

>>> import mpmath
>>> mpmath.polylog(1, (1+1j)/mpmath.sqrt(2))
mpc(real='0.26739999836978523', imag='1.1780972450961724')
>>> import numpy
>>> for x in numpy.linspace(0.75, 1.0, 11): print mpmath.polylog(x, mpmath.sqrt(1j))
... 
(0.13205298513153 + 1.22168160828646j)
(0.147070653672424 + 1.21783198879211j)
(0.161744021745461 + 1.21384547018142j)
(0.17607856182228 + 1.20973163332567j)
(0.190079791459314 + 1.20549967730407j)
(0.2037532636373 + 1.20115843051452j)
(0.217104557681381 + 1.19671636161604j)
(0.230139270736483 + 1.19218159029609j)
(0.242863009773464 + 1.18756189785626j)
(0.25528138410235 + 1.18286473761121j)
(0.267399998369785 + 1.17809724509617j)
>>> 

Answer 4

If you plot the partial sums for the series, you get the following result:

This is the real part for $\alpha=0.91$ for partial sums of 10 to 100 terms. Clearly, there is a pattern here, and the imaginary part follows a similar pattern. It seems that you can probably take 8 successive partial sums and average them to get a much better approximation. Or at least you can get the asymptotic numeric behavior from these partial sums.

More analytically, you can break the sum into 8 separate sums over the different roots of unity:

$$ \sum_{j=1}^8 \left( \frac{i+1}{\sqrt{2}} \right)^j \sum_{k=0}^\infty (8k+j)^{-\alpha} $$

Each of the inner sums is a Hurwitz zeta function, for which there are special evaluation routines (like in GSL). Plotting out the sum in this way as a function of $\alpha$ gives

The blue line is the real part, red is the imaginary part. This seems pretty well behaved.

Answer 5

The series is divergent for $\alpha < 1$. Its convergence is suspicious using the ratio test. For $\alpha = 1$, the sum is $-\log(1 - \sqrt{i})$, which is 0.2673999983697852 + 1.1780972450961725$i$, close to what is given in other answers. Results from Mathematica. By the way, $\frac{1+i}{\sqrt{2}}$ is $\sqrt{i}$, which also explains that its powers attain only 8 distinct values.