The exploration of quantum phenomena in complex materials such as moir\'e superlattices is limited by the O(N^3) scaling of conventional electronic structure methods. Here we introduce a high-performance tight-binding framework that reduces the complexity to O(N^1.5) by transforming the Hamiltonian into a real symmetric form and combining Sylvester's inertia law with LDL decomposition. This approach enables efficient band structure calculations for large systems: solving magic angle twisted bilayer graphene in minutes on a laptop and scaling to 1.5 million atoms within days on a workstation. We apply it to the previously inaccessible ultra-low twist-angle regime (less than 0.16 degree) with mechanical strain relaxation and find robust flat bands persisting down to 0.09 degree. Our framework bridges density functional theory accuracy with large-scale quantum simulation, opening a route to systematic data-driven exploration of mesoscale quantum materials.