% BERNSTEIN Sharpness of Bernstein's inequality for
%   Jacobi polynomials P_n(a,b,;.) with b>=a>=-1/2.
%   The output is rho=c_n*||F_n||.
%
function rho=bernstein(n,a,b,q)
if a<-1/2 | b<-1/2 | b<a
  disp('parameters a and/or b not in range')
  return
end
tol=1e2*eps;
pnum=1; pden=1; p2=1;
for nu=1:n
  pnum=(1+(n+a+b)/nu)*pnum;
  pden=(1+q/nu)*pden;
  p2=(1-1/(2*nu))*p2;
end 
c0=pnum/2^(n+(a+b+1)/2);
c1=(n+(a+b+1)/2)^(q+1/2)/(gamma(1+q)*pden);
c2=sqrt(pi)*c1*p2;
c=sqrt(pi)*c0*c1;
extr=zeros(n+1,1);
ab=r_jacobi(n,a,b);
xw=gauss(n,ab);
x=xw(:,1);
x1=[-1 x' 1]';
k0=1; k1=n+1;
if a==-1/2
  if b>-1/2
    k1=n;
  else
    rho=1;
    return
  end
end
for k=k0:k1
  t0=0; t1=(x1(k)+x1(k+1))/2;
  while abs(t1-t0)>tol
    t0=t1;
    t1=t0-fbern(t0,a,b,x)/f1bern(t0,a,b,x);
  end
  p=prod(t1-x);
  extr(k)=(1-t1)^(a/2+1/4)*(1+t1)^(b/2+1/4)*abs(p);
end 
rho=c*max(extr);
if a==-1/2
  if c2>rho
    rho=c2;
  end
end