2 corresponding pure versions are exactly known. This is also the case for the present model under consideration, namely the triangular Ising model, called hereafter as the TrIM, defined as usual by the Hamiltonian where the spin variables s$_{i}$ take on the values − 1 , +1, < ij > indicates summation over all nearest-neighbor pairs of sites, and J > 0 is the ferromagnetic exchange interaction. The TrIM belongs to the same universality class with the corresponding square Ising model (SqIM), sharing the same values of critical exponents and a log- arithmic behavior of the specific heat [60, 61]. Addition- ally, the critical temperature of the model and also the critical amplitude A$_{0}$ of Ferdinand and Fisher’s [62] spe- cific heat’s logarithmic expansion [see also the discussion in Sec. III and Eq. (7)] are exactly known from the early work of Houtappel [63] to be T$_{c}$ = 4 / ln 3 = 3 . 6409 · · · and A$_{0}$ = 0 . 499069 · · · , respectively. Nevertheless, it ap- pears that for the TrIM a verification of the finite-size scaling (FSS) properties of the model using high qual- ity data from MC simulation is still lacking. Thus, the first part of this work is devoted to the investigation of the FSS behavior of the model, especially the estimation of the amplitudes and other relevant coefficients in the specific heat’s logarithmic expansion and also to the es- timation of the critical exponents. In this sense, the aim of this first part is twofold: First, to provide the first detailed FSS analysis of the pure model and, second, to present a concrete reliability test of the proposed numer- ical scheme. Our main focus, on the other hand, is the case with bond disorder given by the bimodal distribution so that r reflects the strength of the bond randomness and we fix 2 k$_{B}$/ ( J$_{1}$ + J$_{2}$ ) = 1 to set the temperature scale. The value of the disorder strength considered throughout this paper is r = 1 / 3. The resulting quenched disor- dered (random-bond) version of the Hamiltonian defined in Eq. (1) reads now as and will be referred in the sequel as the random-bond triangular Ising model (RBTrIM). The corresponding random-bond SqIM will be denoted hereafter respectively as RBSqIM. The model on the square lattice has the ad- vantage that the critical temperature is exactly known as a function of the disorder strength r by duality re- lations [64]. For the RBTrIM there exist only several approximations for the critical frontier of the site- and bond-diluted cases, obtained via renormalization-group techniques [65] and, to our knowledge, a study of the critical behavior of the model is lacking. The rest of the paper is laid out as follows: In Sec. II we present the necessary simulation details of our numer- ical scheme. In Sec. III we discuss the FSS behavior of the pure model, testing with our high accuracy numeri- cal data the exact expansion of the critical specific heat. Then, in Sec. IV we present a detailed FSS analysis for the random version of the model, including - apart from the classical FSS techniques - concepts from the scaling theory of disordered systems. Our results and the rele- vant discussion clearly favors the scenario of strong uni- versality in marginal disordered systems. Finally, Sec. V summarizes our conclusions. Resorting to large scale MC simulations is often nec- essary [66], especially for the study of the critical be- havior of disordered systems. It is also well known [67] that for such complex systems traditional methods be- come very inefficient and that in the last few years sev- eral sophisticated algorithms, some of them are based on entropic iterative schemes, have been proven to be very effective. The present numerical study of the RB- TrIM will be carried out by applying our recent and ef- ficient entropic scheme [59, 68, 69]. In this approach we follow a two-stage strategy of a restricted entropic sam- pling, which is described in our study of random-bond Ising models (RBIM) in 2d [59] and is very similar to the one applied also in our numerical approach to the 3d random-field Ising model (RFIM) [69]. In these pa- pers, we have presented in detail the various sophisti- cated routes used for the restriction of the energy sub- space and the implementation of the Wang-Landau (WL) algorithm [70]. Further details and an up to date imple- mentation of this approach, especially for the study of disordered systems, is provided in our recent paper on the universality aspects of the pure and random-bond 2d Blume-Capel model [71]. We do not wish to reproduce here the details of our two-stage implementation and the practice followed in our scheme for improving accuracy by repeating the sim- ulations. However, we should like to include a brief dis- cussion on the approximate nature of the WL method. The usual WL recursion proceeds by modifying the den- sity of states G ( E ) according to the rule G ( E ) → fG ( E ) and initially one chooses G ( E ) = 1 and f = f$_{0}$ = e . Once the accumulative energy histogram is sufficiently flat, the modification factor f is redefined as: f$_{j}$$_{+1}$ = √ f$_{j}$ , with j = 0 , 1 , . . . and the energy histogram reset to zero until f is very close to unity (i.e. f = e 10 − 8 ≈ 1 . 000 000 01). As has been reported by many authors in the study of sev- eral models, once f is close enough to unity, systematic deviations become negligible. However, the WL recur- sion violates the detailed balance from the early stages of H = − J ∑ <ij> s$_{i}$s$_{j}$, (1) P ( J$_{ij}$ ) = 1 2 [ δ ( J$_{ij}$ − J$_{1}$ ) + δ ( J$_{ij}$ − J$_{2}$ )] ; (2) J$_{1}$ + J$_{2}$ 2 = 1 ; J$_{1}$ > J$_{2}$ > 0 ; r = J$_{2}$ J$_{1}$ , H = − ∑ <ij> J$_{ij}$s$_{i}$s$_{j}$ (3) II. SIMULATION DETAILS