// Calcul de probabilités de jets de dés par dénombrement

nfaces = 6; // nombre de faces du dé
nfaces2 = nfaces*nfaces;
nfaces3 = nfaces2*nfaces;
nfaces4 = nfaces3*nfaces;

de = 1:nfaces;

// 3d6 parmi 3

de13 = matrix(ones(nfaces2,1)*de, nfaces3, 1);
de23 = matrix(ones(nfaces, 1)*matrix(de'*ones(1, nfaces), 1, nfaces2), nfaces3, 1);
de33 = matrix(de'*ones(1, nfaces2), nfaces3, 1);

de3 = [de13 de23 de33];

freq_3d = tabul(sum(de3, "c"));
freq_3d(:, 2) = 100/nfaces3*freq_3d(:, 2);

// 3d6 parmi 4
de14 = matrix(ones(nfaces3,1)*de, nfaces4, 1);
de24 = matrix(matrix(ones(nfaces2, 1)*de, nfaces3, 1)*ones(1, nfaces), nfaces4, 1);
de34 = matrix(ones(nfaces, 1)*matrix(de'*ones(1, nfaces2), 1, nfaces3), nfaces4, 1);
de44 = matrix(de'*ones(1, nfaces3), nfaces4, 1);

de4=[de14 de24 de34 de44];

[de4min, indexmin] = min(de4, "c");
freq_4d = tabul(sum(de4, "c") - de4min);
freq_4d(:, 2) = 100/nfaces4*freq_4d(:, 2);

scf(0);
clf();

subplot(2, 1, 1)
plot2d(freq_3d(:, 1), freq_3d(:, 2), style=13);
p1 = gce();
p1.children.mark_style=5;
p1.children.mark_foreground=13;
plot2d(freq_4d(:, 1), freq_4d(:, 2), style=2);
p2 = gce();
p2.children.mark_style=6;
p2.children.mark_foreground=2;
xtitle("3d6 (vert) et les trois meilleurs dés de 4d6 (bleu)",...
    "résultat des dés", "probabilité (%)")

subplot(2, 1, 2)
plot2d(freq_3d(:, 1), 100-cumsum(freq_3d(:, 2)), style=13);
p1 = gce();
p1.children.mark_style=5;
p1.children.mark_foreground=13;
plot2d(freq_4d(:, 1), 100-cumsum(freq_4d(:, 2)), style=2);
p2 = gce();
p2.children.mark_style=6;
p2.children.mark_foreground=2;
p3=gca();
p3.tight_limits(2) = "on";
xtitle("3d6 (vert) et les trois meilleurs dés de 4d6 (bleu)",...
    "résultat des dés", "probabilité (%)")
xtitle("Fréquences cumulées",...
    "résultat des dés", "probabilité (%)")