We present an auto-differentiable spectral modeling of exoplanets and brown dwarfs. This model enables a fully Bayesian inference of the high-dispersion data to fit the ab initio line-by-line spectral computation to the observed spectrum by combining it with the Hamiltonian Monte Carlo in recent probabilistic programming languages. An open source code, exojax, developed in this study, was written in Python using the GPU/TPU compatible package for automatic differentiation and accelerated linear algebra, JAX (Bradbury et al. 2018). We validated the model by comparing it with existing opacity calculators and a radiative transfer code and found reasonable agreements of the output. As a demonstration, we analyzed the high-dispersion spectrum of a nearby brown dwarf, Luhman 16 A and found that a model including water, carbon monoxide, and $mathrm{H_2/He}$ collision induced absorption was well fitted to the observed spectrum ($R=10^5$ and $2.28-2.30 mumathrm{m}$). As a result, we found that $T_0 = 1295 pm 14 mathrm{K}$ at 1 bar and $mathrm{C/O} = 0.62 pm 0.01$, which is slightly higher than the solar value. This work demonstrates the potential of full Bayesian analysis of brown dwarfs and exoplanets as observed by high-dispersion spectrographs and also directly-imaged exoplanets as observed by high-dispersion coronagraphy.