# Evaluate the sum

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.

• Welcome, Qomo! Can you give us some background on what you have already done? – Deathbreath Mar 8 '13 at 18:59
• Also some context might be helpful. – Deathbreath Mar 8 '13 at 19:27
• @Deathbreath. I actually have no lead so far. Since this is a infinite series, then the series must converge at some point up to the digits of accuracy I want. I used the wolfram alpha, and I found the sum is actually called Polylogarithm. From Wikipedia, I found that I can change the sum to the integral of the Bose–Einstein distribution. However, I don't find anything that could help me solve it. – Qomo Mar 9 '13 at 5:11
• "Since this is a infinite series, then the series must converge at some point [...]" Uhm...that does not follow without some more conditions. Just looking at your series is I think it might, but I haven't proven it. – dmckee --- ex-moderator kitten Mar 9 '13 at 16:13
• What you have is more commonly known as a periodic zeta function, which, as you note, is a special case of the polylogarithm. If push comes to shove, you might be able to use a sequence transformation method on partial sums of the periodic zeta function for computational purposes. – J. M. May 14 '13 at 18:01

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.

• You clearly get this, but ... "using the fact that $(1+i)^k$ can only attain 8": The base there has non-unit length, it is the whole expression raised to the power $k$ that only takes on 8 values. – dmckee --- ex-moderator kitten Mar 9 '13 at 16:27
• Sorry, it's actually $\left(frac{1+i}{\sqrt{2}}\right)^k$ that can only attain 8 different values. – Wolfgang Bangerth Mar 10 '13 at 21:34

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.

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)
>>>


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.

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.

• The series is not divergent, it is the definition of the polylogarithm. – Kirill May 15 '13 at 20:38
• You're correct. I believed Mathematica when it said "Sum::div: Sum does not converge." for $\alpha < 1$. – user335 May 18 '13 at 16:41