In this work, we determine the short-distance coefficients for gamma inclusive decay into a charm pair through relative order v(4) within the framework of nonrelativistic QCD ( NRQCD) factorization formula. The short-distance coefficient of the order-v(4) color-singlet NRQCD matrix element is obtained through matching the decay rate of b (b) over bar(S-3(1))-> c (c) over bar gg in full QCD to that in NRQCD. The double and single IR divergences appearing in the decay rate are exactly canceled through the next-to-next-to-leading-order renormalization of the operator O(S-3(1)) and the next-to-leading-order renormalization of the operators O(P-3(J)). To investigate the convergence of the relativistic expansion arising from the color-singlet contributions, we study the ratios of the order-v(2) and -v(4) color-singlet short-distance coefficients to the leading-order one. Our results indicate that though the order-v(4) color-singlet short-distance coefficient is quite large, the relativistic expansion for the color-singlet contributions in the process gamma -> c (c) over bar + X converges well due to a small value of v. In addition, we extrapolate the value of the mass ratio of the charm quark to the bottom quark, and find the relativistic corrections rise quickly with increase of the mass ratio.