In this paper, we advance an efficient algorithm of electromagnetic spatial dyadic Greens function in a horizontal-layered anisotropic medium through high order window function and so on. First, we use the high order window function to transform the Sommerfeld integrals of dyadic Greens function into integrals along a generalized steep descent path. And we give a new and higher order of power series of expansion expression of Hankel transforms of the window function and an accurate formula of Lommel function so that the dyadic Greens function can be precisely computed either near or far from the transmitter. Then, we divide the integral path into a series of subintervals based on the zeros of Bessel function and compute integral per subinterval using an improved adaptive Gauss quadrature. Furthermore, an algorithm of the finite continued fraction expansion is used to sum up the integrals on each subinterval to greatly accelerate the convergence velocity of the numerical integration along the descent path. Finally, our numerical results validated the efficiency of the algorithm.