  The Cervera et al. formula, the best known approximate formula of neutrino
oscillation probability for long-baseline experiments, can be regarded as a
second-order perturbative formula with small expansion parameter epsilon \equiv
Delta m^2_{21} / Delta m^2_{31} \simeq 0.03 under the assumption s_{13} \simeq
epsilon. If theta_{13} is large, as suggested by a candidate nu_{e} event at
T2K as well as the recent global analyses, higher order corrections of s_{13}
to the formula would be needed for better accuracy. We compute the corrections
systematically by formulating a perturbative framework by taking theta_{13} as
s_{13} \sim \sqrt{epsilon} \simeq 0.18, which guarantees its validity in a wide
range of theta_{13} below the Chooz limit. We show on general ground that the
correction terms must be of order epsilon^2. Yet, they nicely fill the mismatch
between the approximate and the exact formulas at low energies and relatively
long baselines. General theorems are derived which serve for better
understanding of delta-dependence of the oscillation probability. Some
interesting implications of the large theta_{13} hypothesis are discussed.
