function ab=r_dblexp_freud(N,al,be,ga,a,b)
%   ab=R_DBLEXP(N,al,be,,ga,a,b) generates the first N recurrence
%   coefficients for the monic polynomials orthogonal on [0,Inf]
%   relative to the weight function w(t)=t^ga*exp(-a*t^(-al)-b*t^be),
%   ga>0, a>0, b>0, al>0, be>0, using a discretization procedure
%   based on Gauss-Freud quadrature.

global mc mp iq idelta irout eps0 Mmax ab1 ab2
mc=2; mp=0; iq=1; idelta=2; irout=1;
t0=(a*al/(b*be))^(1/(al+be));
c=b^(1/be)*t0; mu=ga; nu=be; 
%dig0=72; dd=4; nofdig=12;
%[ab,dig]=dig_u1sf(100,c,mu,nu,dig0,dd,nofdig); 
tic
ab=sr_u1sf(80,12,100,c,mu,nu);
toc
ab2=double(ab)
c=a^(1/al)/t0; mu=-(ga+2); nu=al;
%[ab,dig]=dig_u1sf(100,c,mu,nu,dig0,dd,nofdig);
tic
ab=sr_u1sf(80,12,100,c,mu,nu);
toc
ab1=double(ab)
[ab,Mcap,kount]=mcdis(N,eps0,@quaddblexp_freud,Mmax)
