code: //==============================================================================
  //	
  //	Copyright (c) 2002-
  //	Authors:
  //	* Dave Parker <david.parker@comlab.ox.ac.uk> (University of Oxford, formerly University of Birmingham)
  //	* Vojtech Forejt <forejt@fi.muni.cz> (Masaryk University)
  //	
  //------------------------------------------------------------------------------
  //	
  //	This file is part of PRISM.
  //	
  //	PRISM 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 2 of the License, or
  //	(at your option) any later version.
  //	
  //	PRISM 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 PRISM; if not, write to the Free Software Foundation,
  //	Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  //	
  //==============================================================================
  
  package explicit;
  
  import prism.PrismException;
  
  public final class FoxGlynn
  {
  	// constructor parameters
  	private double underflow, overflow, accuracy;
  	private double q_tmax;
  
  	// returned values
  	private int left, right;
  	private double totalWeight;
  	private double[] weights;
  
  	public FoxGlynn(double qtmax, double uf, double of, double acc) throws PrismException
  	{
  		q_tmax = qtmax;
  		underflow = uf;
  		overflow = of;
  		accuracy = acc;
  		run();
  	}
  
  	public final double[] getWeights()
  	{
  		return weights;
  	}
  
  	public final int getLeftTruncationPoint()
  	{
  		return left;
  	}
  
  	public final int getRightTruncationPoint()
  	{
  		return right;
  	}
  
  	public final double getTotalWeight()
  	{
  		return totalWeight;
  	}
  
  	private final void run() throws PrismException
  	{
  		if (q_tmax == 0.0) {
  			throw new PrismException("Overflow: TA parameter qtmax = time * maxExitRate = 0");
  		}
  
  		if (accuracy < 1e-10) {
  			throw new PrismException("Overflow: Accuracy is smaller than Fox Glynn can handle (must be at least 1e-10)");
  		}
  
  		if (q_tmax < 400)
  		{ //here naive approach should have better performance than Fox Glynn
  			final double expcoef = Math.exp(-q_tmax); //the "e^-lambda" part of p.m.f. of Poisson dist.
  			int k; //denotes that we work with event "k steps occur"
  			double lastval; //(probability that exactly k events occur)/expcoef
  			double accum; //(probability that 0 to k events occur)/expcoef
  			double desval = (1-(accuracy/2.0)) / expcoef; //value that we want to accumulate in accum before we stop
  			java.util.Vector<Double> w = new java.util.Vector<Double>(); //stores weights computed so far.
  			
  			//k=0 is simple
  			lastval = 1;
  			accum = lastval;
  			w.add(lastval * expcoef);
  			
  			//add further steps until you have accumulated enough
  			k = 1;
  			do {
  				// TODO: catch case where lastval gets so small that
  				// accnum never reaches desval due to rounding/floating point precision errors (infinite loop)
  
  				lastval *= q_tmax / k; // invariant: lastval = q_tmax^k / k!
  				accum += lastval;
  				w.add(lastval * expcoef);
  				k++;
  			} while (accum < desval);
  
  			//store all data
  			this.left=0;
  			this.right=k-1;
  			this.weights = new double[k];
  
  			for(int i = 0; i < w.size(); i++)
  			{
  				this.weights[i] = w.get(i);			
  			}
  
  			//we return actual weights, so no reweighting should be done
  			this.totalWeight = 1.0;
  		}
  		else
  		{ //use actual Fox Glynn for q_tmax>400
  			final double factor = 1e+10; //factor from the paper, it has no real explanation there
  			final int m = (int) q_tmax; //mode
  			//run FINDER to get left, right and weight[m]
  			{
  				final double sqrtpi = 1.7724538509055160; //square root of PI
  				final double sqrt2 = 1.4142135623730950; //square root of 2
  				final double sqrtq = Math.sqrt(q_tmax);
  				final double aq = (1.0 + 1.0/q_tmax) * Math.exp(0.0625) * sqrt2; //a_\lambda from the paper			
  				final double bq = (1.0 + 1.0/q_tmax) * Math.exp(0.125/q_tmax); //b_\lambda from the paper
  
  				//use Corollary 1 to find right truncation point
  				final double lower_k_1 = 1.0 / (2.0*sqrt2*q_tmax); //lower bound on k from Corollary 1
  				final double upper_k_1 = sqrtq / (2.0*sqrt2); //upper bound on k from Corollary 1
  				double k;
  
  				//justification for increment is in the paper:
  				//"increase k through the positive integers greater than 3"
  				for(k=lower_k_1; k <= upper_k_1;
  					k=(k==lower_k_1)? k+4 : k+1 )
  				{
  					double dkl = 1.0/(1 - Math.exp(-(2.0/9.0)*(k*sqrt2*sqrtq+1.5))); //d(k,\lambda) from the paper
  					double res = aq*dkl*Math.exp(-k*k/2.0)/(k*sqrt2*sqrtpi); //right hand side of the equation in Corollary 1
  					if (res <= accuracy/2.0)
  					{
  						break;
  					}
  				}
  
  				if (k>upper_k_1)
  					k=upper_k_1;
  
  				this.right = (int) Math.ceil(m+k*sqrt2*sqrtq + 1.5); 
  
  				//use Corollary 2 to find left truncation point
  				//NOTE: the original implementation used some upper bound on k,
  				//      however, I didn't find it in the paper and I think it is not needed
  				final double lower_k_2 = 1.0/(sqrt2*sqrtq); //lower bound on k from Corollary 2
  
  				double res;
  				k=lower_k_2;
  				do
  				{
  					res = bq*Math.exp(-k*k/2.0)/(k*sqrt2*sqrtpi); //right hand side of the equation in Corollary 2
  					k++;			
  				}
  				while (res > accuracy/2.0);
  				
  				this.left = (int) (m - k*sqrtq - 1.5);
  				
  				//According to the paper, we should check underflow of lower bound.
  				//However, it seems that for no reasonable values this can happen.
  				//And neither the original implementation checked it
  				
  				double wm = overflow / (factor*(this.right - this.left));
  
  				this.weights = new double[this.right-this.left+1];
  				this.weights[m-this.left] = wm;
  			}
  			//end of FINDER
  
  			//compute weights
  			//(at this point this.left, this.right and this.weight[m] is known)
  			
  			//Down from m
  			for(int j=m; j>this.left; j--)
  				this.weights[j-1-this.left] = (j/q_tmax)*this.weights[j-this.left];
  			//Up from m
  			for(int j=m; j<this.right; j++)
  				this.weights[j+1-this.left] = (q_tmax/(j+1))*this.weights[j-this.left];
  
  			//Compute totalWeight (i.e. W in the paper)
  			//instead of summing from left to right, start from smallest
  			//and go to highest weights to prevent roundoff
  			this.totalWeight = 0.0;
  			int s = this.left;
  			int t = this.right;
  			while (s<t)
  			{
  				if(this.weights[s - this.left] <= this.weights[t - this.left])
  				{
  					this.totalWeight += this.weights[s-this.left];
  					s++;
  				}
  				else
  				{
  					this.totalWeight += this.weights[t-this.left];
  					t--;
  				}
  			}
  			this.totalWeight += this.weights[s-this.left];
  		}
  	}
  
  	public static void test()
  	{
  		double[] weights;
  		double totalWeight = 0.0;
  		int left, right;
  
  		FoxGlynn w = null;
  		try {
  			// q = maxDiagRate, time = time parameter (a U<time b)
  			double q = 1, time = 1;
  			w = new FoxGlynn(q * time, 1.0e-300, 1.0e+300, 1.0e-6);
  		} catch (PrismException e) {
  			// ...
  		}
  		weights = w.getWeights();
  		left = w.getLeftTruncationPoint();
  		right = w.getRightTruncationPoint();
  		totalWeight = w.getTotalWeight();
  		w = null;
  	}
  
  }
  