We present guaranteed, robust and computable a posteriori error bounds for approximate solutions of the equation ∆∆u + κ²u = f by classical and mixed Ciarlet-Raviart finite element methods. We concentrate on the case when the reaction coefficient κ² is subdomain (finite element) wise constant and chaotically varies between subdomains in the sufficiently wide range. It is proved that the bounds for the classical FEM's are robust with respect to κ ∈ [0, ch⁻²], where c = const and h is the maximal size of finite elements, and possess additional useful features. The coefficients in fronts of two typical norms in their right parts only insignificantly worse than those for κ ≡ const, and the bounds can be calculated without resorting to the equilibration procedures. Besides, they are sharp at least for low order methods, if the testing moments and deflection in their right parts are found by accurate recovery procedures. The technique of derivation of the bounds is based on the approach similar to one used in our preceding papers for simpler problems.