  The decay B_s -> mu^+ mu^- in the Standard Model is generated by the
well-known W-box and Z-penguin diagrams that give rise to an effective
quark-lepton operator Q_A at low energies. We compute QCD corrections of order
alpha_s^2 to its Wilson coefficient C_A. It requires performing three-loop
matching between the full and effective theories. Including the new corrections
makes C_A more stable with respect to the matching scale mu_0 at which the
top-quark mass and alpha_s are renormalized. The corresponding uncertainty in
|C_A|^2 gets reduced from around 1.8% to less than 0.2%. Our results are
directly applicable to all the B_{s(d)} -> l^+ l^- decay modes.
