The answer to your very last question is yes. As for the rest:
Computing the parity up to N is, in theory, quasi linear by using Pentagonal Number Theorem: compute the inverse of the pentagonal number power series mod $x^N$. This is of course faster than the quadratic-complexity "en masse" method of using simply the recurrence relation. I am not sure how inverse_mod is actually computed in sage, but it ran fine for the following code:
def get_parity(n):
x = GF(2)['x'].0
f = 1
for k in xrange(1,sqrt(2*n/3)+10):
f += x^((k*(3*k-1))/2) + x^((-k*(-3*k-1))/2)
return f.inverse_mod(x^n).coeffs()
def get_best(l):
n = len(l)
highest, lowest = [0,0,0], [0,0,1]
for a in [1..sqrt(n)]:
for b in [0..a-1]:
score = (l[b::a].count(1) + 0.0)*a/(n-b)
if score > highest[2]:
highest = [a, b, score]
elif score < lowest[2]:
lowest = [a, b, score]
return highest, lowest
l = get_parity(2^22)
best = get_best(l[:2^21])
for x in best:
print x[2], (l[x[1]::x[0]].count(1)+0.0)*x[0]/(len(l)-x[1])
What this does is calculate for every a.p. $an+b$ the density of odd values, and finds the two a.p.'s with most and least odd values.
get_parity took more than 20 minutes (not sure exactly since I was watching TV), and get_best took almost an hour (again, I think). I'm running a macbook pro with 4gb ram.
The results were:
(1442, 766) 0.557846694263366 0.531611381990907
(1389, 357) 0.440522320970815 0.468967411597574
This is to say that the most special a.p.'s up to $2^{21}$ become a bit less special when going up to $2^{22}$, and quite close to equidistribution. Hence, I would say that yes, the computer does suggest that the analogous result holds when we restrict n to lie in a fixed arithmetic progression.
Edit 1: If we change the "score" of an a.p. to:
$$\frac{\\\#\{odd\\ in\\ a.p.\} - \\\#\{even\\ in\\ a.p.\}}{\sqrt{length}}$$
Then up to $2^{21}$ the largest values are:
(712, 254) 4.84633129231862
(1389, 357) -4.60795692951670