function [p, mNa, mNd]=risk(Na,Nd,Nsim)

## RISK GNU OCtave/MATLAB function för calculation of the battle outcome probabilities of the game of risk by stochastic simulation, or that shows the outcome of one random battle (one sequence of dice rolls).
## Call: [p, mNa, mNd]=risk(Na,Nd,Nsim)
## Example: p = risk(1:3,1:3,5000)
## Example: risk(4, 5, 1)
##
## Input parameters: 
## Na is the number of attacking armies. (Scalar, or vector of several cases).
## Nd is the number of defending armies. (Scalar, or vector of several cases).
## Nsim is the number of randomized simulations. 
## If Nsim is 1, the outcome of that battle is shown (number of remaining armies). 
## This is useful in view to make a real game faster, but maybe not as fun... 
%%
## Output parameters: p is the probability that the attacker wins. 
## If several cases, i.e. if Na and/Nd Nd are vectors, p is a matrix,
## where the element p(Ndindx,Naindx) is the probabilitity that the attacker wins 
## when Na(Naindx) attackers and Nd(Ndindx) försvarare.
%% mNa(Ndindx,Naindx) is the average number of remaning attackers.
%% mNd(Ndindx,Naindx) is the average number of remaning defenders.

## Author: Mange01 <magnus.eriksson@miun.se> 
## Created: 2010-07-04

if length(Na)>1, % Several cases
	for k=1:length(Na)
		disp(['Antal anfallare: ' int2str(Na(k))])
		[p(:,k), mNa(:,k), mNd(:,k)]=risk(Na(k),Nd,Nsim); % Självanrop för varje fall.
	end

elseif length(Nd)>1
	for k=1:length(Nd)
		disp(['Antal försvarare: ' int2str(Nd(k))])
		[p(k,:), mNa(k,:), mNd(k,:)]=risk(Na,Nd(k),Nsim); % Självanrop för varje fall.
	end
else % Only one Na and Nd case

	for k=1:Nsim
		[Vinnare(k),NA(k),ND(k)] = risk1(Na,Nd);
	end
	p=mean(Vinnare);
	mNa=mean(NA);
	mNd=mean(ND);
	disp('')
	if Nsim==1, % Display the result
		if Vinnare(1)==0
			disp(['The defender wins and keep ' int2str(ND(1))])
		else
			disp(['The attacker wins and keep ' int2str(NA(1))])
		end
	else
%	disp(['Probability that the attacker wins:  ' num2str(p*100) ' %.'])
%	disp('')
%	disp(['Frequency chart for the number of remaining nodes:'])
%	values = unique([NA ND]);
%	fNA = hist(NA, values);
%	fND = hist(ND, values);
%	disp([values;[fNA;fND]/Nsim])
	
	end % else
end % else
endfunction

function [winflag,Na,Nd] = risk1 (Na,Nd)
## Sub-function that is called for one randomized simulation run of one battle.
while Na>0 && Nd>0, % Repete tills en har slut på armeer.
	a=dice(min(Na,3));
	f=dice(min(Nd,2));
	if Na>1 && Nd>1, % Om båda slår med minst två tärningar
		if a(2)>f(2), % Jämför de näst största tärningsvärdena.
			Nd = Nd-1; % Om Anfallaren fick högre resultat, minska försvararens armeer.
			else
			Na = Na-1; % Om anfallaren fick lika eller lägre, minska försvararens armeer.
		end

	end
	if a(1)>f(1) % Jämför de största tärningsvärdena.
		Nd = Nd-1;
	else
		Na = Na-1;
	end
end % while

if Na>0
	winflag=1; % Attacker wins
elseif Nd>0
	winflag=0; % Defender wins
else
	error('Programfel')
end
endfunction

function y=dice(N)
## Subfunctions that return N random die values. The result is sorted - the largest dice comes first.
	y=sort(ceil(rand(1,N)*6));
	y=y(end:-1:1);
endfunction