// CS 1538 Fall 2009
// Program to solve Exercise 6.21 in the Banks text
// Note how even with a program there are some difficult issues in solving
// these problems -- most notably with number size and precision.

public class Ex621
{
	double lambda = 1000.0;
	double mu = 1.0/3.0;
	
	public Ex621()
	{
		// Approach 1:
		// Initially, we don't really have any idea of how many spaces (servers) we need.
		// Thus, we need a fast way of closing in on the answer.  One reasonable approach
		// is to double the value until we exceed our goal, then use a binary search to
		// narrow the answer down.  Since lambda is 1000 and mu is 1/3, a reasonable 
		// place to start is with lambda/mu = 3000 spots
		int c = 3000;
		int old = 1500;
		double p_no_fit = 1 - Pn(c);
		System.out.println("1 - P(" + c + ") = " + p_no_fit);
		while (p_no_fit > 0.001)
		{
			old = c;
			c = 2 * c;
			p_no_fit = 1 - Pn(c);
			System.out.println("1 - P(" + c + ") = " + p_no_fit);
		}

		boolean done = false;
		int low = old, high = c, mid = low;
		while (low < high)
		{
			mid = (low + high)/2;
			p_no_fit = 1 - Pn(mid);
			System.out.println("1 - P(" + mid + ") = " + p_no_fit);
			if (p_no_fit > 0.001)
				low = mid;
			else
			{
				high = mid;
				if (high - low <= 1)
					break;
			}

		}
		System.out.println("A total of " + mid + " spaces are needed ");

		// Approach 2:
		// Another (arguably better) approach is to start with 1.0 and subtract
		// the P value for each i (starting at 0) until we get below our desired
		// value (of 0.001).  See the P2 method below.
		int ans2 = P2(0.001);
		System.out.println("\nAnother way " + ans2 + " spaces are needed");
	}


	public static void main (String [] args)
	{
		Ex621 temp = new Ex621();
	}

	// Calculate F(n) (prob. of there being n or fewer "customers" 
	// in the system)
	public double Pn(int n)
	{
		double p_tot = 0;
		for (int i = 0; i <= n; i++)
		{
			double pcurr = Pi(i);
			//System.out.println("P " + i + ":" + pcurr);
			p_tot += pcurr;
		}
		//System.out.println("Pi = " + p_tot);
		return p_tot;
	}

	// Calculate number of customers required such that
	// 1 - F(n) ~= pct (the parameter)
	public int P2(double pct)
	{
		double p_tot = 1;
		int i;
		for (i = 0; p_tot > pct; i++)
		{
			double pcurr = Pi(i);
			p_tot -= pcurr;
			System.out.println("P " + i + ":" + pcurr +
				"  P(N <= i): " + p_tot);
		}
		return i-1;
	}

	// Both methods for determining the answer above require the calculation
	// P(i), the probability of having i customers in an M/G/infinity system.
	
	// This is not trivial to do.
	
	// Since many of the numbers involved in this calculation are very large or
	// very small -- we must make sure we can represent the answer without having
	// any precision problems.  For example, with lambda/mu = 3000, part of our
	// equation will be 3000^n where n is the number of customers.  Clearly this
	// number will very quickly go beyond the resolution of a Java double.
	// Similarly, (with lambda/mu = 3000) another part of our equation is e^(-3000),
	// which will also quickly go beyond the resolution of a Java double (because
	// it is too small to represent).  Thus, though the overall product of these
	// numbers may have a reasonable value, the individual parts are in themselves
	// non-representable.  After some amount (!!!) of work, I found a way to
	// calculate this without requiring non-primitive types.  The idea is that
	// we incrementally calculate each of the values and temporarily stop if any
	// gets too large or too small.  This way, our overall answer will always be
	// representable.
	public double Pi(int i)
	{
		double lamOverMu = lambda/mu;
		int ePower = (int) lamOverMu;
		double left = lamOverMu - ePower;
		int eNum = 1;
		int factNum = 1;
		int laMuNum = 1;
		double temp = 1.0;
		boolean done = false;
		while (!done)
		{
			// Divide by e if the number is not too small
			if (temp > 0.00000000001 && eNum <= ePower)
			{
				temp = temp / Math.E;
				eNum++;
			}
			// Multiply by (lambda/mu) if the number is not too large
			if (temp < 1.0E9 && laMuNum <= i)
			{
				temp = temp * lamOverMu;
				laMuNum++;
			}
			// Divide by the next part of the factorial if the number
			// is not too small
			if (temp > 0.00000000001 && factNum <= i)
			{
				temp = temp / factNum;
				factNum++;
			}
			// If the number is very small and there is nothing left
			// that will increase it, simply return 0
			if (temp < 0.00000000001 && laMuNum > i)
			{
				temp = 0.0;
				done = true;
			}
			if (eNum > ePower && factNum > i && laMuNum > i)
				done = true;
		}
		// Complete the non-integral part of the exponent on e
		temp = temp*Math.exp(-left);
		return temp;
	}

}