Boolean networks have been used successfully in modeling biological networks and provide a good framework for theoretical analysis. However, the analysis of large networks is not trivial. In order to simplify the analysis of such networks, several model reduction algorithms have been proposed; however, it is not clear if such algorithms scale well with respect to the number of nodes. The goal of this paper is to propose and implement an algorithm for the reduction of AND-NOT network models for the purpose of steady state computation. Our method of network reduction is the use of "steady state approximations" that do not change the number of steady states. Our algorithm is designed to work at the wiring diagram level without the need to evaluate or simplify Boolean functions. Also, our implementation of the algorithm takes advantage of the sparsity typical of discrete models of biological systems. The main features of our algorithm are that it works at the wiring diagram level, it runs in polynomial time, and it preserves the number of steady states. We used our results to study AND-NOT network models of gene networks and showed that our algorithm greatly simplifies steady state analysis. Furthermore, our algorithm can handle sparse AND-NOT networks with up to 1000000 nodes.