Hi just an update on this. The reason that the computation wasn't symmetric was that the previous chosen convention was only asymptotically symmetric. I have now modified this such that the case for AlternativeHypothesis->"Unequal" is symmetric for any Method. I have also added an "Exact" method following original sources (even though they are missing boundary conditions for the derivation), which matches the tables in the original sources. I have chosen a convention to be as consistent as possible with the preexisting "Asymptotic" computation. Unfortunately it won't make it into 14.0.

Thank you very much. One more request: as long as you will fix the MannWhitneyTest function, it would be great also to fix the behavior of Method->"Automatic". Now it always calculates p-value from asymptotic normal distribution for U. At the same time for small samples (<20) an exact p-value could be calculated based on permutations (not Monte Carlo, a direct formula/table instead). For instance, in the reported case of data1, data2 the exact p-value is the probability of getting 2 extreme cases: all ranks of sample 1 below the sample 2 and all ranks of sample 1 above the sample 2, which is exactly 2(6!6!)/(6+6)! = 1/462 ~ 0.00216.The asymptotic normal approximation for U gives 0.003053. This difference is not a big deal if a single M-W test is performed as both results are significant. But if multiple comparisons are required, then multiple M-W tests are applied with subsequent correction of p-values by for instance Holm–Bonferroni method. And here the difference between 0.002 and 0.003 might become important.

That sucks!

https://support.wolfram.com/12507?src=mathematica

Is there any way to report this bug? It is extremely influential. I can't imagine how many datasets from animal studies I've processed with M-W in Mathematica in last 5 years... It will be a nightmare to revisit all the data once again, taking into account that a lot of data is already published. I'm shocked.