jStat beta.cdf/ibeta wrong (even negative) for a,b > ~1e5: CF capped at 100
jStat 1.9.6 (latest; also used by @formulajs BETA.DIST/F.DIST) returns jStat.beta.cdf(0.5, 1e7, 1e7) = -0.17 and beta.cdf(0.5, 1e6, 1e6) = 0.44; both should be 0.5. jStat.betacf runs the Lentz continued fraction for at most 100 iterations and uses the result whether or not it converged. Near x ~ (a+1)/(a+b+2) the fraction needs O(sqrt(max(a,b))) terms. The error sits in a narrow band around the mean (at +-2 sd it is still ~1e-9), so spot checks at the tails miss it. This bites Bayesian A/B-test code with Beta(s+1, n-s+1) posteriors at n >= 1e6. Fix: maxIter = max(100, ceil(10*sqrt(max(a,b)))), which gives <5e-7 up to a=b=1e9. Or use scipy/@stdlib betainc as a cross-check.
jStat: beta.cdf(0.5, 1e7, 1e7) = -0.1736 and (0.5, 1e6, 1e6) = 0.4404, while 1e5 gives 0.49989 and 1e4 gives 0.50000006. At +2 sd (a=b=1e6) the relative error is 2.4e-9. I patched the loop bound in betacf to your max(100, ceil(10*sqrt(max(a,b)))). Absolute errors at x=0.5 then dropped to 6-7e-8 (1e5 and 1e6), 6.5e-8 (1e7) and 1.4e-7 (1e9), and beta.cdf(0.3, 3e6, 7e6) to 7e-8. The +2 sd value did not change, though (still 2.4e-9), so the iteration cap isn't what limits accuracy there. A floor of roughly 1e-9 to 1e-7 remains with the patch, probably from the prefactor computed via large log-gamma values; I haven't traced it. For posterior probabilities near 0.5 that is enough, but don't compare jStat to SciPy at rtol=1e-12 even after patching.
Also applies to HyperFormula 3.5.0, which bundles a copy of jStat: =BETA.DIST(0.5,10000000,10000000,TRUE()) returns -0.1736 and =F.DIST(1,1000000,1000000,TRUE()) returns 0.4826 (both should be 0.5). Patching the npm jstat package does not fix HyperFormula's copy.