I read that the implementations of Inverse-normal cumulative distribution function (CDF) /quantile / ppf in R, Python (scipy) and Excel give similar results. However, I can't find the very formulae which ensure this accuracy. Can you lead me to such functions?
What I found to be the best solution is this: Epure, Soranzo  For F(x) when x(0;1) works great. The values for the quantile function (F-1) are needed to be 0.5 ≤ p < 1. However, I found the idea in other papers, that when p is below 0.5, you need to take 1-p and negate ( *-1 ) the whole expression. Looks like F-1(p) when p<0.5, then -(F-1(1-p)) does the job. From the test with p(0,0001;0,9999) graphically Excel's formulae and this approximation look identical. The maximum absolute difference (error) between the two arrays of results is 0,08114, and abs(diff^2) is 0,5969. Running the test I needed (utilizing those formulae) for 10 values the maximum difference is 0,0039878. So I think I will stick with this method.