(Maxima) load(basic)$ a:[]$ p:-1$ j:0$ m:0$
chli():= block (if w>341550071728321 then
(n:11000, print("# over limit") ), return)$
for n:1 thru 10000 step 0 do
(p:p+6, q:p+1, r:p+2, if (primep(p) and primep(r)) then
(push(q, a), n:n+1, j:j+1) else
(w:3*q, chli(), if (primep(w-1) and primep(w+1)) then
(push(q, a), n:n+1, m:m+1
) ) )$
a:reverse(a)$ d:length(a)$ k:float(m*100/d)$ h:", "$
y:last(a)$ b:float(d*100/(y/6))$
print("# Real ATPP = ", j, h, " Pseudo ATPP = ", m, h, " Percent pseudo = ", k)$
print("# First ", d, " of sequence are ", b, "% of ints. up to ", y, " cong. to 0 mod 6.")$ for i:1 thru d do (s:pop(a), print(i, h, s) )$
