In this paper we present a novel approach towards variance reduction for discretised diffusion processes. The proposed approach involves specially constructed control variates and allows for a significant reduction in the variance for the terminal functionals. In this way the complexity order of the standard Monte Carlo algorithm ($\varepsilon^{-3}$) can be reduced down to $\varepsilon^{-2}\sqrt{\left|\log(\varepsilon)\right|}$ in case of the Euler scheme with $\varepsilon$ being the precision to be achieved. These theoretical results are illustrated by several numerical examples.