MooseBoys is right when he says
integrate the hemisphere about the normal
And to some extend you are also right when you say that a low mipmap is doing that, but not 100% mathematically correct.
What you need to have physically speaking, is a convolution with a cosine lobe in linear space. What you get with the mipmap sampling, is a convolution with a box window in gamma space.
Now, iI'm not saying it will give completely different results, its a good way to approximate the physically accurate model, with hardware support and ease of implementation ! You can keep it this way if you want. If you want super crazy physically correct precision, then I suggest you download AMD CubeMapGen modified by Sebastien Lagarde, and check the GUI box "irradiance", then you can export the result as .dds and instead of sampling the mipmap you can sample the irradiance map.
Now, you need to take the Sun OFF your skycube before you generate the irradiance map, because you WILL, like you feared, apply the sun 2 times !
To do that, just down-tone the sun corolla area to a brightness that has similar values with the surrounding (e.g. manually with Gimp). Its enough. Physically speaking, a sun as a directional light should get values of about 20/25, with an exposure of 1, it would simply saturate everything in an SDR renderer. If you have a range-compressor tonemapper you could get away with this kind of values. If you have an SDR renderer, you want to pre-tone-map manually the intensity of your directional light, to something near 1 so that you don't saturate things. But in that case, the contribution from your irradiance map will be unusually high, it would make sense to premultiply it with the same down-tone filter after its baked by cubemapgen. (again, in a full HDR pipeline, you care not about those things, its handled in real time).
EDIT:
Here is the values we find in a classic HDR sky map:
