Generally speaking functions of operators are defined as spectral functions, making use of the spectral theorem machinery.
The approaches based of power series usually have not a rigorous basis: Even identities like $$e^A = \sum_n \frac{1}{n!}A^n$$ are generally false if $A$ is an unbounded operator.
Nevertheless in most relevant physical cases some of the results obtained by formal (non-rigorous) manipulations can be obtained following alternative rigorous ways.
However the case you are considering is very easy. Suppose that $F$ is a smooth, generally non-analytic, complex-valued function bounded by some polynomial of $\vec{x}$.
In this case both $F(\vec{x})$ and $p_i$ admit the space of Schwartz functions ${\cal S}(\mathbb R^3)$ as invariant space, so that $[p_i, F(\vec{x})]\psi$
is well defined for $\psi \in {\cal S}(\mathbb R^3)$. On that space $p_i$
coincides to $-i \hbar\frac{\partial}{\partial x_i}$. By direct computation:
$$[p_i, F(\vec{x})]\psi = -i \hbar\frac{\partial}{\partial x_i} F(\vec{x}) \psi(x) + i \hbar F(\vec{x})\frac{\partial \psi}{\partial x_i}= -i\hbar
\frac{\partial F}{\partial x_i} \psi\:.$$
We have obtained that
$$\left([p_i, F(\vec{x})] + i\hbar
\frac{\partial F}{\partial x_i} \right) \psi =0 \quad \forall \psi \in {\cal S}(\mathbb R^3)\:.$$
In other words, at least on the domain ${\cal S}(\mathbb R^3)$:
$$[p_i, F(\vec{x})] = -i\hbar
\frac{\partial F}{\partial x_i} $$
In general, this identity holds on larger domains and, in fact, it can be extended exploiting some further know property of $F$ and some other properties like the fact that $\cal S(\mathbb R^3)$ is a core for $p_i$...
One could assume weaker hypotheses. The self-adjointness domain of $p_j$ is so made.
$$D(p_j) := \{\psi \in L^2(\mathbb R^3)\:|\: \exists\: w\mbox{-}\partial_{x_j}\psi \in L^2(\mathbb R^3)\}$$
where $w\mbox{-}\partial_{x_j}\psi$ denotes the weak partial $j$-derivative of $\psi$. On that domain $p_j = w\mbox{-}\partial_{x_j}$ as expected.
If $F$ is, for instance, just $C^1$ and compactly supported, all above reasoning can be re-implemented with $\psi \in D(p_j)$, obtaining exactly the same result, using the fact that the weak derivative of a compactly supported $C^1$ function $F$ coincides with the standard one.
