Most efficient way of defining the following sets for every step $n$

Question

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}$ .

(Edit 2: Note I define $A_n$ and $B_n$ as half-open intervals since I want to partition the reals into sets $A$ and $B$ that are dense (with positive measure) in every sub-interval of $(a,b)$ of $\mathbb{R}$. For more info see this and this post. Thus, it wouldn't be mathematically appropriate if we have that all intervals in $A_n$ and $B_n$ are closed.)

---

Here is my attempt. We start with:

(Edit 1: Below, I made changes to the code.)

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$

(Edit 3: I made $1/2^{2+1}$ and $1/2^{2}$, $1/2^{1+1}$ and $1/2^{1}$. I apologize for my stupidity.)

A[2] = Reduce[a[1, 1] <= x <= (a[1, 1] + a[1, 2])/2 - 1/2 (1/2^(1 + 1)) (a[1, 2] - a[1, 1]) || 
       (a[1, 1] + a[1, 2])/2 + 1/2 (1/2^(1 + 1)) (a[1, 2] - a[1, 1]) <= x < a[1, 2] || 
       (b[1, 1] + b[1, 2])/2 - 1/2 (1/2^(1)) (b[1, 2] - b[1, 1]) <= x < (b[1, 1] + b[1, 2])/2+
       1/2 (1/2^(1)) (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^(1)) (b[1, 2] - b[1, 1]) || 
       (b[1, 1] + b[1, 2])/2 + 1/2 (1/2^(1)) (b[1, 2] - b[1, 1]) <= x < b[1, 2] ||
       (a[1, 1] + a[1, 2])/2 - 1/2 (1/2^(1 + 1)) (a[1, 2] - a[1, 1]) <= x < (a[1, 1] + a[1, 2])/2 + 
     1/2 (1/2^(1 + 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^(s)) 
             (a[s, 2, i] - a[s, 1, i]) || (a[s, 1, i] + a[s, 2, i])/2 + 1/2 (1/2^(s)) 
             (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^(s-1)) (b[s, 2, i] - b[s, 1, i]) <= x < (b[s, 1, i] + b[s, 2, i])/2 
             + 1/2 (1/2^(s-1)) (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^(s-1)) 
              (b[s, 2, i] - b[s, 1, i]) || (b[s, 1, i] + b[s, 2, i])/2 + 1/2 (1/2^(s-1)) 
              (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^(s)) (a[s, 2, i] - a[s, 1, i]) <= x < (a[s, 1, i] + a[s, 2, i])/2
              + 1/2 (1/2^(s)) (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?

Answer 1

Here is code that generates the endpoints of the intervals

Clear[a, b]
a[1] = {{0, 2/3}}; b[1] = {{2/3, 1}};
a[n_] := a[n] = Join[Flatten[{{#1, #1 + (#2 - #1)/2 - (#2 - #1)/2^(n + 2)}, {#1 + (#2 - #1)/2 + (#2 - #1)/2^(n + 2), #2}} & @@@a[n - 1], 1], {#1 + (#2 - #1)/2 - (#2 - #1)/2^(n + 1), #1 + (#2 - #1)/2 + (#2 - #1)/2^(n + 1)} & @@@ b[n - 1]] // Sort
b[n_] := b[n] = Join[Flatten[{{#1, #1 + (#2 - #1)/2 - (#2 - #1)/2^(n + 1)}, {#1 + (#2 - #1)/2 + (#2 - #1)/2^(n + 1), #2}} & @@@b[n - 1], 1], {#1 + (#2 - #1)/2 - (#2 - #1)/2^(n + 2), #1 + (#2 - #1)/2 + (#2 - #1)/2^(n + 2)} & @@@ a[n - 1]] // Sort

a[2]
b[2]
(* {{0, 7/24}, {3/8, 2/3}, {19/24, 7/8}} *)
(* {{7/24, 3/8}, {2/3, 19/24}, {7/8, 1}} *)

These functions are memoized so that we eliminate any recursion depth problem. (That means that every time you changed the definitions, you'll need to Clear[a, b], as I've shown.)

To convert this to the interval form desired by the OP, we do

Prepend[#1 <= x < #2 & @@@ Rest@#, #1 <= x <= #2 & @@ First@#] &@a[2]
Prepend[#1 <= x < #2 & @@@ Rest@#, #1 <= x <= #2 & @@ First@#] &@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} *)