Question: Simplifying a large loop

Hey folks, I'll give you the actual problem I'm trying to solve before I show you my code in case anyone can think of a better idea...

Long read but you can ignore the examples and additional info if you don't need it.

Take the doubling map,

f(x) = 2x if 0 <= x <= 1/2

f(x) = 2x - 1 if 1/2 < x <= 1.

Or in maple code, f := x -> piecewise(x<=1/2, 2*x, 2*x - 1):

This looks something like this when graphed...

http://img413.imageshack.us/img413/2126/29431684.jpg

Now, divide the x-axis of the above graph into n equally spaced intervals where n is of the form 2^k (i.e. n=4,8,16,32, etc...). What I'm trying to do is find a periodic orbit such that each of these intervals contains exactly one point from this orbit. Perhaps you can think of it as being a histogram with uniform distribution if it helps..?

An example. Taking n=4, we have intervals [0,1/4], [1/4, 1/2], [1/2, 3/4], [3/4, 1]

If we take the initial point to be x=2/15.

Then if we repeatedly iterate this point using our map f we get...

x=2/15, f(x) = 4/15, f(f(x)) = 8/15, f(f(f(x))) = 1/15,  then back to 2/15.

It's easy to see that we have two points in the first interval (2/15 and 1/15) so this is not an orbit I'm after.

However take the point x=3/15.

Again, repeated iteration gives...

3/15, 6/15, 12/15, 9/15 then back to 3/15.

This orbit does have one point in each of the intervals above so it is one I'm looking for.

Now the question is how to find these orbits as 2^k gets pretty large quickly!

The way I've done this is to take my initial point as being a point in the first interval, this means I can narrow down where the point must be contained for it form this orbit. So for the 4 interval case, 0.01 will get mapped to 0.02 so we already have 2 points in the first interval, hence 0.01 would not be part of an orbit I want. Therefore the first point I take is in the second half of the first interval (i.e [1/8, 1/4] for the n=4 case).

Here is my code for this...

f := x -> piecewise(x<=1/2, 2*x, 2*x - 1):
n := 16:
K := NULL:
S := Array(1..n):

for i from (2^n)/(2*n) + 1 to (2^n)/n - 1 by 2 do
A := 1:
x[0] := i/(2^n - 1);
for j to n-1 do
S[1] := x[0]:
  x[j] := f(x[j-1]):
  S[j+1] := x[j]:
end do:
S := sort(S):
  for k from 2 to n do
    if S[k]<k/n and S[k]>(k-1)/n then
      A := A+1:
    end if:
  end do:
if A = n then
  K := K, x[0]
end if:
end do:

K := [K]:

for t in K do
  convert((2^n - 1)*t, binary)
end do;

Some explanation...

for k from 2 to n do
    if S[k]<k/n and S[k]>(k-1)/n then
      A := A+1:
    end if:
  end do:

This part is the part that's frankly a bit of a mess.

After S has been sorted this just checks each value in turn to make sure it's in the desired interval, if it is, then A := A+1. If every point is in the desired interval, A will be equal to n after the loop has finished and I store the value of x[0] for that cycle of the loop.

This may be pretty ugly but it does work. However, it doesn't work for n=32 (and above obviously) as it says the is/in object is too large. (I think that's what it says, didn't note it down when it failed earlier tonight)

So can anyone make this code any more efficient or think of another way to solve this?

Some more information that you may find useful.

If we write a point in [0,1] as x = a[1]*(1/2) + a[2]*(1/4) + a[3]*(1/8) + a[4]*(1/16) + ...

Then write this as 0.a[1]a[2]a[3]a[4]... where the a[k] = 0 or 1.

Then the mapping f (as defined above) is equivalent to a symbol shift of 0.a[1]a[2]a[3]a[4]...

I.e. if x = 0.a[1]a[2]a[3]a[4]...

f(x) = 0.a[2]a[3]a[4]...

f(f(x)) = 0.a[3]a[4]... etc

(Basically a binary symbol shift)

Any point of this interval orbit I'm after must have an equal number of 0s and 1s.

Dunno if this helps! If not disregard it but it might...