[betaresult,a,b]=beta(a,b);

function [betaresult,a,b]=beta(a,b);
betaresult=[];
persistent alnsml beta firstCall xmax xmin ; if isempty(firstCall),firstCall=1;end; 

if isempty(alnsml), alnsml=0; end;
if isempty(betaresult), betaresult=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
%***BEGIN PROLOGUE  BETA
%***PURPOSE  Compute the complete Beta function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C7B
%***TYPE      SINGLE PRECISION (BETA-S, DBETA-D, CBETA-C)
%***KEYWORDS  COMPLETE BETA FUNCTION, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% BETA computes the complete beta function.
%
% Input Parameters:
%       A   real and positive
%       B   real and positive
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  ALBETA, GAMLIM, GAMMA, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   900727  Added EXTERNAL statement.  (WRB)
%***end PROLOGUE  BETA
if firstCall,   xmax =[0.];  end;
if firstCall,  alnsml=[0.];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BETA
if( alnsml==0.0 )
[xmin,xmax]=gamlim(xmin,xmax);
alnsml = log(r1mach(1));
end;
%
if( a<=0. || b<=0. )
xermsg('SLATEC','BETA','BOTH ARGUMENTS MUST BE GT 0',2,2);
end;
%
if( a+b<xmax )
betaresult = gamma(a).*gamma(b)./gamma(a+b);
end;
if( a+b<xmax )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
[betaresult ,a,b]=albeta(a,b);
if( betaresult<alnsml )
xermsg('SLATEC','BETA','A AND/OR B SO BIG BETA UNDERFLOWS',1,2);
end;
%
betaresult = exp(betaresult);
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK BETAI