Maple Procedure to plot graph of prime number twin primes smaller than x at x

Question

I'm in my first year studying for my BSc in maths and I have the following problem to solve by Sunday: Write a Maple procedure with the entry of a real number x that will count the number of twin primes (p,p+2) with p+2<=x. Plot the results in a procedure between 1 and 1 000 000 I have the procedure to find the number of twin primes, but I can't seem to put it into a graph. Here's my progress:

primzahl:=proc(n)

local a,Q;

Q:=0;

for a from 2 to n-2 do

if isprime(a) then

if isprime(a+2) then

Q:=Q+1;

fi;

fi;

od;

end:

Thanks in advance,

emj

Answer 1

Your routine calls isprime on every a between 2 and n-2. It could be done faster by using nextprime to get to the next candidate a to test.

Also, it's much better to have the primzahl procedure keep a running total, so that it can store the values where new twin primes are found. It looked like you intended on calling your original procedure many times, to figure out the totals separately. That'd be unnecessarily costly.

restart;
primzahl:=proc(n)
  local a,prev,Q;
  Q[2]:=0;
  prev:=2;
  for a from 2 to n-2 do
    if isprime(a) and isprime(a+2) then
      Q[a+2]:=Q[prev]+1;
      prev:=a+2;
    fi;
  od;
  Q[2]:='Q[2]';
  eval(Q);
end:

R:=primzahl(100);

st:=time():
  R:=primzahl(1000000): # suppress this large output using colon
time()-st;

plot([seq([j,R[j]], j=[indices(R,nolist)])],
     style=point, symbolsize=1);

restart;

primzahl2:=proc(n)
  local a,Q,total;
  a:=2;
  total:=0;
  while a <= n-2 do
    if a+2<n-2 and isprime(a+2) then
      total:=total+1;
      Q[a+2]:=total;
    fi;
    a:=nextprime(a);
  od;
  eval(Q);
end:

R2:=primzahl2(100);

st:=time():
  R2:=primzahl2(1000000): # suppress this large output using colon
time()-st;

plot([seq([j,R2[j]], j=[indices(R2,nolist)])],
     style=point, symbolsize=1);

[edited] And another minor adjustment skips the nextprime call in the case that a twin prime pair was found (as the next candidate a can just be taken as a+2).

restart;

primzahl3:=proc(n)
  local a,Q,total;
  a:=2;
  total:=0;
  while a <= n-2 do
    if a+2<n-2 and isprime(a+2) then
      total:=total+1;
      Q[a+2]:=total;
      a:=a+2;
    else;
      a:=nextprime(a);
    end if;
  od;
  eval(Q);
end:

R3:=primzahl3(100);

st:=time():
  R3:=primzahl3(1000000): # suppress this large output
time()-st;

plot([seq([j,R3[j]], j=[indices(R3,nolist)])],
     style=point,symbolsize=1);