184 Reputation

2 Badges

16 years, 113 days

MaplePrimes Activity

These are answers submitted by tom_

Good to know, thanks again alec.

Ah thanks alec, much appreciated.  It worked out nicely in Maple also.  I'll spare you the plots.

Nevermind, seems as though Maple is happy with just defining the interval twice as long as the array.  But I also created a function 'fn:=x -> array[x/2]' and passed that as the function to plot.  I didn't think this would work but it has.

OK, with a lot of effort I think I have it.  I got the Hans Riesel book "prime numbers and computer methods for factorization" out of my university library and I found a function he uses to correct the error -

where rho_k = 0.5 + i*alpha_k is the k-th non-trivial zero.

I just tried this in C using the GNU GSL library and it seems to do the trick.  Basically, I just did R_n(x) = R(x) + [sum, k=1..n] C_k(x) and plotted the values I got  (I'm not sure if I was supposed to use the trivial zeroes as well. I probably was. Ah well). It was nice and quick too.  If you plot it against pi(x) it gets better the higher n is. 

Here is my program - .  Ugly code, but it does the trick. It calculates Rn(x) for 0<x<100 and just spits out the values.  I dumped these values alongside corresponding values for pi(x) in a datafile for gnuplot to do its work. Here are some graphs for  n=10 and n=100 respectively:

Cool huh? I'll try it in Maple on Thursday.

May I ask how you arrived at that -0.5 + arctan()/pi formula? I can't find this anywhere, are you sure it is correct.  My C program is not going well either.  I have written what I think should do the trick but I am getting strange results. *sigh*

Your suggestion to use add() instead of sum() was a good one as R_n(x) now is real, whereas I was getting some non-zero imaginary part before.  However, it is just far too slow and inprecise that it seems too fruitless to continue.  I am getting negative values for R_n(x) too, so something is going wrong.  For example, pi(50)=15, while R_n(50,10) = -2.238047486 10^22.

I like the idea of writing it in C. At home I run linux, no windows at all, so if I was just going to calculate values using a C program, then I could just give it to gnuplot, and discard Maple.  I downloaded the GNU scientific library last night - and so I will try and see how that goes.

Sorry it just dawned on me you may have meant so that you could test it yourself. So here is 'Zeroes' if you care


Zeroes := [14.134725,21.022039,25.010857,30.424876,

Originally they were to 1000 decimal places, so I wrote a C program to chop them down to 6 and place them in a list like the above for me.

Thanks for the replies, appreciated.

> It would be easier to help if you directly provided what you call "Zeroes".

'Zeroes' is a list of the imaginary part of the first 100 non-trivial zeroes.  So literally, a list with 100 elements in it.  I got these from some website.

> the sum to infinity has to be used. It is very slow in Maple though.

Precisely why I didn't use it.  I only want to produce a plot for small x anyway.  The point to all this is I am writing a paper on prime numbers and I need to show that I have produced my own plots.

I assumed I would have to use complexplot since we have the non-trivial zeroes in there which are complex?

> Rx := x -> evalf(add( ( ln(x) )^k / ( k * k! * evalf( Zeta( k+1 ) ) ) , k = 1..100 ));

I will try this tomorrow as I don't have access to Maple at the moment.  Thanks.

> The sum over trivial zeros is......

Are you sure?  Is it worth replacing this with sum( R( x^(-2*m) ), m=1..whatever )?

Anyway, thanks for your reponses, much appreciated.

Create your own thread.

Thanks for that.  The question still stands though.

That's great but I do not have access to Mathematica.

When I say non-trivial zeroes, I mean zeroes of the Riemann zeta function which have real part 1/2.

Aha, I have it.



I have tried the following but it gives me an empty graph


plot([Zeta(0.5 + I*t),t,t=0..50],coords=polar);

I found a solution:


plot( pi(floor(x)), x=1..100, numpoints = 100, axesfont=[HELVETICA, SYMBOL,TIMES], thickness=2, axes=boxed );

1 2 Page 1 of 2