program Gamma;

type data = array [1..3] of real;
var 	results : text;
	alpha, beta, max, min, eps, G1, G2, R, TSO, ESO, Q, FR : real;
	i, numper : integer;
	a,lt,prob: data;
{$H+}

function LOGAMMA (X: real) : double;                                            	                                                                        
        var A1, A2, A3, A4, A5, F, Z : double;
        begin
        A1:=0.918938533204673; A2:=0.000595238095238; A3:=0.000793650793651;
        A4:=0.002777777777778; A5:=0.083333333333333;
        F:=0.0;
        if X < 7.0 then
        	begin
        		F:=X;
        		while X < 7.0 do
        		begin
        			X:=X+1.0; 	if X < 7.0 then F:=F*X;
        		end; 							
        	F:=-ln(F);
        end; 									
        Z:=1.0/(X*X);
        LOGAMMA:=F + ((X - 0.5)*ln(X)) - X + A1 +
        (((((((-A2*Z) + A3)*Z) - A4)*Z) + A5)/X);
        end;

function AREA (y, p: real) : double;                                            

var a, c, eps, err, F, G : double;

begin

        eps:=0.00001; 				                                
        c:=1.0; G:=1.0;
        F:=exp(p*ln(y)-LOGAMMA(p+1)-y);
        a:=p+1.0; c:=c*y/a; G:=G+c; err:=c/G;

while err > eps do
      begin
           a:=a+1.0; c:=c*Y/a; G:=G+c; err:=c/G;
      end;
      AREA:=G*F;
end;

begin 							                	

	assign(results, 'c:result.txt'); 			                
	rewrite(results);
	numper:=3; 						                
	alpha:=2.0; beta:=2.0; 					                
			                                                        			                                                        
lt[1]:=1.0; 	a[1]:=alpha; 		prob[1]:=0.25;
lt[2]:=2.0; 	a[2]:=alpha*2.0;	prob[2]:=0.45;
lt[3]:=3.0; 	a[3]:=alpha*3.0; 	prob[3]:=0.3;

Q:=20; FR:=0.98; 					                        
TSO:=Q*(1.0-FR); 					                                                                                               

max:=a[3]*10.0; min:=0.0; eps:=0.001; 			                        
R:=0.0; ESO:=0.0;

writeln(results,'lead time data:');                                             
writeln(results,' ');
writeln(results,'      duration    prob.    avg.demand');

for i:=1 to numper do
	begin
		R:=R+((a[i]/beta)*prob[i]);			                

writeln(results, lt[i]:12:2, prob[i]:10:2, (a[i]/beta):12:2);

end; 								                

if R > TSO then 						                
	begin
		while abs (ESO-TSO) > eps do 			                

	begin
		ESO:=0.0;
		for i:=1 to numper do 			                        
			begin
				G2:=1.0 - AREA(beta*R,a[i]);                    
				G1:=1.0 - AREA(beta*R,a[i]+1);
				ESO:=ESO + prob[i]*(((a[i]/beta)*G1)-(R*G2)); 	
			end; 							

		if abs(ESO-TSO) > eps then
			if ESO > TSO then begin 		                
				min:=R;
				R:=(R+max)/2.0;
			end 					                
			else begin 				                
				max:=R;
				R:=(R+min)/2.0;
			end; 					                
		end; 						                
	end; 							                

writeln(results,'******************************************');
writeln(results,' ');
writeln(results,'Parameters of demand distribution',
		' alpha =',alpha:8:3,', beta =',beta:8:3);
writeln(results,' ');
writeln(results,'Optimal reorder point R =',R:8:2);
writeln(results,'for',FR*100.0:8:2,' % fill rate')

end.