As a preparation step to compute Jacobian elliptic functions efficiently, we created a fast method to calculate the complete elliptic integral of the first and second kinds, K (m) and E(m), for the standard domain of the elliptic parameter, 0 < m < 1. For the case 0 < m < 0.9, the method utilizes 10 pairs of approximate polynomials of the order of 9-19 obtained by truncating Taylor series expansions of the integrals. Otherwise, the associate integrals, K (1 − m) and E(1 − m), are first computed by a pair of the approximate polynomials and then transformed to K (m) and E(m) by means of Jacobi's nome, q, and Legendre's identity relation. In average, the new method runs more-than-twice faster than the existing methods including Cody's Chebyshev polynomial approximation of Hastings type and Innes' formulation based on q-series expansions. Next, we invented a fast procedure to compute simultaneously three Jacobian elliptic functions, sn(u|m), cn(u|m), and dn(u|m), by repeated usage of the double argument formulae starting from the Maclaurin series expansions with respect to the elliptic argument, u, after its domain is reduced to the standard range, 0 ≤ u < K (m)/4, with the help of the new method to compute K (m). The new procedure is 25-70% faster than the methods based on the Gauss transformation such as Bulirsch's algorithm, sncndn, quoted in the Numerical Recipes even if the acceleration of computation of K (m) is not taken into account.