Ionic solutions exhibiting multiple association are described within the binding mean spherical approximation (BiMSA). This model is based on the Wertheim formalism, in the framework of the primitive model at the McMillan-Mayer level. The cation and the anion form the various complexes according to stepwise complexation-equilibria. Analytic expressions for the Helmholtz energy, the internal energy, the speciation, and for the osmotic and activity coefficients are given considering a binary solution with an arbitrary number of association sites on one type of ion (polyion) and one site on the ions of opposite sign (counterions). As an alternative, mean field expressions, as developed in SAFT-type theories, are also presented. The result obtained from the latter approximate method exhibits a reasonable agreement with those from BiMSA for the speciation, and a remarkable one for the osmotic coefficient.