A common approach for studying a solid solution or disordered system within a periodic ab-initio framework is to create a supercell in which a certain amount of target elements is substituted with other ones. The key to generating supercells is determining how to eliminate symmetry-equivalent structures from the large number of substitution patterns. Although the total number of substitutions is on the order of trillions, only symmetry-inequivalent atomic substitution patterns need to be identified, and their number is far smaller than the total. A straightforward solution would be to classify them after determining all possible patterns, but it is redundant and practically unfeasible. Therefore, to alleviate this drawback, we developed a new formalism based on the {\it canonical augmentation}, and successfully applied it to the atomic substitution problem. Our developed \verb|python| software package, which is called \textsc{SHRY} (\underline{S}uite for \underline{H}igh-th\underline{r}oughput generation of models with atomic substitutions implemented by p\underline{y}thon), enables us to pick up only symmetry-inequivalent structures from the vast number of candidates very efficiently. We demonstrate that the computational time required by our algorithm to find $N$ symmetry-inequivalent structures scales {\it linearly} with $N$ up to $\sim 10^9$. This is the best scaling for such problems.