We propose a novel parallel numerical algorithm for calculating the smallest eigenvalues of highly ill-conditioned matrices. It is based on the LDLT decomposition and involves finding a k × k sub-matrix of the inverse of the original N × N Hankel matrix H −1 N . The computation involves extremely high precision arithmetic, message passing interface, and shared memory parallelisation. We demonstrate that this approach achieves good scalability on a high performance computing cluster (HPCC) which constitute a major improvement of the earlier approaches. We use this method to study a family of Hankel matrices generated by the weight w(x) = e −x β , supported on [0, ∞) and β > 0. Such weight generates Hankel determinant, a fundamental object in random matrix theory. In the situation where β > 1/2, the smallest eigenvalue tend to 0, exponentially fast as N gets large. If β < 1/2, the situation where the classical moment problem is indeterminate, the smallest eigenvalue is bounded from below by a positive number for all N , including infinity. If β = 1/2, it is conjectured that the smallest eigenvalue tends to 0 algebraically, with a precise exponent. The algorithm run on the HPCC producing fantastic match between the theoretical value of 2/π and the numerical result.