Superhard cubic boron-carbon nitride (c-BNC) in terms of bulk modulus along a composition range of (BN)((1-x))(C(2))(x) (0 <= x <= 1) is systematically explored by Monte Carlo simulations and cluster expansion techniques based on first-principles calculation. Bulk moduli for the c-BNC ordered structures are reasonably expanded up to quadruplet clusters, indicating that dependence of the bulk modulus on atomic arrangements is not simply attributed to pairwise interactions. A negative correlation can be seen between bulk modulus and formation energies, which is consistent with previous theoretical works. Monte Carlo simulation reveals that all the ordered structures with the highest bulk modulus at each composition exhibit a strong preference of neighboring B-N and C-C atoms, which is consistent with the bond counting rule previously suggested. A composition dependence of these ordered structures can be observed. At a BN-rich composition of x = 0.25, C atoms form a nearest-neighbor network with a hexagonal cluster shape, while at equiatomic and diamond-rich compositions of x = 0.5 and 0.75, B and N atoms form nearest-neighbor networks with a planar shape. At x = 0.875, c-BNC ordered structure with neighboring B and N atoms forming a stereoscopic shape exhibit the highest bulk modulus of 459.3 GPa, which is similar to 0.6% smaller than that of diamond.