With growing popularity, unmanned aerial vehicles (UAVs) are pivotally extending conventional terrestrial Internet of Things (IoT) into the sky. To enable high-performance two-way communications of UAVs with their ground pilots/users, cellular network-connected UAV has drawn significant interests recently. Among others, an important issue is whether the existing cellular network, designed mainly for terrestrial users, is also able to effectively cover the new UAV users in the three-dimensional (3D) space for both uplink and downlink communications. Such 3D coverage analysis is challenging due to the unique air-ground channel characteristics, the resulted interference issue with terrestrial communication, and the non-uniform 3D antenna gain pattern of ground base station (GBS) in practice. Particularly, high-altitude UAV often possesses a high probability of line-of-sight (LoS) channels with a large number of GBSs, while their random binary (LoS/Non-LoS) channel states and (on/off) activities give rise to exponentially large number of discrete UAV-GBS association/interference states, rendering coverage analysis more difficult. This paper presents a new 3D system model to incorporate UAV users and proposes an analytical framework to characterize their uplink/downlink 3D coverage performance. To tackle the above exponential complexity, we introduce a generalized Poisson multinomial (GPM) distribution to model the discrete interference states, and a novel lattice approximation (LA) technique to approximate the non-lattice GPM variable and obtain the interference distribution efficiently with high accuracy. The 3D coverage analysis is validated by extensive numerical results, which also show effects of key system parameters such as cell loading factor, GBS antenna downtilt, UAV altitude and antenna beamwidth.