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`

).

Are you sure there are any such pairs to find? Unless I'm missing something, you can't fit $2p^2$ into $(a+1) (b+1) (a+b+2)$ successfully.

Ah @DSM this may indeed be the case. Now I suppose I would like to find points where $d(a,b) = n$ for $n \in \mathbb{N}$. Does Sage have a toolkit that allows me to extract integral points on the surface $F(a,b,n) = d(a,b) - n = 0$?