Final update on 11/29/2019: I have worked on this a bit more, and wrote an article summarizing all the main findings. You can read it here.
I am looking for a reference regarding the following. Every real number $x\in [1, 2]$ is uniquely representable as a product of distinct factors of the form $1+2^{-k}$, that is
$$x=\prod_{k=1}^\infty \Big(1+\frac{b_k}{2^k} \Big) \mbox{ with } b_k\in \{0, 1\}.$$
The algorithm to find the $b_k$'s is a simple version of the greedy algorithm. I could not find any reference to this fact, except a small note in the Wikipedia entry for logarithms. It is said to be related to Feynman's algorithm, see here. If all $b_k$'s are equal to one, then $x=2.384231029… = \frac{1}{1-z}$ where $z$ is the Pell constant (see here, here and here.)
The reason for my interest is because I am studying numeration systems (see my numerous questions on Math.StackExchange and Stats.StackExchange, most recently this one), and in particular, numeration systems where the digits are random numbers. In this case, it means I am interested in the case where the $b_k$'s are independent Bernouilli random variables of parameter $p_k$ (for instance $p_k=\frac{1}{2}$ or $p_k = 2^{-k}$.)
Update
It turns out that this system is ill-conditioned. A number can actually have two different representations, for instance
$$ \frac{3}{2} = 1 + \frac{1}{2^1} = \Big(1+\frac{1}{2^2}\Big)\Big(1+\frac{1}{2^3}\Big)\Big(1+\frac{1}{2^4}\Big)\Big(1+\frac{1}{2^8}\Big)\Big(1+\frac{1}{2^{16}}\Big)\Big(1+\frac{1}{2^{32}}\Big)\cdots$$
More precisely,
$$\frac{3}{2} = 1 + \frac{1}{2^1} =\Big(1+\frac{1}{2^3}\Big).\prod_{k=1}^\infty\Big(1+\frac{1}{2^{2^k}} \Big).$$
If you use the greedy algorithm to find the $b_k$'s, it will lead to the first (leftmost) representation: $\frac{3}{2}= 1+\frac{1}{2^1}$. Yet if $b_k$ has a Bernouilli distribution of parameter $\frac{1}{2}$ (that is 50% of zeros, 50% of ones), then $x$ also has some statistical distribution, pictured below (the chart below shows the percentile distribution):
The distribution (CDF) in question is very well approximated by $F(x) = \log_{\lambda} x$, with $x \in [1, \lambda]$. Here $\lambda =2.384231029…$ is the constant discussed earlier. See below the chart picturing $F(x)$ in blue, and its logarithmic approximation in red. The CDF is the inverse of the percentile distribution.
Now this is becoming very interesting: compare this chart with the one obtained for continued radicals in section 2.2 in this article. How similar! I plan on sharing my spreadsheet with the computations required to produce the percentile and PDF distribution.
The chart below shows the approximation error $\log_\lambda x – F(x)$ between the logarithmic approximation and the exact CDF, with $x\in [1, \lambda]$.
If $b_k$ has a Bernoulli distribution of parameter $p=\frac{1}{6}$ (rather than $p=\frac{1}{2}$), then $x$ has a more chaotic distribution, see percentile distribution below:
To the contrary, if $b_k$ is uniform on $\{0, 1, 2, 3\}$ then the distribution is much smoother, see percentile distribution below:
Best Answer
Not an answer but too long for a comment.
$$\log x=\sum_{k=1}^\infty\log\left(1+\frac{b_k}{2^k}\right)\le\sum_{k=1}^\infty\frac{b_k}{2^k}$$ is able to represent all numbers in $[0,\log x_1]$ (where $x_1$ is your constant). The gap between $\log\left(1+\dfrac{b_k}{2^k}\right)$ and $\dfrac{b_k}{2^k}$ quickly decreases as $k$ gets larger, and the two numeration systems become virtually identical.
The successive "weights" of the bits are
$$0.40546510\cdots\\0.22314355\cdots\\0.11778303\cdots\\0.06062462\cdots\\0.03077165\cdots\\0.01550418\cdots\\\cdots$$
to be compared to
$$0.5\\0.25\\0.125\\0.0625\\0.03125\\0.015625\\\cdots$$
The non-uniqueness is explained by the convexity of the function, implying
$$\log\left(1+\frac1{2^k}\right)<2\log\left(1+\frac1{2^{k+1}}\right)$$ so that a term can be reconstructed by higher order terms.
To obtain approximations of the $\text{cdf}$ of the random variable, using $n$ bits, you associate to the bits of all integers $k$ in $[0,2^n-1]$ the sum of the corresponding weights; these sums form the abscissas and the ordinates are $\dfrac k{2^n}$, but there are inversions.
The higher $n$, the more the portion of the curve between two points tends to a line segment.