By a recursive method numerically exact free energies are calculated for square L×L Ising lattices, with 6<~L<~18, for several kinds of frozen-in bond disorder: (i) bonds ±J with various concentrations of negative bonds; (ii) bonds distributed according to a Gaussian distribution. Ground states of these systems are identified, the response to "ordering fields" is studied, and the correlation function 〈S0SR〉T2 is calculated as a function of temperature for various distances R in the lattice. This correlation is found to decay strongly (presumably exponentially) with increasing R even at temperatures distinctly below the apparent freezing temperature Tf of previous Monte Carlo simulations; this "freezing transition" is hence unambiguously identified as a nonequilibrium effect. However, the correlation length is found to become long ranged at low temperatures, and it is suggested that a phase transition still occurs at T=0; while in the Gaussian model the spin-glass order parameter q(T=0)=1, it is found that q≡0 in the ±J model where rather a power-law decay of correlations 〈S0SR〉T=02 occurs. Performing Monte Carlo simulations for precisely the same systems, the cooling times necessary to reach the true ground states of the system are identified, as well as the simulation times necessary to reach thermal equilibrium for the correlation functions. These times are found to increase so strongly with L that for systems of macroscopic size the correct thermal equilibrium is probably irrelevant for experimental purposes. Rather a statistical mechanics based on the many long-lived metastable states would be required.