Large scale density-functional theory calculations have been performed to understand the spatial distribution and magnetic coupling of Cr-doped GaN, in which exhaustive structural and magnetic configurations have been investigated by doping of up to five Cr atoms in large supercells. Our results provide direct evidence that the distribution of the doped magnetic ions is neither homogeneous nor random as widely assumed previously. Rather, under both Ga-rich and N-rich growth conditions, the Cr atoms have a strong tendency to form substitutional, embedded clusters with short-range magnetic interactions maintaining the wurtzite structure. Significantly, while the ferromagnetic state is favored for pair doping, for more than two-Cr-atom clustering configurations, states containing antiferromagnetic or ferrimagnetic coupling with net spins in the range of 0.06-1.47µB/Cr are preferred. The formation of embedded clusters leads to notable local structural distortions and considerable magnetic moments on the Cr-bonded N atoms. Also importantly, the electrical properties (metallic, half-metallic, or semiconducting) are found to strongly depend on the dopant concentration. We propose a picture where various cluster configurations coexist and the statistical distribution and associated magnetism depend sensitively on sample growth details. The results obtained are in agreement with recent experiments. Such a view can explain many hitherto puzzling experimental observations, e.g., the much lower value of the measured mean saturation magnetic moment on Cr as compared to the theoretically predicted value for the isolated dopants; the anomalous lattice constant change in relation to the dopant concentration and temperature; and the strong dependence of the magnetization on the Cr concentration, growth temperature, and annealing. We find a similar behavior for Mn in GaN and Cr and Mn in AlN and argue that such a scenario may also hold for other dilute magnetic semiconductor systems.