---
title: "Dyadic expansion with R"
date: "2016-08-08"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, collapse=TRUE, fig.path="./assets/fig/dyadic-", fig.align="center")
fscale <- 1 # scale for figures
```

We provide a function that computes the dyadic representation of a real number in the interval $[0,1]$. 
Then we give an implementation of two transformations of the set ${\{0,1\}}^\mathbb{N}$ which are well-known in ergodic theory: the dyadic odometer and the Pascal transformation. For each of these transformations, we plot the graph of the conjugate transformation of $[0,1]$ obtained by the dyadic representation.


## Dyadic expansion

Every real number $u \in [0,1]$ has a *dyadic expansion* (or *binary expansion*):
$$
u = \frac{\epsilon_1}{2} + \frac{\epsilon_2}{2^2} + \frac{\epsilon_3}{2^3} + \ldots
$$
where $\epsilon_i=0$ or $1$. 
We say that the sequence $(\epsilon_1, \epsilon_2, \ldots)$ is the *dyadic representation* of $u$. 

The `num2dyadic` function below returns the dyadic representation of $u \in [0,1]$.

```{r}
num2dyadic <- function(u, nmax=1024L){ 
  out <- integer(nmax)
  i <- j <- 0L
  while(u>0 && i<nmax){
    j <- 1L + max(0L,floor(-log2(u*(1+.Machine$double.eps^.5)))) 
    if(i+j <= nmax){
      i <- i + j
      out[i] <- 1L
      u <- 2L^j*u - 1L
    }else{
      i <- nmax
    }
  }
  return(out[1:i])
}
```

The `dyadic2num` function below does the reverse action:

```{r}
dyadic2num <- function(d) sum(d/2L^(seq_along(d))) 
```

Let us check that the dyadic representation of $0.75 = \frac{1}{2}+\frac{1}{4}$ is $(1,1)$:

```{r}
num2dyadic(1/2+1/4)
```

The real number $u=0.2$ has the infinite periodic dyadic representation $(0, 0, 1, 1, 0, 0, 1, 1, \ldots)$. 
The `num2dyadic` function applied to $0.2$ returns only the first $54$ digits of its dyadic representation:

```{r}
( d <- num2dyadic(0.2) )
length(d)
```

But it makes no difference for R:

```{r}
dyadic2num(d) == 0.2
```


## The dyadic odometer 

The dyadic odometer is the transformation $O$ of the set ${\{0,1\}}^{\mathbb{N}}$ defined by $O(d) = d + (1, 0, 0, \ldots)$, where "$+$" is the addition $\bmod\, 2$ with carry to the right. 

The `odometer` function below is an implementation of the dyadic odometer and its inverse (option `image="backward"`).

```{r odometer}
odometer <- function(d, image=c("forward", "backward")){
  image <- match.arg(image)
  if(image=="forward"){
    if(all(d==1L)){
      d <- c(rep(0L, length(d)), 1L)
    }else{
      k <- which.min(d)
      d[1:k] <- 1L-d[1:k]
    }
  }
  if(image=="backward"){
    if(all(d==0L)){
      d <- c(rep(1L, length(d)), 0L)
    }else{
      k <- which.max(d)
      d[1:k] <- 1L-d[1:k]
    }
  }
  return(d)
}
```

Using the dyadic representation, the odometer also defines a map from the interval $[0,1]$ to itself. We plot its graph below:

```{r plot_odometer, fig.width=fscale*4, fig.height=fscale*4}
par(mar=c(4,4,2,2))
u <- seq(0, 0.995, by=0.005)
Ou <- sapply(u, function(u) dyadic2num(odometer(num2dyadic(u))))
plot(u, Ou, xlab="u", ylab="O(u)", 
     xlim=c(0,1), ylim=c(0,1),
     pch=19, cex=.25, 
     pty="s", xaxs="i", yaxs="i")
grid(nx=10)
```


## The Pascal transformation 

The Pascal transformation $P$ is defined for every $d \in {\{0,1\}}^{\mathbb{N}}$ 
except when $d=000\ldots$ or when $d$ has the form $d=0^i111\ldots$ ($i\geq 0$). Such a $d$ has the form $d= 0^m1^k10x$ where $m,k \geq 0$ and $x \in {\{0,1\}}^{\mathbb{N}}$, and then the image of $d$ by the Pascal transformation is
$$
P(0^m1^k10x) = 1^k0^m01x.
$$
The case when $d=0^i111\ldots$ does not occur for us since we deal with finite sequences only. One naturally extends the Pascal transformation to include the case $d=000\ldots$ by setting $P(000\ldots) = 000\ldots$.

```{r pascal}
pascal <- function(d){
  if(all(d==0L)) return(0L)
  i <- which.max(d)
  m1 <- i-1L
  d0 <- c(d, 0L)
  k1 <- which.min(d0[-(1:i)])-1L
  begin <- c(rep(1L, k1), rep(0L, m1+1L), 1L)
  if(length(d)==m1+k1+1L) d <- begin else d[1L:(m1+k1+2L)] <- begin
  return(d)
}
```

By the dyadic representation, the Pascal transformation also defines a map from the interval $[0,1)$ to itself, whose graph is plotted below:

```{r plot_pascal, fig.width=fscale*4, fig.height=fscale*4}
par(mar=c(4,4,2,2))
u <- seq(0, 1-1/2^10, by=1/2^10)
Pu <- sapply(u, function(u) dyadic2num(pascal(num2dyadic(u))))
plot(u, Pu, xlab="u", ylab="P(u)", 
     xlim=c(0,1), ylim=c(0,1),
     pch=19, cex=.25, 
     pty="s", xaxs="i", yaxs="i")
grid(nx=10)
```