How to find instances where $d(a,b) = p^2$ for $p$ a prime

edit

I am not sure i understand your question correctly. If you want to iterate over the pairs of integers, you can have a look at multidimensional enumeration:

sage: f = lambda a,b : (a+1)*(b+1)*(a+b+2)
sage: n = 720

sage: for a,b in xmrange([n,n]):
sage:     if f(a,b) == 2*n:
sage:         print a,b      
7 9
9 7

And actually you can be more restrictive and ask a and b to be less than floor(sqrt(2*n)) (instead of n):

sage: bound = floor(sqrt(2*n))
sage: for a,b in xmrange([bound,bound]):
sage:     if f(a,b) == 2*n:
sage:         print a,b      
7 9
9 7

Hence, you get all answers in time $O(n)$ instead of $O(n^2)$. But why not solving directly your equation ? Unfortunately, it seems that the solver is not able to solve the diophantine equation f(a,b) == 2*n directly:

sage: var('a','b')
(a, b)
sage: solve([f(a,b) == 2*n ], a, b)
([a == -1/2*(4*b + sqrt(b^4 + 4*b^3 + 6*b^2 + 5764*b + 5761) + b^2 + 3)/(b + 1), a == -1/2*(4*b - sqrt(b^4 + 4*b^3 + 6*b^2 + 5764*b + 5761) + b^2 + 3)/(b + 1)],
 [1, 1])

But you can still use Sage solver when a is fixed:

sage: var('b')
b
sage: assume(b, 'integer')
sage: assume(b >= 0)

sage: for a in xrange(bound):
sage:     sol = solve([f(a,b) == 2*n],b)
sage:     if len(sol) != 0:
sage:         print a, sol[0].rhs()
sage:         
7 9
9 7

Now, you get an answer in time $O(\sqrt{n})$, which is not completely satisfactory (one could expect a polynomial in $\log(n)$), but not that bad.

That said, since the solver is quite slow, the multiplicative constant in the complexity is quite huge and the previous method is faster for small n (this is clear for n=6216 : 42.7ms vs 1.11s). The last method becomes faster only later (this is clear for n=3065451 : 94.43s vs 29.77s).