In this paper, a number theoretic method is introduced to calculate the high order correlation functions of the Ulam map. In this method, the calculation is firstly transformed into solving a class of exponential Diophantine equations with variable coefficients. Then thses equations are simplified to the Diophantine equations with strictly monotonic exponentials. Finally, the equations are solved by means of an order reduction method. As an application, the first five order correlation functions of the Ulam map are calculated.