--- a/main/queueing/doc/references.txi	Tue Jan 07 21:23:53 2014 +0000
+++ b/main/queueing/doc/references.txi	Wed Jan 08 14:18:42 2014 +0000
@@ -129,7 +129,11 @@
 @item [ZaWo81]
 J. Zahorjan and E. Wong, @cite{The solution of separable queueing
 network models using mean value analysis}. SIGMETRICS
-Perform. Eval. Rev. 10, 3 (Sep. 1981), 80-85. DOI
+Perform. Eval. Rev. 10, 3 (Sep. 1981), 80-85.
 DOI @uref{http://doi.acm.org/10.1145/1010629.805477, 10.1145/1010629.805477}
 
+@item [Zeng03]
+G. Zeng, @cite{Two common properties of the erlang-B function, erlang-C function, and Engset blocking function}, Mathematical and Computer Modelling, Volume 37, Issues 12-13, June 2003, Pages 1287-1296 DOI
+@uref{http://dx.doi.org/10.1016/S0895-7177(03)90040-9, 10.1016/S0895-7177(03)90040-9}
+
 @end table

--- a/main/queueing/doc/singlestation.txi	Tue Jan 07 21:23:53 2014 +0000
+++ b/main/queueing/doc/singlestation.txi	Wed Jan 08 14:18:42 2014 +0000
@@ -114,6 +114,10 @@
 
 @GETHELP{qsmmm}
 
+@GETHELP{erlangb}
+
+@GETHELP{erlangc}
+
 @noindent @strong{REFERENCES}
 
 @noindent G. Bolch, S. Greiner, H. de Meer and K. Trivedi, @cite{Queueing Networks
@@ -125,6 +129,10 @@
 @auindex de Meer, H.
 @auindex Trivedi, K.
 
+@noindent G. Zeng, @cite{Two common properties of the erlang-B function, erlang-C function, and Engset blocking function}, Mathematical and Computer Modelling, Volume 37, Issues 12-13, June 2003, Pages 1287-1296
+
+@auindex Zeng, G.
+
 @c
 @c M/M/inf
 @c
@@ -156,11 +164,12 @@
 @node The M/M/1/K System
 @section The @math{M/M/1/K} System 
 
-In a @math{M/M/1/K} finite capacity system there can be at most
-@math{k} jobs at any time. If a new request tries to join the system
-when there are already @math{K} other requests, the arriving request
-is lost. The queue has @math{K-1} slots. The @math{M/M/1/K} system is
-always stable, regardless of the arrival and service rates
+In a @math{M/M/1/K} finite capacity system there is a single server
+and there can be at most @math{k \geq 1} jobs at any time (including
+the job currently in service). If a new request tries to join the
+system when there are already @math{K} other requests, the arriving
+request is lost. The queue has @math{K-1} slots. The @math{M/M/1/K}
+system is always stable, regardless of the arrival and service rates
 @math{\lambda} and @math{\mu}.
 
 @GETHELP{qsmm1k}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/erlangb.m	Wed Jan 08 14:18:42 2014 +0000
@@ -0,0 +1,116 @@
+## Copyright (C) 2014 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{B} =} erlangb (@var{A}, @var{m})
+##
+## @cindex Erlang-B formula
+##
+## Compute the value of the Erlang-B formula @math{E_B(A, m)} giving the
+## probability that a request arriving to a pool of @math{m} identical
+## servers is rejected. 
+## 
+## @iftex
+##
+## @math{E_B(A, m)} is defined as:
+##
+## @tex
+## $$
+## E_B(A, m) = \displaystyle{{A^m \over m!} \left( \sum_{k=0}^m {A^k \over k!} \right) ^{-1}}
+## $$
+## @end tex
+##
+## @end iftex
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item A
+## Offered load, defined as @math{A = \lambda / \mu} where
+## @math{\lambda} is the mean arrival rate and @math{\mu} the mean
+## service rate of each individual server (real, @math{A > 0}).
+##
+## @item m
+## Number of identical servers (integer, @math{m @geq{} 1}). Default @math{m = 1}
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item B
+## The value @math{E_B(A, m)}
+##
+## @end table
+##
+## @var{A} or @var{m} can be vectors, and in this case, the results will
+## be vectors as well.
+##
+## @seealso{qsmmm}
+##
+## @end deftypefn
+
+## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
+## Web: http://www.moreno.marzolla.name/
+function B = erlangb(A, m)
+  if ( nargin < 1 || nargin > 2 )
+    print_usage();
+  endif
+
+  ( isnumeric(A) && all( A(:) > 0 ) ) || error("A must be positive");
+  
+  if ( nargin == 1 )
+    m = 1;
+  else
+    ( isnumeric(m) && all( fix(m(:)) == m(:)) && all( m(:) > 0 ) ) || error("m must be a positive integer");
+  endif
+
+  [err A m] = common_size(A, m);
+  if ( err )
+    error("parameters are not of common size");
+  endif
+
+  B = arrayfun( @__erlangb_compute, A, m);
+endfunction
+
+## Compute E_B(A,m) recursively, as described in Guoping Zeng (June
+## 2003), Two common properties of the erlang-B function, erlang-C
+## function, and Engset blocking function, Elsevier Science, retrieved
+## 2011-02-03
+##
+## To improve numerical stability, the recursion is based on the inverse
+## of E_B rather than E_B itself.
+function B = __erlangb_compute(A, m)
+  Binv = 1.0;
+  for k=1:m
+    Binv = 1.0 + k/A*Binv;
+  endfor
+  B = 1.0 / Binv;
+endfunction
+
+%!test
+%! fail("erlangb(1, -1)", "positive");
+%! fail("erlangb(-1, 1)", "positive");
+%! fail("erlangb(1, 0)", "positive");
+%! fail("erlangb(0, 1)", "positive");
+%! fail("erlangb('foo',1)", "positive");
+%! fail("erlangb(1,'bar')", "positive");
+%! fail("erlangb([1 1],[1 1 1])","common size");
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/erlangc.m	Wed Jan 08 14:18:42 2014 +0000
@@ -0,0 +1,102 @@
+## Copyright (C) 2014 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{C} =} erlangc (@var{A}, @var{m})
+##
+## @cindex Erlang-C formula
+##
+## Compute the steady-state probability that a request finds all
+## @math{m} servers musy in a @math{M/M/m/\infty} queueing system with
+## @math{m} identical servers and infinite queue. This quantity is
+## also called Erlang-C formula @math{E_C(A, m)}.
+## 
+## @iftex
+##
+## @math{E_C(A, m)} is defined as:
+##
+## @tex
+## $$
+## E_C(A, m) = \displaystyle{ {A^m \over m!} {1 \over 1-\rho} \left( \sum_{k=0}^{m-1} {A^k \over k!} + {A^m \over m!} {1 \over 1 - \rho} \right) ^{-1}}
+## $$
+## @end tex
+##
+## where @math{\rho = A / m = \lambda / (m \mu)}.
+##
+## @end iftex
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item A
+## Offered load, defined as @math{A = \lambda / \mu} where
+## @math{\lambda} is the mean arrival rate and @math{\mu} the mean
+## service rate of each individual server (real, @math{0 < A < m}).
+##
+## @item m
+## Number of identical servers (integer, @math{m @geq{} 1}). Default @math{m = 1}
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item B
+## The value @math{E_C(A, m)}
+##
+## @end table
+##
+## @var{A} or @var{m} can be vectors, and in this case, the results will
+## be vectors as well.
+##
+## @seealso{qsmmm}
+##
+## @end deftypefn
+
+## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
+## Web: http://www.moreno.marzolla.name/
+function C = erlangc(A, m)
+  if ( nargin < 1 || nargin > 2 )
+    print_usage();
+  endif
+
+  ( isnumeric(A) && all(A(:) > 0) ) || error("A must be positive");
+   
+  if ( nargin == 1 )
+    m = 1;
+  else
+    ( isnumeric(m) && all( fix(m(:)) == m(:)) && all( m(:) > 0 ) ) || error("m must be a positive integer");
+  endif
+  
+  [err A, m] = common_size(A, m);
+  if ( err )
+    error("parameters are not of common size");
+  endif
+  
+  all( A(:) < m(:) ) || error("A must be < m");
+   
+  rho = A ./ m;
+  B = erlangb(A, m);
+  C = B ./ (1 - rho .* (1 - B));
+endfunction
+
+%!test
+%! fail("erlangc('foo',1)", "positive");
+%! fail("erlangc(1,'bar')", "positive");

--- a/main/queueing/inst/qnjackson.m	Tue Jan 07 21:23:53 2014 +0000
+++ b/main/queueing/inst/qnjackson.m	Wed Jan 08 14:18:42 2014 +0000
@@ -21,89 +21,9 @@
 ## @deftypefnx {Function File} {[@var{U}, @var{R}, @var{Q}, @var{X}] =} qnjackson (@var{lambda}, @var{S}, @var{P}, @var{m} )
 ## @deftypefnx {Function File} {@var{pr} =} qnjackson (@var{lambda}, @var{S}, @var{P}, @var{m}, @var{k})
 ##
-## @cindex open network, single class
-## @cindex Jackson network
-##
-## This function is deprecated. Please use @code{qnopensingle} instead.
-##
-## With three or four arguments, this function computes the steady-state
-## occupancy probabilities for a Jackson network. With five arguments,
-## this function computes the steady-state probability
-## @code{@var{pi}(j)} that there are @code{@var{k}(j)} requests at
-## service center @math{j}.
-##
-## This function solves a subset of Jackson networks, with the
-## following constraints:
-##
-## @itemize
-##
-## @item External arrival rates are load-independent.
-## 
-## @item Service center @math{i} consists either of @code{@var{m}(i) @geq{}
-## 1} identical servers with individual average service time
-## @code{@var{S}(i)}, or of an Infinite Server (IS) node.
-##
-## @end itemize
-##
-## @strong{INPUTS}
-##
-## @table @var
-##
-## @item lambda
-## @code{@var{lambda}(i)} is
-## the external arrival rate to service center @math{i}. @var{lambda}
-## must be a vector of length @math{N}, @code{@var{lambda}(i) @geq{} 0}.
-##
-## @item S
-## @code{@var{S}(i)} is the average service time on service center @math{i}
-## @var{S} must be a vector of length @math{N}, @code{@var{S}(i)>0}.
+## This function is deprecated. Please use @code{qnos} instead.
 ##
-## @item P
-## @code{@var{P}(i,j)} is the probability
-## that a job which completes service at service center @math{i} proceeds
-## to service center @math{j}. @var{P} must be a matrix of size
-## @math{N \times N}.
-##
-## @item m
-## @code{@var{m}(i)} is the number of servers at center
-## @math{i}. If @code{@var{m}(i) < 1}, center @math{i} is an
-## infinite-server node. Otherwise, it is a regular FCFS queueing center with
-## @code{@var{m}(i)} servers. Default is 1.
-##
-## @item k
-## Compute the steady-state probability that there are @code{@var{k}(i)}
-## requests at service center @math{i}. @var{k} must have the same length
-## as @var{lambda}, with @code{@var{k}(i) @geq{} 0}.
-## 
-## @end table
-##
-## @strong{OUTPUT}
-##
-## @table @var
-##
-## @item U
-## If @math{i} is a FCFS node, then
-## @code{@var{U}(i)} is the utilization of service center @math{i}.
-## If @math{i} is an IS node, then @code{@var{U}(i)} is the
-## @emph{traffic intensity} defined as @code{@var{X}(i)*@var{S}(i)}.
-##
-## @item R
-## @code{@var{R}(i)} is the average response time of service center @math{i}.
-##
-## @item Q
-## @code{@var{Q}(i)} is the average number of customers in service center
-## @math{i}.
-##
-## @item X
-## @code{@var{X}(i)} is the throughput of service center @math{i}.
-##
-## @item pr
-## @code{@var{pr}(i)} is the steady state probability 
-## that there are @code{@var{k}(i)} requests at service center @math{i}.
-##
-## @end table
-##
-## @seealso{qnopen}
+## @seealso{qnos}
 ##
 ## @end deftypefn
 
@@ -115,7 +35,7 @@
   if (!warned)
     warned = true;
     warning("qn:deprecated-function",
-	    "qnjackson is deprecated. Please use qnopensingle insgead");
+	    "qnjackson is deprecated. Please use qnos insgead");
   endif
   if ( nargin < 3 || nargin > 5 )
     print_usage();
@@ -194,35 +114,3 @@
 
   endif
 endfunction
-%!test
-%! # Test various error conditions
-%! fail( "qnjackson( [0.5 0.5], [0 1], [0 0; 0 0], [1 1])", "S must be" );
-%! fail( "qnjackson( [-1 1], [1 1], [0 0; 0 0], [1 1])", "lambda must be" );
-%! fail( "qnjackson( [0.5 0.5], [1 1], [1 1; 0 0], [1 1])", "P is not" );
-%! fail( "qnjackson( [0.5 0.5], [1 1], [1 0; -1 0], [1 1])", "P is not" );
-%! fail( "qnjackson( [0.5 0.5], [1 1 1], [0 0; 0 0], [1 1])", "lambda and S" );
-%! fail( "qnjackson( [0.5 0.5], [1 1], [0 0; 0 0], [1 1 1])", "m and lambda" );
-%! fail( "qnjackson( [0.5 0.5], [1 1], [0 0; 0 0], [1 1], [1 1 1])", "k must be" );
-%! fail( "qnjackson( [0.5 0.5], [1 1], [0 0; 0 0], [1 1], [1 -1])", "k must be" );
-%! fail( "qnjackson( [0 1], [2 2], [0 0.9; 0 0] )", "not ergodic" );
-
-%!test
-%! # Example 7.4 p. 287 Bolch et al.
-%! S = [ 0.04 0.03 0.06 0.05 ];
-%! P = [ 0 0.5 0.5 0; 1 0 0 0; 0.6 0 0 0; 1 0 0 0 ];
-%! lambda = [0 0 0 4];
-%! k = [ 3 2 4 1 ];
-%! [U R Q X] = qnjackson( lambda, S, P, 1 );
-%! assert( X, [20 10 10 4], 1e-4 );
-%! assert( U, [0.8 0.3 0.6 0.2], 1e-2 );
-%! assert( R, [0.2 0.043 0.15 0.0625], 1e-3 );
-%! assert( Q, [4, 0.429 1.5 0.25], 1e-3 );
-
-%!test
-%! # Example 7.4 p. 287 Bolch et al.
-%! S = [ 0.04 0.03 0.06 0.05 ];
-%! P = [ 0 0.5 0.5 0; 1 0 0 0; 0.6 0 0 0; 1 0 0 0 ];
-%! lambda = [0 0 0 4];
-%! k = [ 3 2 4 1 ];
-%! p_i = qnjackson( lambda, S, P, 1, k );
-%! assert( p_i, [0.1024, 0.063, 0.0518, 0.16], 1e-4 );

--- a/main/queueing/inst/qsmmm.m	Tue Jan 07 21:23:53 2014 +0000
+++ b/main/queueing/inst/qsmmm.m	Wed Jan 08 14:18:42 2014 +0000
@@ -1,4 +1,4 @@
-## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013 Moreno Marzolla
+## Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2014 Moreno Marzolla
 ##
 ## This file is part of the queueing toolbox.
 ##
@@ -22,7 +22,10 @@
 ##
 ## @cindex @math{M/M/m} system
 ##
-## Compute utilization, response time, average number of requests in service and throughput for a @math{M/M/m} queue, a queueing system with @math{m} identical service centers connected to a single FCFS queue.
+## Compute utilization, response time, average number of requests in
+## service and throughput for a @math{M/M/m} queue, a queueing system
+## with @math{m} identical service centers connected to a single FCFS
+## queue.
 ##
 ## @iftex
 ## The steady-state probability @math{\pi_k} that there are @math{k}
@@ -93,7 +96,7 @@
 ## @var{lambda}, @var{mu} and @var{m} can be vectors of the same size. In this
 ## case, the results will be vectors as well.
 ##
-## @seealso{qsmm1,qsmminf,qsmmmk}
+## @seealso{erlangc,qsmm1,qsmminf,qsmmmk}
 ##
 ## @end deftypefn
 
@@ -106,38 +109,33 @@
   endif
   if ( nargin == 2 )
     m = 1;
+  else
+    ( isnumeric(lambda) && isnumeric(mu) && isnumeric(m) ) || ...
+	error( "the parameters must be numeric vectors" );
   endif
-  ( isvector(lambda) && isvector(mu) && isvector(m) ) || ...
-      error( "the parameters must be vectors" );
   [err lambda mu m] = common_size( lambda, mu, m );
   if ( err ) 
     error( "parameters are not of common size" );
   endif
   lambda = lambda(:)';
   mu = mu(:)';
+  m = m(:)';
   all( m>0 ) || ...
       error( "m must be >0" );
   all( lambda>0 ) || ...
       error( "lambda must be >0" );
-  all( lambda < m .* mu ) || ...
-      error( "Processing capacity exceeded" );
   X = lambda;
   U = rho = lambda ./ (m .* mu );
+  all( U < 1 ) || ...
+      error( "Processing capacity exceeded" );
   Q = p0 = pm = 0*lambda;
   for i=1:length(lambda)
-#{
-    k=[0:m(i)-1];
-    p0(i) = 1 / ( ...
-		 sum( (m(i)*rho(i)).^ k ./ factorial(k)) + ...
-                 (m(i)*rho(i))^m(i) / (factorial(m(i))*(1-rho(i))) ...
-                 );
-#}
     p0(i) = 1 / ( ...
                  sumexpn( m(i)*rho(i), m(i)-1 ) + ...		 
 		 expn(m(i)*rho(i), m(i))/(1-rho(i)) ...
                  );
-    pm(i) = expn(m(i)*rho(i),m(i))*p0(i)/(1-rho(i));
   endfor
+  pm = erlangc(lambda ./ mu, m);
   Q = m .* rho .+ rho ./ (1-rho) .* pm;
   R = Q ./ X;
 endfunction