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?