The zoning of massive orebodies at Cu–Ni sulfide deposits such as Noril’sk and Sudbury is commonly explained by fractional crystallization of magmatic sulfide melt. On the theoretical description of fractionation of its components, the results of mineralogical studies of orebodies are usually interpreted using the Rayleigh equation or its modification. But this equation is not applicable to describe crystallization of multicomponent melt and cocrystallization of several phases. In this work we present strict equations describing the distribution of components in a directly crystallized sample. We analyzed the influence of phase reactions on the successive formation of phases during crystallization and on the formation of primary zoning in the sample. This approach permits one to compute the component distribution curves and the crystallization paths by the quantitative phase diagram model. An experimental study of fractionation in the systems Fe–Ni–S and Cu–Fe–S was carried out. They can be regarded as systems modeling the formation of Ni- or Cu-rich sulfide ores. Such studies also yield qualitative and quantitative information about the phase diagrams of geochemical systems. We demonstrated that directed crystallization can be applied to determine the equations of phase reactions and the dependence of partition coefficients on the melt composition and to construct the paths of crystallization and evolution of the tie-line position during one-phase and cotectic crystallization. By the example of the system Fe–Ni–S, all possible types of sample zoning after fractional crystallization are shown. The main regularities of fractionation have been formulated, which are also applicable to multicomponent systems, e.g., Cu–Fe–Ni–S, which is widely used on the modeling of formation of zonal Cu–Ni sulfide ores.