We give algorithms for approximating the partition function $Z_N$ of the Sherrington–Kirkpatrick (SK) model and finite mixed Ising $p$-spin models with $p\ge 3$. The algorithms apply at a fixed margin below the second-moment threshold $β_{\mathrm f}$ with homogeneous weak external fields $|h|\leqρ_0N^{-1/4}$, for fixed $ρ_0>0$. For any accuracy $\varepsilon\in(0,1)$ and failure parameter $δ\in(0,1)$, they approximate $\log Z_N$ to additive error $\varepsilon$. The number of exact operations is polynomial jointly in $N$, $\varepsilon^{-1}$, and $δ^{-1}$. The disorder variables have i.i.d symmetric distribution with unit variance. The failure probability is at most $δ$ plus an explicit disorder tail bound. This tail bound tends to zero under a finite fourth moment condition. Our approach builds on the cluster expansions of Aizenman–Lebowitz–Ruelle \cite{ALR87} and Dey–Wu \cites{DW23,DW26}. For higher-order mixtures, we retain hypergraphs using few vertices and evaluate their total weight deterministically through local spin averages and inclusion–exclusion. For SK, we impose cutoffs on the number of vertices of degree at least three and on the total number of edges. We estimate the retained sum using randomized color coding \cite{AYZ95}. A lower tail control based on reverse hypercontractivity \cite{MOS13} converts the mean-square approximation of the normalized partition function into logarithmic accuracy.