Python Scipy multinomial.pmf function different results?

Question

I was working on figuring out an exact multinomial test of goodness-of-fit with Python (3.9) in JupyterLab via Anaconda on Windows 10. To do this, I need the probability mass function of the multinomial distribution. Scipy has a nice function for this. I've set the following:

from scipy.stats import multinomial
multinomial.pmf(x=[16, 17, 22], n=55, p=[1/3, 1/3, 1/3])

The result of this is 0.008700756210332266

Now, note that the expected proportions are equal for each category, so the order of the observed counts should not matter, however:

from scipy.stats import multinomial
multinomial.pmf(x=[22, 17, 16], n=55, p=[1/3, 1/3, 1/3])

This results in: 0.00870075621033202

A very small difference. However, for the exact multinomial test, I'm iterating over a lot of combinations and checking if the value is less than or equal to this. It turns out these quickly add up.

I've solved the issue myself by sorting each option, and another option would of course be to simply round the decimals to 8 or so.

The question: Why does the multinomial.pmf function from Scipy give different results when they should be the same?

Thanks in advance for any comments/suggestions.