The structural and thermodynamical properties of GaxIn1-xP solid solutions are studied using Monte Carlo simulations based on density-functional perturbation theory. The free-standing alloy displays a model random behavior, with a miscibility gap below approximately 820 K. In samples grown on a GaAs (001) substrate, the resulting epitaxial strain stabilizes various ordered phases which give rise to a rich phase diagram below approximately 230 K. Among these phases, the chalcopyrite is dominant for x = 0.5, and two new, more complex, stable phases have been identified at different concentrations.