This problem's solution ends up at

nCr(A+B-i, i) * nCr(A+B-2i, A-i) * pow(x, A-i) * pow(y, B-i) * pow(z, i) ; for every 0 <= i <= min(A, B)

All inputs are between 1 and 10^17 inclusive, and result should be printed as (result % M); where M = 1e6 + 3

I've seen a solution like this

while(A || B)
{
        u = A % M;
        v = B % M;
        tmp = 0;
        for ( i = 0; i <= min(u, v); i++)
        {
                tmp = (tmp + nCr(u+v-i, i) * nCr(u+v-2i, u-i) * pow(x, u-i) * pow(y, v-i) * pow(z, i)) % M;
        }
        ans = (ans * tmp) % M;
        A /= M;
        B /= M;
}

Full Source Code: link

Can you please explain how did he derive this?