The +999 will make the result go up by 1 except when xy is exactly 0 mod 1000.
It's easy to calculate xy mod 1000 without overflow, as (x mod 1000)*(y mod 1000). So we have two cases:
xy = 0 mod 1000: then your f(x,y) is the same as xy div 1000.
xy != 0 mod 1000: then your f(x,y) is the same as (xy div 1000) + 1
So we can focus on just calculating xy div ...