tirage <- function(){
# renvoie 1 si c'est gagné, 0 sinon
  return(runif(1)<1/2)
}

joueMM <- function(N, a){
  # joue en appliquant la martingale montante
  # N: fortune visée
  # a: fortune initiale
  X <- c(a);
  traj <- X;
  mise <-1;
  while (X>0 & X<N){
    if (tirage()){
      X <- X +  mise;
      mise <- 1;
    }
    else{
      X <- X-mise;
      mise <- min(X, 2*mise);
    }
    traj <- c(traj, X);
  }
  return(traj);
}

probaGagnerMM <- function(N, a, joue){
  #estime la proba de gagner avec la fonction joue
  # N: fortune visée
  # a: fortune initiale
  # joue: fonction qui joue une partie et renvoie la trajectoire des fortunes
  n = 1000;
  succes = 0;
  for (k in 1:n){
    traj <- joue(N,a);
    gagne <- (traj[length(traj)])==N;
    succes = succes + gagne;
  }
  return(succes / n);
}

probaGagnerRel <- function(N,a){
  # calcule la proba de gagner en utilisant la relation
  # N: fortune visée
  # a: fortune initiale
  M <- matrix(0, nrow = N+1, ncol = N+1);
  b <- array(0, N+1);
  b[N+1] <- 1;
  for (j in 2:N){
    x = j-1;
    k <- floor(log2(x+1)); #nb de coups qu'on peut perdre à la suite
    y <- x - (2^(k) - 1);
    M[j, j+1] <- 1-2^(-k);
    M[j, 1+y] <- 2^(-k);
  }
  p = solve(diag(N+1)-M, b);
  # print(p*N) # on remarque que p(a) = a/N !!
  return(p[1+a]);
}

N <- 126; # fortune visée
a <- 63; # fortune initiale

traj <- joueMM(N, a);
plot(traj, type='l');

p <- probaGagnerMM(N, a, joueMM);
print(p)
p <- probaGagnerRel(N,a)
print(p)