How do I properly do gamut mapping for an image?
Wikipedia page Gamut puts it nicely:


  Generally, the color gamut is specified in the hue–saturation plane,
  as a system can usually produce colors over a wide intensity range
  within its color gamut; [...]


But you don't have the general case. Without black for example, you can't produce darker shades of all colors even if they are within the gamut projected onto the XY plane. So you will need to determine the 3-d gamut polyhedron. The gamut is a polyhedron in an intensity-linear (gamma=1.0, linear gamma) RGB or XYZ color space, because intensities of color components mix additively in additive color mixing.

You have additive color mixing in the kind of dithering that you do. If someone looks at the picture from afar, they only see the average intensities of the color components of the pixels. This is not the same as for example average squared intensities, so gamma matters. In dithering, the error calculation and the residual error diffusion arithmetic should be done in an intensity-linear color space. 

3-d gamut mapping is not going to be trivial. First you need to find the convex hull of the set of color coordinate points in your palette, in an intensity-linear color space. Then you need to check if a wanted color is inside or outside the convex hull. If it's outside, you need to map it to the surface of the convex hull, which is (the union of) a set of triangles in 3-d space. Here the easiest mapping might be to find the nearest point on the convex hull.

Better results might be had by using a perceptual color difference metric like CIEDE2000. However, this makes finding the nearest point a non-trivial optimization problem, possibly non-convex. Possibly the perceptual color difference metric would also give better results in dithering, just at the step where you choose from the palette the color to use.