A surface integral equation (SIE) formulation under the magneto-quasi-static assumption is proposed to efficiently and accurately model arbitrarily shaped interconnects in packages. Through decently transferring all electromagnetic quantities into circuit elements, the loop analysis is used to carefully construct matrix equations with an independent and complete set of unknowns based on graph theory. In addition, an efficient preconditioner is developed, and the proposed formulation is accelerated by the pre-corrected fast Fourier transform (pFFT). Four practical examples, including a rectangular metallic interconnect, bonding wire arrays, interconnects in a real-life circuit, and the power distribution network (PDN) used in packages, are carried out to validate its accuracy, efficiency, and scalability. Results show that the proposed formulation is accurate, efficient, and flexible to model complex interconnects in packages.