In this paper, a highly accurate algorithm for computation of complex valued Bessel, Neumann and Hankel functions of integer order is given. The algorithm enables the computation of these functions in the entire complex plane with quadruple precision, which can be reduced to double precision. The complex values of the Bessel and Neumann functions of the zeroth and first order can be computed in a special way for small, medium-sized and large arguments in the first quadrant of the complex plane. The mapping of functions from the first quadrant to the other quadrants is described by simple formulas. Bessel and Neumann functions of higher positive integer order can be computed using forward and backward recurrence relations. Two types of Hankel functions are linear combinations of the Bessel and Neumann functions. Bessel, Neumann and Hankel functions of negative integer order are equal to positive order functions up to the sign.