How do I get an accurate probability when my program that tests in the millions returns significantly different results each time? Apologies in advance, I have never (seriously) studied statistics and I'm out of my depth, so this may not be a simple question.
I wrote a program in java to determine the likelihood of a Tenhou in Mahjong - it's basically just a very rare hand.  I set it to run a series of tests simulating a hand each time, a million at a time.  Basically, it deals out tiles - if it's a Tenhou, that's a pass, otherwise it fails.  Sometimes, there is only one Tenhou out of a million.  Other times, there are 10 or more out of a million.  If I run it for say, 50 million, the results still vary quite a bit.
I'd like to get an accurate probability.  Is there some kind of formula or principle that I can apply?  Even if I could just get a link to an article to read, that would be helpful - I just don't even know where to start.
Edit: Imagine you are playing 5 card poker.  You are dealt five cards.  If it's not a royal flush, shuffle the cards and try again.  That's essentially all my program is doing, but it's doing it one million times, several times.
To be more specific - it shuffles the 136 tiles in a Mahjong set, deals out 14 of them.  If those 14 don't constitute a winning hand, shuffle them all again and deal out a new hand of 14.  Do that a million times, several times.  So, it's not simulating a real game with many different variables - there's just the one variable (the hand that is dealt) and one parameter (condition for a winning hand).
My issue is that when running that same process over and over, the outcome sometimes varies by a factor of 10 or more - that is, there might be one winning hand out of one million dealt, or there may be ten or more winning hands out of one million dealt.  My question is, how do I reconcile a spread that large? Please let me know if more clarification is needed.
 A: Suppose we want to know the probability of getting 4 Aces in a 5-card poker hand
dealt from a standard American deck. By easy combinatorics that's
$48/{52\choose 5} = 48/2\,598\,960 = 1.846893e-05$ or a little less than $2$ chances in $100\,000.$
dhyper(4, 4,48, 5)
[1] 1.846893e-05

It is easy to write a program in R to simulate a million hands, but
difficult to get an approximation by simulation that has small relative error.
set.seed(2021)
deck = c(1,1,1,1, rep(0,48))
nr.a = replicate(10^6, sum(sample(deck,5)))
mean(nr.a==4)
[1] 2e-05        # aprx 0.000018
2*sd(nr.a==4)/1000
[1] 8.944187e-06 # aprx 95% margin of sim error

So we get $0.00002 \pm 0.00001.$ This interval contains the correct value,
but the margin of simulation error is about half as large as the quantity being simulated.
Two additional runs give:
nr.a = replicate(10^6, sum(sample(deck,5)))
mean(nr.a==4)
[1] 1.9e-05
nr.a = replicate(10^6, sum(sample(deck,5)))
mean(nr.a==4)
[1] 2e-05

My simulations are as accurate as expected, but are not giving
useful approximations. If I had the patience to do simulations
with 10 or 100 million iterations, I would do somewhat better,
but the probability of getting 4 aces in a poker hand is not
something I would choose to simulate.
I don't know the probability of a Tenhou in Mahjong, but I
suspect you may have the same kind of difficulty trying to approximate
that probability by simulation.
Note: A Google search fetched the sentence: "A quick division of the two values shows that the probability of Tenhou or Chihou is around 3.982 × 10–12 or somewhere around 1 in two trillion." (Like anything on the Internet, it may or may not be true.) Presumably this is from one
of several articles listed. Perhaps this one.
