We propose a ``blending" algorithm that projects the all-to-all fermion propagator onto spatial low-frequency modes (LFM) combined with a stochastic estimate of spatial high-frequency modes (SHFM) at each time slice. This approach enables the calculation of arbitrary-point correlation functions for arbitrary hadron states in strongly interacting quantum field theories (QFT) with fermions, such as quantum chromodynamics (QCD). Specifically, LFM allows the construction of spatially extended hadron states below a certain energy threshold by diagonalizing multi-fermion interpolation fields. Meanwhile, the local interactions required for N-point correlation functions in QFT can be approximated in an unbiased manner through a reweighted summation of both LFM and SHFM contributions. To demonstrate the efficiency of this algorithm, we obtained {\color{black} $g_A^u=0.895(15)$, $g_A^d=-0.338(15)$, $g_A^s=-0.0245(72)$, $g_A^{u+d+s}=0.533(28)$ and $g_A^{u-d}=1.2339(43)$ } for nucleon at $m_{\pi}=300$ MeV and $a=0.077$ fm using 40 configurations. The consistency check of the pion electric form factor and charge radius derived from 3-point and 4-point correlation functions is also provided.
Comment: 20 pages, 15 figures, some of statements updated