Hi! Let f[i={1,...,N}][j={1,...,M}] be the number of ways if you have the first i children and the i-th has j candies. That is M*N=10^12 states in total. My approach uses a state f[i={1,...,N}][j={1,...sqrt(M)}][k={0,1}] that is the answer if we have the first i children and if k=0 j is the number of the candies that were given to the i-th child. If k=1 j is the maximum number of candies that you can give to the (i-1)-th child, so if you are in state (i,j,1) you know that you will give M/j to the i-th child.

f[1][1,...,sqrt(M)][0,1]=1

f[i][j][0]=number of ways (i-1)th child can take 1,...,M/j candies.

f[i][j][1]=number of ways (i-1)th child can take 1,...,j candies.

So:

f[i][j][0]=f[i-1][j][0]+...+f[i-1][min(sqrt(M),M/j)][0]+f[i-1][1][1]+...+f[i-1][j][1]

and

f[i][j][1]=f[i-1][j][0]+...+f[i-1][min(sqrt(M),j)][0]+f[i-1][1][1]+...+f[i-1][M/j][1]

You can compute each state in O(1) if you have the sums from the last step. I hope I didn't make a mistake, but if I did you have the idea!