Programming PARI/GP to do a sum

Question

I'm trying to compute the following sum in PARI/GP $C=\sum_{n=1}^{\infty} \frac{g(n)}{n^2}$ where $g(n)$ defined as $$g(n)=(-1)^r, \qquad r=\text{number of even indexed prime factors of $n$}$$ By even indexed primes I mean the primes 3, 7, 13, 19, 29... corresponding to $p_2, p_4, p_6, p_8, p_{10} \ldots$ where $p_n$ is the $n$th prime number. @YiyuanLee suggested the following pari code to compute $g(n)$

compute(n) =
{
    g = 1;
    index = 0;
    p = 0;
    while(p < n, 
        p = nextprime(p + 1); 
        index += 1;
        while(n % p == 0, n = n / p; if (index % 2 == 0, g *= -1));
    );
    printf("%d\n", g);
}

Although it gives individual values of $g(n)$ nicely, I cannot seem to use it to compute $C$. I tried the following in PARI/GP,

sum(n=1, N, compute(n)/n^2)

But all I get is a list of values of $g(n)$ all the way till the value for $N$ and a final output of $0$. So how do I ask PARI to do the sum.

Answer 1

Your procedure 'compute' must first become a function so replace printf("%d\n", g); by g (the function's result or return g)

Your sum is in fact a second function of N so define :

S(N)= sum(n=1, N, compute(n)/n^2)

After that try for example

> S(10)
%1 = 1563689/1270080
> vector(10, n, S(n))
%2 = [1, 5/4, 41/36, 173/144, 4469/3600, 4369/3600, 210481/176400, 852949/705600, 7754941/6350400, 1563689/1270080]

If you want real values multiply compute(n) by 1.0. and try (as indicated by Derek) suminf or rather sumalt (with more precision asked by \p) :

> \p 100
> sumalt(n=1, compute(n)*1.0/n^2)
%3 = 1.230439635196936722876166915075052259785033661462657803428111193673870265619650907934776407938143274

> \p 1000
> sumalt(n=1, compute(n)*1.0/n^2)
%4 = 1.230393982736387671131352032...

> \p 8000
> sumalt(n=1, compute(n)*1.0/n^2)
%5 = 1.23040307381058333073

Increasing the precision will give you more significative digits! The result is near $1.230403$.

Here it seems easier to compute directly some values :

> \p 19
> s(10^4)
%6 = 1.230403171755051041
> s(4*10^4)
%7 = 1.230403247788137930

(pari/gp tutorial should help you to learn more)