I think that after my comment I should provide some sort of answer to avoid appearing as a pedantic mathematician (mathematical physicist). Actually my PhD is in theoretical physics and I'm proud of it ...
First of all let us get rid of issues with domains and of the subtle distinctions between Hermitian, symmetric, selfadjoint and so on.
If we assume to deal with operators $A,B$ defined everywhere on a Hilbert space ${\cal H}$, then $[A,B]$ is well defined and all the considered operators are also bounded (continuous) as a consequence of a well-known theorem on everywhere defined selfadjoint operators.
In this case, also the distinction between Hermitian, symmetric and selfadjoint becomes immaterial.
In this case, the spectra $\sigma(A)$ and $\sigma(B)$ are closed bounded subsets of $\mathbb{R}$. I stress that the spectra may include both continuous and point parts.
Unfortunately, the number of situations where this happens in physics is really negligible. The fundamental reason is that the spectrum is made of the observed values of the observables $A$ and $B$ and bounded spectra means that these observables attain finite ranges of values, which is not the case in the most frequent situations in physics.
The only really interesting case is when ${\cal H}$ is finite-dimensional, for instance when dealing with the spin of particles.
In that case
$$A = \sum_{j=1}^n \lambda_j P_j$$
and, by definition,
$$f(A) := \sum_{j=1}^n f(\lambda_j) P_j$$
where $\lambda_1,\ldots, \lambda_n$ are the eigenvalues of $A$ and $P_j$ the corresponding orthogonal projectors onto the respective eigenspace.
You see that the used functions $f$ are defined on a discrete set and thus any issue on their regularity (analyticity) make no sense or is ``artificial''.
So, forgetting physical (ir)relevance, we assume that ${\cal H}$ is infinite dimensional, but $A$, $B$, and thus $[A,B]$ are everywhere defined, thus bounded -- $||A||, ||B|| < +\infty$ -- and $A$ and $B$ are Hermitian (= selfadjoint) and their spectra are bounded closed sets.
The definition of $f(A)$, and this definition is valid also if $A=A^*$ has (dense) domain $D(A) \subsetneq {\cal H}$,
is
$$f(A) := \int_{\sigma(A)} f(\lambda) dP^{(A)}(\lambda)\tag{2}$$
where $f: \mathbb{R} \to \mathbb{C}$ is Borel measurable and $P^{(A)}$ is the spectral measure of $A$. Actually, only the restriction of $f$ to $\sigma(A)$ matters here.
If $A$ and $f$ satify the said strong hypotheses
$$\left|\left|\int_{\sigma(A)} f(\lambda) dP^{(A)}(\lambda)\right|\right| \leq ||f|_{\sigma(A)}||_\infty\:,$$
I will use this inequality below.
Within this setup let us consider the various raised questions.
- $[A,B]=0$ implies $[f(A),g(B)]=0$.
YES it is true for every choice of Borel measurable (for instance continuous) $f, g : \mathbb{R} \to \mathbb{C}$ such that $f$ is bounded on the spectrum of $A$ and $g$ is bounded on the spectrum of $B$.
The proof is based on the fact that, in the considered case $[A,B]=0$ is equivalent to the commutativity of the spectral measures of $A$ and $B$ and this fact, in turn, implies the thesis.
Sketch of proof. Consider a sequence of simple functions $s_n(x) = \sum_{k=1}^{N_n} s^{(n)}_{k} \chi_{E_{n,k}}(x)$ such that
converges uniformly to the function $f$ in a compact $[a,b]$ including $\sigma(A)$ and $\Sigma(B)$.
A similar sequence $t_n= \sum_{h=1}^{N_n} t^{(n)}_{h} \chi_{F_{n,h}}(x)$ uniformly converges to $g$ on $[a,b]$.
According to (2)
$$s_n(A) = \sum_{k=1}^{N_n} s^{(n)}_{k} P^{(A)}_{E_{n,k}}$$
and
$$t_n(B) = \sum_{h=1}^{N_n} s^{(n)}_{h} P^{(A)}_{F_{n,h}}$$
where $[P^{(A)}_{E}, P^{(B)}_{F}]=0$ as a consequence (it is true in our case) of $[A,B]=0$.
We also have
$$||f(A) -s_n(A)|| \leq ||f-s_n||_\infty \to 0\:, \quad ||g(B) -t_n(B)|| \leq ||g-t_n||_\infty \to 0$$
as $n \to +\infty$.
And also
$$[f(A),g(B)] = \lim_{n\to +\infty} \left[\sum_{k=1}^{N_n} s^{(n)}_{k} P^{(A)}_{E_{n,k}}, \sum_{h=1}^{N_n} t^{(n)}_{h} P^{(B)}_{F_{n,h}}\right]= \sum_{k_1}^{N_n}\sum_{h=1}^{N_n} s^{(n)}_{k} t^{(n)}_{h} [P^{(A)}_{E_{n,k}}, P^{(B)}_{F_{n,h}}]=0$$
- $[A,B] \neq 0$ implies $[f(A),g(B)]\neq 0$.
NO. This is in general false. An easy counterexample was provided by @Oбжорoв: take $f(x) =g(x) = 1$ for every $x\in \mathbb{R}$. In that case $f(A)=g(B)=I$.
- and 4. How to prove $[A,g(B)]= g'(B)[A,B]$, if $A,B$ are Hermitian everywhere defined in ${\cal H}$ and $[[A,B],B]=0$?
Without this last hypothesis I do not know a proof and I think the identity is false (see the final ADDENDUM). I also assume that $f$ is $C^1$ (no analyticty is necessary).
Proposition.
If $A,B$ are Hermitian, everywhere defined in ${\cal H}$, $[[A,B],B]=0$, and $f$ is a $C^1$ function defined on an bounded interval $[a,b] \supset\sigma(B)$, then
$$[A,g(B)]= g'(B)[A,B]\:.$$
Proof. Using an easy extension of the Stone-Weierstrass theorem, one proves that, since $[a,b]$ is compact and $g$ is $C^1$, there is a sequence of polynomials $p_n$ such that $p_n \to g$ and $p'_n \to g'$ uniformly on $[a,b]$.
As already established in another post in PSE via direct algebraic manipulations, $$[A,p_n(B)]= p_n'(B)[A,B]\tag{1}$$ is true for polynomials when $[[A,B],B]=0$ and assuming that $A$ and $B$ are everywhere defined. We have
$$||[A,p_n(B)] - [A,g(B)]|| = ||[A,p_n(B)-g(B)]|| \leq 2||A||\: ||p_n(B)-g(B)||$$ $$ \leq 2||A||\: ||p_n(x)-g(x)||_\infty \to 0 \quad \mbox{for $n\to +\infty$.}$$
Analogously
$$||p_n'(B)[A,B]- g'(B)[A,B]|| \leq ||p_n'(B)- g'(B)|| 2||A||\:||B||$$
$$\leq ||p'_n-g'||_\infty 2||A||||B|| \to 0\:.$$
Since $[A,p_n(B)] \to [A,g(B)]$ and $p_n'(B)[A,B]\to g'(B)[A,B]$ in the operator norm as $n\to +\infty$ and (1) holds, we conclude that
$$[A,g(B)]= g'(B)[A,B]\:.$$ QED
As a final comment, I stress that one could (should?) investigate what happens when assuming less rigid hypotheses more close to physics, i.e., $A$ and $B$ unbounded, selfadjoint defined in dense domains. In this case $[A,B]=0$ does not make sense in view od problems with domains, but in spirit it means that the observables $A$ and $B$ are compatible, i.e., their spectral measures commute...
ADDENDUM. The hypothesis $[[A,B],B]=0$ cannot be easily relaxed. Indeed, suppose that $[A,g(B)]= g'(B)[A,B]$ for a suitable class of functions $g$ as the function $C^1$, that includes the polynomials. We conclude that $[A,B^2]= 2B[A,B]$. On the other hand $[A,B^2]= B[A,B] + [A,B]B$. Taking the difference, $B[A,B]-[A,B]B=0$ that is $[B,[A,B]]=0$, namely $-[[A,B],B]=0$
