What is the derivation of the exact Maxwell-Boltzmann distribution?

[rabbed]

I would like to see a derivation of the exact Maxwell-Boltzmann distribution shown as (16) in this document: https://www.researchgate.net/publication/222670999_Exact_Maxwell-Boltzmann_Bose-Einstein_and_Fermi-Dirac_Statistics

This is my starting point (f being the function to maximize, g and h being the constraints, a and b being the lagrange multipliers and nk being the number of particles in level/energy k):

f(n0, ..., n6) = ln(6!/(n0!*n1!*n2!*n3!*n4!*n5!*n6!))
g(n0, ..., n6) = 0*n0 + 1*n1 + 2*n2 + 3*n3 + 4*n4 + 5*n5 + 6*n6 = 6
h(n0, ..., n6) = n0 + n1 + n2 + n3 + n4 + n5 + n6 = 6
gradient(f) = a*gradient(g) + b*gradient(h)

I know ln(gamma(x))' = digamma(x)

 

[rabbed]

I believe the most probable macrostate (for 6 particles, constant energy 6) is achieved with occupancy 3,1,1,1,0,0,0 in energy levels 0,1,2,3,4,5,6.
Is the solution to take the outputs from the following calculation (where I've picked suitable values a=-0.06 and b=1.96 to fit the solution, but maybe there are correct values?) and round them to their nearest integers?
https://www.wolframalpha.com/input/?i=evaluate+-x!'/x!-0.06*x+1.96+at+x=0,+1,+2,+3,+4,+5,+6