%%%%%%%%%%%%%%%%%%%%%
%Delhaye Quentin et Garcia Duarte Eduardo
%Polytech BA2 - MATH-H-204
%%%%%%%%%%%%%%%%%%%%%

%%%
%Situation : un serpent du classique Snake s'allonge depuis sa position de
%départ. Toutes les directions son équiprobables et il ne peut pas revenir
%sur ses pas.
%Problème : quelle est la probabilité qu'il fasse plus de 23 mouvements
%sans se rentrer dedans ?
%%%

format long
%On cherche la proba qu'il se touche après 24 mouvements ou plus => on va 
%calculer la probabilité complémentaire (proba de se toucher endéans 23  
%mouvements qu'on va soustraire à  l'unité)
max=23;

%Nombre d'expériences
N=200000;
%histDir(1) = direction opposée au mouvement (n-1) et 
%histDir(2) = direction du mouvement n.
histDir=zeros(2,1);
A=zeros(N,1);
valMoyenne=zeros(N,1);
nbTouche=0;
tic
for i=1:N
    %Coord : On a deux colonnes : x et y, contenant les coordonnées 
    %successives du bout du serpent.
    Coord=zeros(max+1,2);
    touche=0;
    for j=1:max
        %Il ne peut pas retourner sur lui-même. randi(4,1) retourne un
        %chiffre entre 1 et 4 de manière équiprobable.
        while histDir(1) == histDir(2)
            histDir(2) = randi(4,1);
        end
        switch histDir(2)
            case 1 %gauche
                Coord(j+1,1)=Coord(j,1)-1;
                Coord(j+1,2)=Coord(j,2);
                %S'il se déplace à  gauche, le mouvement suivant ne peut 
                %pas se faire vers la droite, droite qui correspond à la
                %direction 3.
                histDir(1)=3;
            case 2 %haut
                Coord(j+1,1)=Coord(j,1);
                Coord(j+1,2)=Coord(j,2)+1;
                histDir(1)=4;
            case 3 %droite
                Coord(j+1,1)=Coord(j,1)+1;
                Coord(j+1,2)=Coord(j,2);
                histDir(1)=1;
            case 4 %bas
                Coord(j+1,1)=Coord(j,1);
                Coord(j+1,2)=Coord(j,2)-1;
                histDir(1)=2;
        end
        %Comparaison des coordonnées actuelles avec celles déjà parcourues
        for k=1:j
            if [Coord(j+1,1) Coord(j+1,2)] == [Coord(k,1) Coord(k,2)]
                touche=1;
            end
        end
        %S'il se touche, c'est terminé
        if touche == 1
            nbTouche=nbTouche+1;
            break
        end
        histDir(2)=histDir(1); %on lui indique d'où il vient pour qu'il 
        %puisse rentrer dans la boucle et éliminer la direction idoine
    end
    %Fréquence relative qu'il ne se touche pas en 23 coups
    valMoyenne(i,1)=1-(nbTouche/i);
end
toc
grid on
plot([1:N],valMoyenne,'r','LineWidth',2)
freqRel = 1 - (nbTouche/N)

%%%
%Output :
% >> snake
% Elapsed time is 84.104369 seconds.
% 
% freqRel =
% 
%    0.136840000000000
%%
%Remarque : Ce résultat est obtenu en une durée relativement longue
%(compte tenu des perfomances de la machine de calcul). Ceci peut
%s'expliquer par la comparaison systématique de la coordonnée courante avec
%les précédentes. Un code plus efficace aurait pu être produit en utilisant
%plutôt un tableau de booléens initialisé false et dont les cases
%parcourues deviendraient true.
%%%