Permutation procedures are essential for hypothesis testing when the distributional assumptions about the considered test statistic are not met or unknown, but are challenging in scenarios with limited permutations, such as complex biomedical studies. P-values may either be zero, making multiple testing adjustment problematic, or too large to remain significant after adjustment. A common heuristic solution is to approximate extreme p-values by fitting a Generalized Pareto Distribution (GPD) to the tail of the distribution of the permutation test statistics. In practice, an estimated negative shape parameter combined with extreme test statistics can again result in zero p-values. To address this issue, we present a comprehensive workflow for accurate permutation p-value approximation that fits a constrained GPD and strictly avoids zero p-values. We also propose new methods that address the challenges of determining an optimal GPD threshold and correcting for multiple testing. Through extensive simulations, our approach demonstrates considerably higher accuracy than existing methods. The computational framework will be available as the open-source R package “permAprox”.