We propose a theoretical framework for analyzing two-body nonleptonic D meson decays, based on the factorization of short-distance (long-distance) dynamics into Wilson coefficients (hadronic matrix elements of four-fermion operators). The parametrization of hadronic matrix elements in terms of several nonperturbative quantities is demonstrated for the D -> PP decays, P denoting a pseudoscalar meson. We consider the evolution of Wilson coefficients with energy release in individual decay modes, and the Glauber strong phase associated with the pion in nonfactorizable annihilation amplitudes, that is attributed to the unique role of the pion as a Nambu-Goldstone boson and a quark-antiquark bound state simultaneously. The above inputs improve the global fit to the branching ratios involving the eta' meson and resolve the long-standing puzzle from the D-0 -> pi(+) pi(-) and D-0 -> K+ K- branching ratios, respectively. Combining short-distance dynamics associated with penguin operators and the hadronic parameters determined from the global fit to branching ratios, we predict direct CP asymmetries, to which the quark loops and the scalar penguin annihilation give dominant contributions. In particular, we predict Delta A(CP) equivalent to A(CP()K(+) K-) - A(CP) (pi(+) pi(-)) = -1.00 x 10(-3), lower than the LHCb and CDF data.