A set of compact and stable formulas are developed for the computation of triaxial electromagnetic (EM) logging response in multilayered formations with biaxial anisotropy. These formulas are based on a new global propagator matrix method which derives the amplitudes of upgoing/downgoing waves at the boundaries rather than the reflection/transmission coefficients of each layer. The kernel of the new method is assembling the local continuity equations into a global linear system, which has addressed the notorious overflow problems commonly encountered by traditional propagator matrix approach. The triaxial EM logging responses are then calculated via an efficient two-fold inverse Fourier transform (IFT) technique. Two main techniques are developed in order to significantly improve the computation accuracy and speed of IFT. On the one hand, the spectral fields at the first integral quadrant can be used to express the fields at other three quadrants. Consequently, the original two-fold infinite integral is analytically simplified to a semi-infinite one and the computation cost can be reduced almost by three fourths. On the other hand, a two-level subtraction scheme is developed to reduce the integral singularity. The first level subtraction is achieved by separating the primary field from the total fields. The second level process includes subtracting the spectral scattered fields of an equivalent transverse isotropic (TI) layered medium, integrating the residual fields, and computing the spatial fields of TI model. The subtraction on the primary and scattered fields can make the integral convergence increase by more than ten orders of magnitude in favorable condition. Numerical experiments demonstrate that the new semi-analytical formulas are well suited for any dipping angle and arbitrary number of layers, as well as for formations with complete thicknesses. The modeling method is efficient and robust, which can provide a fast simulator for the interpretation of triaxial EM logging data.