The main purpose of this paper is to establish the Hormander-Mihlin type theorem for Fourier multipliers with optimal smoothness on k-parameter Hardy spaces for k >= 3 using the multi-parameter Littlewood-Paley theory. For the sake of convenience and simplicity, we only consider the case k = 3, and the method works for all the cases k >= 3: T(m)f(x(1), x(2), x(3)) = 1/(2 pi)(n1+n2+n3) integral(Rn1 x Rn2 x Rn3) m(xi)(f) over cap(xi)e(2 pi ix.xi)d xi where x = (x(1), x(2), x(3)) is an element of R-n1 x R-n2 x R-n3 and xi=(xi(1), xi(2), xi(3)) is an element of R-n1 x R-n2 x R-n3. One of our main results is the following: Assume that m(xi) is a function on Rn1+n2+n3 satisfying sup(j, k, l is an element of Z) parallel to m(j, k, l) parallel to(W(s1,s2,s3)) < infinity with s(i) > n(i)(1/p - 1/2) for 1 <= i <= 3. Then T (m) is bounded from H-p (R-n1 x R-n2 x R-n3) to H-p(R-n1 x R-n2 x R-n3) for all 0 < p <= 1 and parallel to T-m parallel to H-p -> H-p less than or similar to sup(j, k, l is an element of Z) parallel to m(j, k, l) parallel to (W(s1, s2, s3)). Moreover, the smoothness assumption on s(i) for 1 <= i <= 3 is optimal. Here we have used the notations m(j,k,l) (xi) = m(2 (j)xi(1), 2(k)xi(2), 2(l)xi(3))Psi(xi(1))Psi(xi(2))Psi(xi(3)) and Psi(xi (i) ) is a suitable cut-off function on R-ni for 1 <= i <= 3, and is a three-parameter Sobolev space on R-n1 x R-n2 x R-n3. Because the Fefferman criterion breaks down in three parameters or more, we consider the L-p boundedness of the Littlewood-Paley square function of T(m)f to establish its boundedness on the multi-parameter Hardy spaces.