Interesting representation of $e^x$ So I discovered the following formula by using the Taylor series for $\ln (x+1)$ $$x= \ln (x+1)+\frac{1}{2}\ln(x^2+1)-\frac{1}{3}\ln(x^3+1)+\frac{1}{2}\ln(x^4+1)-\frac{1}{5}\ln(x^5+1)-\frac{1}{6}\ln(x^6+1)$$$$-\frac{1}{7}\ln(x^7+1)+\frac{1}{2}\ln(x^8+1)-\frac{1}{10}\ln(x^{10}+1)-\frac{1}{11}\ln(x^{11}+1)-\frac{1}{6}\ln(x^{12}+1)$$$$-\frac{1}{13}\ln(x^{13}+1)-\frac{1}{14}\ln(x^{14}+1)+\frac{1}{15}\ln(x^{15}+1)+\frac{1}{2}\ln(x^{16}+1)+...$$
I then realized that this means that $$e^x=\prod_{n=1}^{\infty}(x^n+1)^{a_n},$$ where $a_1 =1$ and $a_n$ is defined by the recurrence relation $$a_n=\sum_{d|n \land d > 1} \frac{a_{\frac{n}{d}}(-1)^d}{d}$$ when $n > 1$.
So far, I've noticed the following properties about $a_n$:$$a_{2^m}=\frac{1}{2}$$$$a_p=-\frac{1}{p}$$$$a_{p^k}=0$$ for prime $p>2$, positive integer $m$, and integer $k>1$.
When I plugged in the denominators of $a_n$ into OEIS, the closest sequence I got was the sequence of integer radicals.

I am interested in the following things:
How can we prove my conjecture about $a_n$ (look below, under "EDIT")? For what values of $x$ does this product converge? Has this formula for $e^x$ been documented anywhere (I'm sure it has, but I'd like to read the paper)? How does this representation of $e^x$ relate to other representations?

Thanks in advance!
EDIT
After plugging in more values for $a_n$, I realized that my original conjecture that $$\lvert a_n \rvert = \begin{cases}
  \displaystyle\operatorname{rad}(n)^{-1}, & \mbox{if } n \neq p^k \\
  0,  & \mbox{if } n=p^k 
\end{cases}$$ for prime $p>2$ and integer $k>1$, was incorrect when I found that $a_{18}=0$. However, I was able to formulate a new conjecture about $a_n$ based on my findings: $$ a_n = \frac{\mu \left( \operatorname{Od}(n) \right) }{\operatorname{rad}(n)} = \frac{\mu \left( \frac{n}{2^{\nu_2 (n)}} \right) }{\operatorname{rad}(n)} $$
where $\mu(n)$ is the Möbius function, $\nu_p(n)$ is the p-adic order of $n$, $\operatorname{Od}(n)$ is the odd part of $n$, and $\operatorname{rad}(n)$ is the radical of $n$. I had a friend test values for this and it holds for all values up to at least $1024$. Unfortunately, though, I have no idea how to prove this conjecture.
NOTE
@ZhenhuaLiu has answered my biggest question by proving my conjecture. However, if you do have an answer to any of my other questions, feel free to leave an answer about it.
 A: Following Michael's suggestion,
if
$x
=\sum_{n=1}^{\infty} a_n \ln(x^n+1)
$,
then
$1
=\sum_{n=1}^{\infty} a_n \frac{nx^{n-1}}{x^n+1}
$,
so
$\begin{array}\\
1
&=\sum_{n=1}^{\infty} a_n \frac{nx^{n-1}}{x^n+1}\\
&=\sum_{n=1}^{\infty} a_n nx^{n-1}\sum_{m=0}^{\infty} (-1)^m x^{mn}\\
\text{or}\\
x
&=\sum_{n=1}^{\infty} a_n nx^{n}\sum_{m=0}^{\infty} (-1)^m x^{mn}\\
&=\sum_{n=1}^{\infty} a_n n\sum_{m=0}^{\infty} (-1)^m x^{mn+n}\\
&=\sum_{n=1}^{\infty} a_n n\sum_{m=0}^{\infty} (-1)^m x^{n(m+1)}\\
&=\sum_{n=1}^{\infty} a_n n\sum_{m=1}^{\infty} (-1)^{m-1} x^{nm}\\
&=\sum_{k=1}^{\infty}\sum_{d|k} a_d d (-1)^{k/d-1} x^{k}\\
&=\sum_{k=1}^{\infty}x^k\sum_{d|k} a_d d (-1)^{k/d-1}\\
\text{so that}\\
a_1
&=1\\
\text{and}\\
0
&=\sum_{d|k} a_d d (-1)^{k/d-1}
\qquad\text{for } k > 1\\
&=ka_k+\sum_{d|k, d<k} a_d d (-1)^{k/d-1}
\text{or}\\
k a_k
&=\sum_{d|k, d<k} a_d d (-1)^{k/d}\\
\text{or}\\
 a_k
&=\dfrac1{k}\sum_{d|k, d<k} a_d d (-1)^{k/d}\\
&=\dfrac1{k}\left((-1)^k+\sum_{d|k, 1<d<k} a_d d (-1)^{k/d}\right)\\
\end{array}
$
In particular,
if $p$ is prime,
$a_p
=\dfrac{(-1)^p}{p}
$,
so
$a_2 = \frac12$
and
$a_p = \dfrac{-1}{p}$
if $p$ is an odd prime.
If
$k = 2^m$,
$\begin{array}\\
a_{2^m}
&=\dfrac1{2^m}\left(1+\sum_{d|2^m, 1<d<2^m} a_d d (-1)^{2^m/d}\right)\\
&=\dfrac1{2^m}\left(1+\sum_{j=1}^{m-1} a_{2^j} 2^j \right)\\
\end{array}
$
Therefore
$a_4
=\frac14(1+2a_2)
=\frac12
$.
If
$a_{2^j}
=\frac12
$
for
$1 < j < m$,
then
$\begin{array}\\
a_{2^m}
&=\dfrac1{2^m}\left(1+\sum_{j=1}^{m-1} \frac12 2^j \right)\\
&=\dfrac1{2^m}\left(1+\sum_{j=0}^{m-2}  2^j \right)\\
&=\dfrac1{2^m}\left(1+2^{m-1}-1 \right)\\
&=\frac12\\
\end{array}
$
for all $m$.
If
$k = p^2$,
$a_{p^2}
=\dfrac1{p^2}(-1+pa_p(-1)^p)
=\dfrac1{p^2}(-1+p\frac{-1}{p}(-1))
=0
$.
If
$k = p^m$
where $m > 2$,
suppose
$a_{p^j}
= 0
$
for
$1 < j < m$.
This is true for
$m=3$.
Then
$\begin{array}\\
a_{p^m}
&=\dfrac1{p^m}\left(1+\sum_{d|p^m, 1<d<p^m} a_d d (-1)^{p^m/d}\right)\\
&=\dfrac1{p^m}\left(1+\sum_{j=1}^{m-1} a_{p^j} p^j \right)\\
&=\dfrac1{p^m}\left(1+a_p p +\sum_{j=2}^{m-1} a_{p^j} p^j \right)\\
&=0\\
\end{array}
$
If
$k=2p$
where $p$ is an odd prime,
$\begin{array}\\
a_{2p}
&=\dfrac1{2p}\left(1+2a_2(-1)^p+a_pp(-1)^2\right)\\
&=\dfrac1{2p}\left(1-2\dfrac12+\dfrac{-1}{p}p\right)\\
&=\dfrac1{2p}\left(1-1-1\right)\\
&=\dfrac{-1}{2p}\\
\end{array}
$
If
$k=pq$
where $p$ and $q$ are distinct odd primes,
then
$a_{pq}
=\dfrac1{pq}\left((-1)^{pq}+a_pp(-1)^p+a_qq(-1)^q\right)
=\dfrac1{pq}\left(-1+1+1\right)
=\dfrac{1}{pq}
$.
I'll leave it at this.
A: We'll use $D(f(n),s)$ to denote the Dirichlet series of arithmetic function $f(n)$.
The recurrence relation obtained by marty cohen is $$\sum_{d|n}da_{d}(-1)^{\frac{n}{d}-1}=
\begin{cases}
1 &\mbox{if } n=1,\\
0&\mbox{if } n>1.
\end{cases}
$$
We can rewrite it as Dirichlet convolution
$$
na_{n}*(-1)^{n-1}=I,
$$
where $I$ is the multiplicative identity of Dirichlet convolution.
The convolution can be translated into the language of Dirichlet series as
$$
D(na_{n},s)D((-1)^{n-1},s)=1.
$$
On the other hand, we have
$$
D((-1)^{n-1},s)=\eta(s)=(1-2^{1-s})\zeta(s),
$$
where $\eta(s)$ is the Dirichlet eta function.
Thus we have
$$\begin{align*}
D(na_n,s)&=\frac{1}{1-2^{1-s}}\frac{1}{\zeta(s)}\\
&=\left(\sum_{n\ge 0}\frac{2^n}{(2^n)^s}\right)D(\mu(n),s).\\
\end{align*}.$$
Translating the product of Dirichlet series into the language of Dirichlet convolution, we have
$$\begin{align*}
a_n&=\frac{1}{n}\sum_{d|n}d[\log_{2}d\in \mathbb{N}]\mu(\frac{n}{d})\\
&=\frac{1}{n}\sum_{k=0}^{\nu_{2}(n)}2^k\mu(\frac{n}{2^k}),\\
\end{align*}
$$
where $[P]$ denotes the Iverson bracket.
If $\mu(\frac{n}{2^{\nu_{2}(n)}})=0$, clearly we have $$a_n=0=\frac{\mu(\frac{n}{2^{\nu_{2}(n)}})}{\mathrm{rad}(n)}.$$
If $\mu(\frac{n}{2^{\nu_{2}(n)}})\neq 0$, we have 
$$
a_n=\frac{1}{n}\mu(n)=\frac{\mu(\frac{n}{2^{\nu_{2}(n)}})}{\mathrm{rad}(n)},$$
when $\nu_{2}(n)=0,$
and
$$\begin{align*}
a_n&=\frac{1}{n}\left(2^{\nu_{2}(n)}\mu(\frac{n}{2^{\nu_{2}(n)}})+2^{\nu_{2}(n)-1}\mu(2\frac{n}{2^{\nu_{2}(n)}})\right)\\
&=\frac{1}{n}\left(2^{\nu_{2}(n)}\mu(\frac{n}{2^{\nu_{2}(n)}})+\frac{1}{2}2^{\nu_{2}(n)}\mu(\frac{n}{2^{\nu_{2}(n)}})\mu(2)\right)\\
&=\frac{2^{\nu_{2}(n)}}{2n}\mu(\frac{n}{2^{\nu_{2}(n)}})\\
&=\frac{\mu(\frac{n}{2^{\nu_{2}(n)}})}{\mathrm{rad}(n)},\\
\end{align*}
$$
when $\nu_{2}(n)\ge 1.$
Combining the results from all cases, we've proved your conjecture $$a_n=\frac{\mu(\frac{n}{2^{\nu_{2}(n)}})}{\mathrm{rad}(n)}.$$
A: This is more a comment
than an answer,
but it is too complex
to enter as a comment.
It turns out that
you do not need to differentiate
to get the recurrence
for the $a_n$.
Here's why.
$\begin{array}\\
x
&=\sum_{n=1}^{\infty} a_n \ln(x^n+1)\\
&=\sum_{n=1}^{\infty} a_n \sum_{m=1}^{\infty}\dfrac{(-1)^{m-1}x^{mn}}{m}\\
&=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty} a_n \dfrac{(-1)^{m-1}x^{mn}}{m}\\
&=\sum_{k=1}^{\infty}\sum_{d|k} a_d \dfrac{(-1)^{k/d-1}x^{k}}{k/d}
\qquad k = mn, d = n, m = k/d\\
&=\sum_{k=1}^{\infty}\dfrac{x^{k}}{k}\sum_{d|k} da_d (-1)^{k/d-1}\\
\text{so}\\
a_1
&=1\\
\text{and}\\
0
&=\sum_{d|k} da_d (-1)^{k/d-1}\\
\text{as before}\\
\end{array}
$
