The $0.001$ upper tail critical value for a $\chi^2_1$ is $10.82756617...$
Since the chi-square with 1 df is the square of a normal, if you have more accuracy in z-tables you can square the 0.0005 upper tail critical value for a Z.
The "formula" for the tail area of a $\chi^2_1$ in effect requires the use of either the error function or the incomplete gamma function, which are generally evaluated by computer algorithms. For tail areas of chi-squares in general, the incomplete gamma function is used. For this problem (given the tail area, find the chi-squared value that gives it), you need the inverse incomplete gamma function.
The incomplete gamma function and its inverse function may be evaluated in a variety of ways -- for example, series approximation, through some iterative or recursive method, split into ranges where accurate functional approximations (possibly using ratios of polynomials for example) might be used. It's impossible to guess exactly what methods were used to produce your table.
I used a computer program to find the value to more digits - specifically, in R, I used
(I think this is a couple more digits than the approximation will be accurate to.)