This might be complicated to describe rigorously, using elementary notions, but the underlying concept is simple: almost all real numbers between $0$ and $1$, when written in binary, have equally many zeros and ones. Because $Z$ assigns zero probability to such numbers--it favors an imbalance of zeros and ones when $p\ne 1-p$--it must have a singular distribution.
The statement in bold needs a clear definition. I am going to sneak up on the definition by developing a series of simpler definitions and notations. But before we proceed, let's look at some of these distributions.
The colored curves plot the CDFs. They are all self-similar in the sense suggested by the underlying rectangles: The curve in the lower left rectangle is obtained by taking the entire curve and multiplying all x-values by $1/2$ and multiplying all y=values by $1-p$. That is, the curve is shrunk towards the origin $(0,0)$ by factors of $1/2$ and $1-p$. Likewise, the curve in the upper right rectangle is obtained by shrinking the entire curve towards the coordinate $(1,1)$ by factors of $1$ and $p$. (This description completely determines each curve!)
Here are the major steps in the analysis.
Let $a_1, a_2, \ldots, a_n$ be any sequence of zeros and ones. It represents a binary number in the interval $[0,1)$ by "putting a binary point in front of it." As a matter of notation, and to define precisely what this means, let $$[a_1\,a_2\,\cdots\,a_n] = a_1 2^{-1} + a_2 2^{-2} + \cdots + a_n 2^{-n}.$$
Since we will be studying the proportions of ones within binary representations, let's invent a notation for that. Define $$\rho([a_1\,a_2\,\cdots\,a_n]) = \frac{\#\{i \in \{1,2,\ldots, n\}\mid a_i=1\}}{n}$$ to be that proportion.
We cannot analyze all the infinitely many binary digits of any real number $x$ all at once: we have to look at finitely many of them and take limits. To this end, extend the function $\rho$ to entire intervals of real numbers determined by their initial sequences. That is, if $$[a_1\,a_2\,\cdots\,a_n] \le x \lt [a_1\,a_2\,\cdots\,a_n] + 2^{-n},$$ then define $$\rho_n(x) = \rho([a_1\,a_2\,\cdots\,a_n]).$$ The functions $\rho_n$ are defined on the entire interval $[0,1)$.
If you like, you could also define the $\rho_n(x)$ numerically in terms of the binary expansion of the integer part of $2^n x$ (which is how the next figure was produced), but I thought the preceding description might give a clearer idea of the underlying concept.
The key idea is to collect all numbers that look, based on the initial (finite!) part of their binary representation, like they have a given proportion of ones. We should allow for that proportion to be in a narrow range, say between $\alpha$ and $\beta$ where $0 \le \alpha \le \beta \le 1$. To this end, define the (Borel measurable) set $R$ by $$R(n,\alpha,\beta) = \{x\in [0,1)\mid \alpha \le \rho_n(x) \le \beta\}.$$ That is, the proportion of ones within the first $n$ binary digits of $x$ lies between $\alpha$ and $\beta$.
In the preceding figure, you may visualize $R(n,\alpha,\beta)$ by referring to the graph of $\rho_n$, drawing the interval $[\alpha,\beta]$ on the vertical axis, identifying all parts of the graph whose heights are within that interval, and looking down at the horizontal ($x$) axis to see which numbers correspond to those heights: that (highly disconnected) set is $R(n,\alpha,\beta)$.
At this point we may leave the world of binary representations and focus on the digits themselves. When those first $n$ binary digits are generated by $X_1, X_2, \ldots, X_n$, each having a chance $p$ of being $1$, then the proportion of ones has a Binomial$(n,p)$ distribution. Thus, unless $X_{n+1}=X_{n+2}=\cdots = 1$, $$\rho_n(Z) = \rho[X_1\,X_2\,\cdots\,X_n] \sim \frac{1}{n}\text{Binomial}(n,p).$$ We may exclude all values of $Z$ that end in all ones because (as is easy to see and show) their chance of occurring is zero.
Consequently, writing $F_{n,p}$ for the Binomial$(n,p)$ distribution function, $$\mathbb{P}(Z \in R(n,\alpha,\beta)) = F_{n,p}(n\beta) - F_{n,p}(n\alpha).$$
Now suppose $p\ne 1/2$. Pick $\alpha,\beta$ such that $p\in[\alpha,\beta]$ but $1/2\notin[\alpha,\beta]$. It is easy to prove (in many ways, ranging from the Weak Law of Large Numbers to Chebyshev's Inequality) that as $n$ grows large, the value of $F_{n,p}(n\beta)-F_{n,p}(n\alpha)$ approaches $1$: that is, the Binomial distribution becomes concentrated around its expected value $p$. This means as $n$ grows large, all the probability of $Z$ becomes concentrated within the (very complicated union of intervals) $R(n,\alpha,\beta)$.
However, the result for the case $p=1/2$ (for which the distribution of $Z$ is Lebesgue measure) shows that the Lebesgue measure of $R(n,\alpha,\beta)$ approaches zero. Together with (7) we conclude that for $p\ne 1/2$, the distribution of $Z$ is singular (with respect to Lebesgue measure).
In fact, the same argument shows that the various distributions of $Z$ for any two different $p$ are all mutually singular: as $n$ grows large, the probability of one concentrates within a region that has vanishingly small probability for the other.
Example R code to plot one of these distribution functions follows. Function f generates the sequence of integers $n\rho_n(0), n\rho_n(2^{-n}), n\rho_n((2)2^{-n}), \ldots, n\rho_n((2^n-1)2^{-n}).$ These count the number of $X_i$ for $i=1,2,\ldots, n$ that are equal to $1$. It then computes the associated probability and accumulates it to approximate the CDF of $Z$ across an equally-spaced set of $2^n$ points in the interval $[0,1]$. This process emulates the self-similarity transformation described at the beginning of this answer.
#
# The distribution function, approximated to an accuracy of 2^(-n) in Z.
#
f <- function(n, p=1/2) {
y <- rep(0, 2^n)
if (n > 0)
for (i in 1:n) {
j <- 1:(2^(i-1)); y[j + 2^(i-1)] <- y[j] + 1
}
cumsum(p^y * (1-p)^(n-y))
}
#
# Example: plot a distribution for given `p`.
#
n <- 11
p <- 1/6
z <- c(0:(2^n) / 2^n)
Probability <- c(0, f(n, p))
plot(z, Probability, type="S")
