I'm now trying to derive the Taylor series for differentiable functions. For simplicity, let $f \in C^2(\mathbb R)$. If $f(a+h)$ is expanded to a polynomial of $h$ with order 1, that is,
$$f(a+h) = f(a) + \alpha h + o(1)$$
for some $\alpha \in \mathbb R$, by the definition of derivative, we have $\alpha = f'(a)$.
Next, we assume $f(a+h)$ is expanded to a polynomial of $h$ with order 2 in the form of
$$f(a+h) = f(a) + f'(a) h + \beta h^2 + o(h^2).$$
If such $\beta$ exists, it must be
$$\beta = \lim_{h \to 0} \frac{f(a+h) - f(a) - f'(a)h}{h^2}$$
but, I don't know why the factor $1/2!$ appears. In general, does this limit exist and coincide with $f''(a)/2!$ if $f$ is smooth enough?
As far as I am concerned, the factor $\dfrac{1}{2!}$ is not that hard to see where it comes from. Also, Taylor series as a whole does contains a lot more, than only one fractional factorial.
Now, let's recall the Taylor series, by definition:
Now, should we try to decode, or derive the series, we will want to make some assumptions:
Then, given an arbitrary function $f$ that fits the assumption we have, then, for us to gain a detailed power series representation of it, we need the coefficient. How do we find it? Let's take a look at the power series we have. To find $x_0$, set $x=a$, then $$f\left( a \right) = \sum\limits_{n = 0}^\infty {{c_n}{{\left( {a - a} \right)}^n}} = {c_0} + {c_1}\left( {a - a} \right) + {c_2}{\left( {a - a} \right)^2} + {c_3}{\left( {a - a} \right)^3} + {c_4}{\left( {a - a} \right)^4} + \cdots$$ which gives us $$f(a) = c_0$$ Next, if we differentiate the power series, we gain $$f'(x) = \dfrac{d}{dx}\left( \sum_{n=0}^{\infty}c_n(x−a)^n \right)=c_1+2c_2(x−a)+3c_3(x−a)^2+\dots.$$ which then again, set $x=a$, we have $$f'(a) = \dfrac{d}{dx}\left( \sum_{n=0}^{\infty}c_n(a−a)^n \right)=c_1+2c_2(a−a)+3c_3(a−a)^2+\dots = c_1$$
Up to this point, you might see the patterns. To find an incremental amount of coefficient of $f$, which itself is $f\in C^{n}(\mathbb{R})$, means you will want to differentiate the power series as a whole to an order of derivative, then set $x=a$, thus gaining the coefficient of that order of expansion. Continuing, we have $$f''\left( x \right) = \dfrac{d^2}{dx^2} \left(\sum_{n=0}^{\infty}c_n(x−a)^n \right)=2c_2+3⋅2c_3(x−a)+4⋅3c_4(x−a)^2+\dots$$ which then, $$f''(a) = 2c_2$$ We only want $c_2$ so we divide it by $2$, gaining $c_2$. Thus $$\dfrac{f''(a)}{2}=c_2$$ For $f'''(x)$: $$f'''(x) = \dfrac{d^3}{dx^3} \left(\sum_{n=0}^{\infty}c_n(x−a)^n\right)=3⋅2c_3+4⋅3⋅2c_4(a−a)+5⋅4⋅3c_5(a−a)^2+\dots =3⋅2c_3$$ Then the coefficient $c_3$ for $x= a$is again $$c_3 = \dfrac{f'''(a)}{3\cdot 2}= \dfrac{f'''(a)}{3!}$$
Do this for $f^{(n)}(x)$ and you will eventually gain the result: $$f(x) = f\left( a \right) + f'\left( a \right)\left( {x - a} \right) + \frac{{f''\left( a \right)}}{{2!}}{\left( {x - a} \right)^2} + \frac{{f'''\left( a \right)}}{{3!}}{\left( {x - a} \right)^3} + \cdots$$ which is shorten to $$f(x)=\sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!}(x-a)^n$$
So as you can see, the factor $\dfrac{1}{2!}$, as well as others, is from differentiation of the power series of the function $f$.