How is the output size of UPFIRDN determined?

I want to understand resampling better, so I am trying to replicate the functionality of Matlab's built-in upfirdn() function. To clarify, I'm not looking for a fast implementation (at first), I just want to see how to get the exact same result. As such, I understand the basic steps (conceptually) are:

          (1) Upsampling by P (zero insertion).
          (2) FIR filtering.
          (3) Downsampling by Q (throwing away samples).

This is fine. However, what is bugging me is that (as stated by Matlab help) the final output size is:

           Ly = ceil( ((Lx-1)*P + Lh)/Q ) where Lx = length(X) and Lh = length(H)

and I don't understand why. With reference to the 3 stages listed above, a naive implementation might give stage-by-stage output sizes:

          (1) P*Lx
          (2) P*Lx
          (3) floor(P*Lx/Q)

Therefore, I really want to know exactly what upfirdn is doing. Can anyone help?

[Bradford Watson]

Whenever two sequences are convolved together, the length of the resulting sequence is Lx + Ly - 1. The inner part of this equation is expressing the length of the convolution of the upsampled sequence (Lx-1)*P and the filter length (Lh), though neglecting the -1 term. This is then divided by the decimation, and the ceiling function is taken to get rid of any fraction.