High-Entropy Alloys (HEAs) exhibit complex atomic interactions, with short-range order (SRO) playing critical role in determining their properties. Traditional methods, such as Monte Carlo generator of Special Quasirandom Structures within the Alloy Theoretic Automated Toolkit (ATAT-mcsqs), Super-Cell Random APproximates (SCRAPs), and hybrid Monte Carlo-Molecular Dynamics (MC-MD)-are often hindered by limited system sizes and high computational costs. In response, we introduce PyHEA, a Python-based toolkit with high-performance C++ core that leverages global and local search algorithms, incremental SRO computations, and GPU acceleration for unprecedented efficiency. When constructing random HEAs, PyHEA achieves speedups exceeding 333,000x and 13,900x over ATAT-mcsqs and SCRAPs, respectively, while maintaining high accuracy. PyHEA also offers a flexible workflow that allows users to incorporate target SRO values from external simulations (e.g., LAMMPS or density functional theory (DFT)), thereby enabling more realistic and customizable HEA models. As a proof of concept, PyHEA successfully replicated literature results 256,000-atom Fe-Mn-Cr-Co system within minutes-an order-of-magnitude improvement over hybrid MC approaches. This dramatic acceleration opens new possibilities for bridging theoretical insights and practical applications, paving the way for the efficient design of next-generation HEAs.