**Edit:** Made changes to the code For each $n\in\mathbb{N}$, how do we compute sets $A_n$ and $B_n$ below: > Let $A_1=[0,2/3)$. Let $B_1=(2/3,1]$. > > If $A_n$ is a union of intervals, then for each interval cut out the > middle $1/2^{n+1}$ of the interval and send it to $B_n$. Similarly, > for each interval in $B_n$ cut out the middle $1/(2^n)$ and send it to > $A_n$. Each set less what’s cut out plus what was transferred > determines $A_{n+1}$,$B_{n+1}$ . Here is my attempt. We start with: Clear["Global`*"] A[1] = 0 <= x < 2/3 (*A1=[0,2/3)*) B[1] = 2/3 <= x < 1 (*B1=[2/3,1)*) Taking the upper and lower bound of intervals $A_1$ and $B_1$ a[1, 1] = A[1][[1]] (*Lowerbound of A1*) a[1, 2] = A[1][[5]] (*Upperbound of A1*) b[1, 1] = B[1][[1]] (*Lowerbound of B1*) b[1, 2] = B[1][[5]] (*Upperbound of B1*) Defining $A_2$ and $B_2$ A[2] = Reduce[a[1, 1] <= x <= (a[1, 1] + a[1, 2])/2 - 1/2 (1/2^(2 + 1)) (a[1, 2] - a[1, 1]) || (a[1, 1] + a[1, 2])/2 + 1/2 (1/2^(2 + 1)) (a[1, 2] - a[1, 1]) <= x < a[1, 2] || (b[1, 1] + b[1, 2])/2 - 1/2 (1/2^(2)) (b[1, 2] - b[1, 1]) <= x < (b[1, 1] + b[1, 2])/2+ 1/2 (1/2^(2)) (b[1, 2] - b[1, 1]), {x}] (*A2*) B[2] = Reduce[b[1, 1] <= x <= (b[1, 1] + b[1, 2])/2 - 1/2 (1/2^(2)) (b[1, 2] - b[1, 1]) || (b[1, 1] + b[1, 2])/2 + 1/2 (1/2^(2)) (b[1, 2] - b[1, 1]) <= x < b[1, 2] || (a[1, 1] + a[1, 2])/2 - 1/2 (1/2^(2 + 1)) (a[1, 2] - a[1, 1]) <= x < (a[1, 1] + a[1, 2])/2 + 1/2 (1/2^(2 + 1)) (a[1, 2] - a[1, 1]), {x}] (*B2*) Repeating the same process for the each of the $3$ intervals of $A_2$, each of the $9$ intervals of $A_3$, and each of the $3^{n-1}$ intervals of $A_n$. a[x_, 1, i_] := A[x][[3^(i - 1)]][[1]] (*If i is a number in [1,3^(x-1)], then a[x,1,i] is the Lowerbound of the 3^(i-1)th interval of A_n*) a[x_, 2, i_] := A[x][[3^(i - 1)]][[5]] (*If i is a number in [1,3^(x-1)], then a[x,1,i] is the Upperbound of the 3^(i-1)th interval of A_n*) b[x_, 1, i_] := B[x][[3^(i - 1)]][[1]] (*If i is a number in [1,3^(x-1)], then b[x,1,i] is the Lowerbound of the 3^(i-1)th interval of B_n*) b[x_, 2, i_] := B[x][[3^(i - 1)]][[5]] (*If i is a number in [1,3^(x-1)], then b[x,1,i] is the Upperbound of the 3^(i-1)th interval of B_n*) A[s_, i_] := Reduce[a[s, 1, i] <= x <= (a[s, 1, i] + a[s, 2, i])/2 - 1/2 (1/2^(2 + 1)) (a[s, 2, i] - a[s, 1, i]) || (a[s, 1, i] + a[s, 2, i])/2 + 1/2 (1/2^(2 + 1)) (a[s, 2, i] - a[s, 1, i]) <=x < a[s, 2, i] || (b[s, 1, i] + b[s, 2, i])/2 - 1/2 (1/2^(2)) (b[s, 2, i] - b[s, 1, i]) <= x < (b[s, 1, i] + b[s, 2, i])/2 + 1/2 (1/2^(2)) (b[s, 2, i] - b[s, 1, i]), {x}] (*Changes the 3^(i-1)th interval of A_n using the description in above the attempt.*) B[s_, i_] := Reduce[b[s, 1, i] <= x <= (b[s, 1, i] + b[s, 2, i])/2 - 1/2 (1/2^(2)) (b[s, 2, i] - b[s, 1, i]) || (b[s, 1, i] + b[s, 2, i])/2 + 1/2 (1/2^(2)) (b[s, 2, i] - b[s, 1, i]) <= x <b[s, 2, i] || (a[s, 1, i] + a[s, 2, i])/2 - 1/2 (1/2^(2 + 1)) (a[s, 2, i] - a[s, 1, i]) <= x < (a[s, 1, i] + a[s, 2, i])/2 + 1/2 (1/2^(2 + 1)) (a[s, 2, i] - a[s, 1, i]), {x}] (*Changes 3^(i-1)th interval of A_n into three intervals using the description above the attempt.*) Combining the $3^{n-1}$ changed intervals into a total of $3^n$ new intervals A[s_] := Reduce[Flatten[Table[A[s, i], {i, 1, 3^(s - 1)}]], {s}] (*Combines all 3^(i-1)th intervals for all i-values*) B[s_] := Reduce[Flatten[Table[B[s, i], {i, 1, 3^(s - 1)}]], {s}] (*Combines all 3^(i-1)th intervals for all i-values*) A[3] (*A3*) B[3](*B3*) The problem is the code works for `A[2]` and `B[2]` 0 <= x <= 7/24 || 3/8 <= x < 2/3 || 19/24 <= x < 7/8 7/24 <= x < 3/8 || 2/3 <= x <= 19/24 || 7/8 <= x < 1 However, the rest of the code returns `$Recursionlimit: Recursion depth of 1024 exceeded...`. Even if we fix the code, there's a more efficient way of writing the programming. **Question:** How should we write the code? Is there a more efficient method?