Infinite discounted sum of betas Let $0 \leq \gamma < 1$, $X_i \sim \text{Beta}(\alpha, \beta)$, and
$$Y \sim \sum_{i = 0}^\infty \gamma^i X_i$$
What is the distribution of $Y$? Does it have a closed form? Can it be sampled efficiently? We know that
\begin{align}
\operatorname{supp}(Y) &= \left[0, \frac{1}{1-\gamma}\right) \\
\operatorname{E}(Y) &= \frac{\operatorname{E}(X)}{1-\gamma} \\
\operatorname{Var}(Y) &= \frac{\operatorname{Var}(X)}{1-\gamma^2} \\
\end{align}
Also, note that
\begin{align}
Y = X + \gamma Y'
\end{align}
where $Y$ and $Y'$ have the same distribution. Therefore, denoting the characteristic function of a random variable $A$ as $\phi_A$,
\begin{align}
\phi_Y(t)
&= \phi_{X + \gamma Y'}(t) \\
&= (\phi_X \cdot \phi_{\gamma Y'})(t) \\
&= \phi_X(t) \cdot \phi_{\gamma Y'}(t) \\
&= \phi_X(t) \cdot \phi_{Y'}(\gamma t) \\
&= \phi_X(t) \cdot \phi_Y(\gamma t) \\
\frac{\phi_Y(t)}{\phi_Y(\gamma t)} &= \phi_X(t) \\
&= {}_1F_1(\alpha; \alpha + \beta; i t)
\end{align}
where the last RHS is the characteristic function of the beta distribution. Is there a way to solve for $\phi_Y(t)$ or at least sample efficiently from $Y$'s distribution?
 A: Note
I've used slightly different notation to you.
Firstly, let's consider the underlying random variables we are dealing with. We have that all
$$X_{i}\sim\text{Beta}(\alpha,\beta)$$
where all random variables are independent of each other. If we now consider some scaled form of these random variables:
$$Y=a X$$
where we have that $a<1$ because $\gamma<1$. Now, for a continuous random variable
$$\begin{align}
f_{Y}(y)&=\frac{f_{X}(y/a)}{a}\\
&=\frac{(y/a)^{\alpha-1}(1-y/a)^{\beta-1}}{a\cdot\text{B}(\alpha,\beta)}
\end{align}$$
Notice that the support has changed from
$$0<x<1$$
to
$$0<y<a$$
If we define
$$\text{B}'(\alpha,\beta,a)=\int_{0}^{a}t^{\alpha-1}(a-t)^{\beta-1}\,dt$$
as a variant on the usual Beta function (this function equals the Beta function for $a=1$), then we can show that
$$\text{B}'(\alpha,\beta,a)=a^{\alpha+\beta-1}\cdot\text{B}(\alpha,\beta)$$
This means that
$$\begin{align}
f_{Y}(y)&=\frac{(y/a)^{\alpha-1}(1-y/a)^{\beta-1}}{a\cdot\text{B}(\alpha,\beta)}\\
&=\frac{y^{\alpha-1}(a-y)^{\beta-1}}{\text{B}'(\alpha,\beta,a)}\cdot\frac{\text{B}'(\alpha,\beta,a)}{a^{\alpha+\beta-1}\cdot\text{B}(\alpha,\beta)}\\
&=\frac{y^{\alpha-1}(a-y)^{\beta-1}}{\text{B}'(\alpha,\beta,a)}
\end{align}$$
For the purposes of this answer, let's denote this random variable as
$$Y\sim\text{Beta}'(\alpha,\beta,a)$$
This is essentially a scaled Beta distribution, where the support has been limited to $[0,a]$. The important point here is that random values generated from this random variable will be bound between $[0,a]$.
Now, in your summation, each Beta distribution $X_{i}$ is scaled by $\gamma^{i}$ for increasing $i$. This means that each summand will have an increasingly reduced support $[0,\gamma^{i}]$ and each subsequent summand's probability mass will be increasingly narrow.
We know this series will converge for $0\leq\gamma<1$ because, for each $i$, the maximum value that each subsequent random value can take is $\gamma^{i}$. So, at most, the series is
$$1+\gamma+\gamma^{2}+\gamma^{3}+\gamma^{4}+\ldots$$
which is a simple geometric series converging to
$$\frac{1}{1-\gamma}$$
Interestingly, this tells us the upper bound of the support of $Y$ for a given $\gamma$. If we define
$$Z=\sum_{i=0}^{\infty}Y_{i}$$
then the support for $Z$ would then be $[0,1/(1-\gamma)]$.
Now, you can derive the moment generating function of $Y$ as
$$M_{Y_{i}}(t)=1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{(\gamma^{i}\cdot t)^{k}}{k!}\Bigg)$$
Then the moment generating function of $Z$ is
$$\begin{align}
M_{Z}(t)&=\prod_{i=0}^{\infty}M_{Y_{i}}(t)\\
&=\Bigg[1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{t^{k}}{k!}\Bigg)\Bigg]\cdot \\
&\Bigg[1+\sum_{k=1}^{\infty}\Bigg(\prod_{r=0}^{k-1}\bigg(\frac{\alpha+r}{\alpha+\beta+r}\bigg)\frac{(\gamma\cdot t)^{k}}{k!}\Bigg)\Bigg]\cdots\\
\end{align}$$
I'm not sure how useful this is. But it's useful to note that for increasing $i$ the subsequent terms will decay in size so that $M_{Y_{i}}(t)\rightarrow 1$ for large $i$.
If we look at some simulations of this distribution, there are some interesting results. First we'll define some functions for $Y$.
    dsbeta=function(x, alpha, beta, gamma) {
      x[x<0]=0
      x[x>gamma]=gamma
      return(x^(alpha-1)*(gamma-x)^(beta-1)/(gamma^(
                alpha+beta-1)*beta(alpha, beta)))
    }
    
    psbeta=function(q, alpha, beta, gamma) {
      q[q<0]=0
      q[q>gamma]=gamma
      return(pbeta(q/gamma, alpha, beta))
    }
    
    qsbeta=function(p, alpha, beta, gamma) {
      p[p<0]=0
      p[p>1]=1
      return(gamma*qbeta(p, alpha, beta))
    }
    
    rsbeta=function(n, alpha, beta, gamma) {
      return(gamma*rbeta(n, alpha, beta))
    }

    #Perform some simulations:
    a=3
    b=2
    gamma=0.2 #Scaling of each X
    N=10000   #Number of samples of Z
    M=100     #Number of Y in each Z
    
    #Generate random values of Z:
    r=as.vector(matrix(rbeta(N*M, a, b), N, M) %*%  
                       gamma^seq(0,M-1))
    
    #Final simulations:
    hist(r, freq = FALSE, breaks = 50)
    
    #Perform fit (here we provide the upper bound gamma as the 
    #             theoretical upper limit):
    f=fitdistrplus::fitdist(r, "sbeta", start = list(alpha = 2, 
         beta = 3), fix.arg = list(gamma = 1/(1-gamma)))
    
    #Visualize the fit:
    x=seq(min(r), max(r), length.out = 1000)
    lines(x, dsbeta(x, f$estimate[1], f$estimate[2], 
             1/(1-gamma)), lwd = 2, col = "blue")

    #GOF tests:
    ks.test(r, "psbeta", f$estimate[1], f$estimate[2],  
               f$fix.arg$gamma)
    goftest::cvm.test(r, "psbeta", f$estimate[1], f$estimate[2], 
                      f$fix.arg$gamma)$p.value
goftest::ad.test(r, "psbeta", f$estimate[1], f$estimate[2], 
                    f$fix.arg$gamma)$p.value
    
    #qq plot:
    plot(sort(r), qsbeta(ppoints(length(r)), f$estimate[1], 
     f$estimate[2], f$fix.arg$gamma), ylab = "Theoretical 
         quantiles", xlab = "Empirical quantiles", col = "blue", 
                     type = "l", lwd = 2)
    abline(0, 1, lty = 2, lwd = 2, col = "red")

It seems from the simulations that $Z$ can be well approximated by the scaled beta
$$Z\sim\text{Beta}'(\cdot,\cdot,1/(1-\gamma))$$
but as per @whuber's comments, for small $\alpha,\beta$ and $\gamma\rightarrow 0$ this breaks down. I think it has something to do with the bi-modality for $\alpha,\beta<1$ and the fewer 'non-degenerate' random variables due to $\gamma$ small.
