  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.
