This is going to be quite formal, but a nice analogy to electrostatics.
The ideal electric dipole
If you take two opposite point charges and at constant (electric) dipole moment $\vec{p}$ reduce their distance to zero you get the charge density
$$
\rho = - \vec{p} \cdot \nabla \delta
$$
with the delta distribution $\delta$. To compute the field created by this charge distributions we need three ingredients:
- Gauß' law
$$
\nabla \cdot \vec{D} = \rho\,,
$$
which relates the dielectric displacement $\vec{D}$ to the charge density,
- Faraday's law (static version here)
$$
\nabla \times \vec{E} = 0\,,
$$
which yields the fact that the electric field can be described as the gradient of a potential, $\vec{E} = - \nabla \phi$, and
- a material law for the vacuum,
$$
\vec{D} = \varepsilon_0 \vec{E}\,.
$$
Together they yield the Poisson equation
$$
-\Delta \phi = \frac{\rho}{\varepsilon_0}\,.
$$
The operator on the left has the fundamental solution
$$
G = \frac{1}{4 \pi r}\,,
$$
i.e.
$$
- \Delta G = \delta
$$
and for a general charge the solution to the Poisson equation is
$$
\phi = \frac{1}{\varepsilon_0} \rho \ast G
$$
where the asterisk denotes the convolution. This way, for example, we get the Coulomb potential for the point charge.
Using the charge distribution of the ideal dipole we get the potential
$$
\phi = \frac{1}{\varepsilon_0} \vec{p} \cdot \nabla G
$$
and the electric field
$$
\vec{E} = - \nabla \frac{1}{\varepsilon_0} \vec{p} \cdot \nabla G\,.
$$
The ideal magnetic dipole
For the ideal magnetic dipole it's exactly the same, only a little different. ;)
Take a current loop of given magnetic dipole moment $\vec{m}$ and reduce its size to zero. The distributional limit will be the current density
$$
\vec{j} = - \vec{m} \times \nabla \delta\,.
$$
To compute the field created by this current distribution we need three ingredients:
- Ampère's law (magnetostatic version)
$$
\nabla \times \vec{H} = \vec{j}\,,
$$
which relates the magnetic field to the current density,
- the equation
$$
\nabla \cdot \vec{B} = 0\,,
$$
which ensures the existence of a vector potential such that $\vec{B} = \nabla \times \vec{A}$, and
- you guessed it, a material law of the vacuum
$$
\vec{B} = \mu_0 \vec{H}\,.
$$
Together they yield the Poisson equation for the vector potential
$$
-\Delta \vec{A} = \mu_0 \vec{j}
$$
We already know how to solve it - the vector potential will be the convolution
$$
\vec{A} = - \mu_0 \vec{j} \ast G\,.
$$
And so, in complete analogy to electrostatics we here find for the magnetic dipole
$$
\vec{A} = \mu_0 \vec{m} \times \nabla G\,,
$$
and the field of the magnetic dipole is
$$
\vec{B} = \nabla \times \left(\mu_0 \vec{m} \times \nabla G\right) \\
= \mu_0 \vec{m} \delta - \mu_0 \nabla \vec{m} \cdot \nabla G \,.
$$
Explicitly calculating the derivatives of the fundamental solution now yields what you are interested in. For including the delta terms you might want to look at Frahm, C.P. Some novel delta-function identities. Am. J. Phys. 1983. You don't have to if you just want to get the classical part of the solution...