Summary: 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 \leq 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.