Let $X_1,\dots, X_n$ be i.i.d. random variables sampled from a normal distribution $N(μ,Σ)$ in ${\mathbb R}^d$ with unknown parameter $θ=(μ,Σ)\in Θ:={\mathbb R}^d\times {\mathcal C}_+^d,$ where ${\mathcal C}_+^d$ is the cone of positively definite covariance operators in ${\mathbb R}^d.$ Given a smooth functional $f:Θ\mapsto {\mathbb R}^1,$ the goal is to estimate $f(θ)$ based on $X_1,\dots, X_n.$ Let $$ Θ(a;d):={\mathbb R}^d\times \Bigl\{Σ\in {\mathcal C}_+^d: σ(Σ)\subset [1/a, a]\Bigr\}, a\geq 1, $$ where $σ(Σ)$ is the spectrum of covariance $Σ.$ Let $\hat θ:=(\hat μ, \hat Σ),$ where $\hat μ$ is the sample mean and $\hat Σ$ is the sample covariance, based on the observations $X_1,\dots, X_n.$ For an arbitrary functional $f\in C^s(Θ),$ $s=k+1+ρ, k\geq 0, ρ\in (0,1],$ we define a functional $f_k:Θ\mapsto {\mathbb R}$ such that \begin{align*} & \sup_{θ\in Θ(a;d)}\|f_k(\hat θ)-f(θ)\|_{L_2({\mathbb P}_θ)} \lesssim_{s, β} \|f\|_{C^{s}(Θ)} \biggr[\biggl(\frac{a}{\sqrt{n}} \bigvee a^{βs}\biggl(\sqrt{\frac{d}{n}}\biggr)^{s} \biggr)\wedge 1\biggr], \end{align*} where $β=1$ for $k=0$ and $β>s-1$ is arbitrary for $k\geq 1.$ This error rate is minimax optimal and similar bounds hold for more general loss functions. If $d=d_n\leq n^α$ for some $α\in (0,1)$ and $s\geq \frac{1}{1-α},$ the rate becomes $O(n^{-1/2}).$ Moreover, for $s>\frac{1}{1-α},$ the estimators $f_k(\hat θ)$ is shown to be asymptotically efficient. The crucial part of the construction of estimator $f_k(\hat θ)$ is a bias reduction method studied in the paper for more general statistical models than normal.