# Nombre de résultats calculés
# de -N à +N
N = 20;
# Taille du dé
M = 10;
# Nombre de résultats affichés
# de -R à +R
R = 20;
# Nombre de dés
J = 2;
# Nombre de résultats affichés
# de -R à +R
R = 20;
# Prob de 1 à N
P1N = 1:N;
# Prob de -N à N
PDN = 0:2*N;

NSG = 0:2*N;

function y = P1_ferme(x,M,J)
  if ((x < 1) || (x > M))
    y = 0;
  else
    y = 1/M;
  end;
end;

function y = P1_limite(x,M,J)
  if (x < 1)
    y = 0;
  elseif (x < M)
    y = 1/M;
  elseif (x < J*M)
    y = P1_limite(x-M,M,J)/M;
  elseif (x == J*M)
    y = (1/M)^J;
  else
    y = 0;	
  end;
end;

function y = P1_ouvert(x,M,J)
  if (x < 1)
    y = 0;
  elseif (x < M)
    y = 1/M;
  else
    y = P1_ouvert(x-M,M)/M;
  end;
end;

function y = P1(x,M,J)
  y = P1_limite(x,M,J);
end;

function y = diff(x,P1N,N)
  if (x < 0)
    c = -x;
  else
    c = x;
  end

  if (c >= N)
    y = 0;
    return;
  end

  y = 0;
  for i = (c+1):N
    y = y + P1N(i) * P1N(i-c);
  end;
end;

function y = prob_real(N, M, J)
# usage: y = prob_real(N, M, J)
# Return a -N:N vector containing
# the compute value of probability
# for J M-side dices.
  for i = 1:N
    p1(i) = P1(i,M,J);
  end;
  for i = -N:N
    y(i+N+1) = diff(i,p1,N);
  end;
end;

function y = prob_theo(N)
# usage: y = prob_theo(N)
# Return a -N:N vector containing
# the theorical value of probability
  for i = -N:N
    if (i < 0)
      y(i+N+1) = 0.5 * 2^(i/3);
    else
      y(i+N+1) = 1 - 0.5 * 2^(-i/3);
    end;
  end;
end;

PDN = prob_real(N,M,J);
IPDN = cumsum(PDN);
NSG = prob_theo(N);

RANGE=-R:R;
IRANGE=N-R+1:N+R+1;

# Affichage des résultats
IPDN(IRANGE)

title("courbe de probabilité");
plot(RANGE, IPDN(IRANGE), "-.", RANGE, NSG(IRANGE), "-");
r = scanf("%s","C");
title("erreur absolue");
plot(RANGE, abs(IPDN(IRANGE)-NSG(IRANGE))*100);
r = scanf("%s","C");
title("erreur relative");
plot(RANGE, abs(IPDN(IRANGE)-NSG(IRANGE))./NSG(IRANGE)*100);
r = scanf("%s","C");

function y = diff_sum(P,R)
  Psize = size(P,1)
  Rsize = size(R,1)
  for i = 1:Rsize-1
    off = R(i+1)-1
    y(i) = sum(P(R(i):off))
  end;
  y(Rsize) = sum(P(R(Rsize):Psize))
end;

X = PDN(N:-1:1);
X = X / sum(X);
S = [1;5;10;15]
diff_sum(X,S)