function coeff = remezp(nomf,wf,a,b,n,iter,np);

% remezp.m

% Meilleure approximation de degre <= n au sens de Tchebycheff

%  executer  remezf.m  dans matlab

% d'apres
%                       the COCA - package.
% Copyright (C) 1992, Bernd Fischer and Jan Modersitzki.
% All rights reserved.
%
disp(nomf);disp(wf);
% points de grille
xg=(a+b)*ones(1,np+1)/2+ (b-a)*cos(0:pi/np:pi+0.00001)/2;
x=xg;fg=eval(nomf);wg=eval(wf);clear x;
axis([a,b,-norm(wg.*fg,inf),norm(wg.*fg,inf)]);plot(xg,fg.*wg);%meta remezpg;

% points de depart
xi=(a+b)*ones(1,n+2)/2+ (b-a)*cos(0:pi/(n+1):pi+0.00001)/2;
signes(1)=1;for i=2:n+2;signes(i)=-signes(i-1);end;

for it=1:iter

disp(' iter '); disp(it); disp(xi);
 % interpolation par pol orthogonaux discrets (tres stable)
  polorth(1,:)=ones(1,n+2) ; normep(1)=n+2;  % degre 0
  for ip=1:n+1;
   aux=xi .* polorth(ip,:) .* polorth(ip,:); alpha(ip)=sum(aux)/normep(ip);
   polorth(ip+1,:)=(xi - alpha(ip)*ones(1,n+2) ) .* polorth(ip,:);
   if ip>1
     beta2(ip)=normep(ip)/normep(ip-1);
     polorth(ip+1,:)=polorth(ip+1,:) - beta2(ip)*polorth(ip-1,:);
   end;
   aux=polorth(ip+1,:) .* polorth(ip+1,:) ; normep(ip+1)=sum(aux);
  end;

 % interpolation de f
  x=xi;fi = eval(nomf);wi=eval(wf); clear x;
                               interpf = (polorth*fi') ./ normep' ;
 % interpolation de (-1)^i/w(xi)
                               interps = (polorth* (signes./wi)' )./normep';
  ee = interpf(n+2)/interps(n+2);
% disp(interpf');
% disp(interps');
  interpf = interpf - ee*interps ;

 % valeurs sur la grille
  polg(1,:)=interpf(n+1)*ones(1,np+1);
  if n>0
  polg(2,:)=xg .* polg(1,:) -alpha(n)*polg(1,:) +interpf(n)*ones(1,np+1);
  for ip=n-1:-1:1
   polg(3,:)=xg .* polg(2,:) -alpha(ip)*polg(2,:) -beta2(ip+1)*polg(1,:) + ...
                                             interpf(ip)*ones(1,np+1);
   polg(1,:)=polg(2,:); polg(2,:)=polg(3,:);
  end;
  polg(1,:)=polg(2,:);
  end;
  aux=wg.*(fg-polg(1,:));
  indp=0; E=0 ;
  for ig=1:np+1;
   if indp < n+2
   if xg(ig) <= xi(indp+1)
      indp=indp+1;
   end;
   end;
   if abs(aux(ig))>E
      E=abs(aux(ig));
      xgi=xg(ig);
      signext=sign(aux(ig));
      indext=indp;
   end;
  end;
  disp([ee E signext indext]);
  auxv=zeros(1,np+1);
  indp=0;
  for ig=1:np+1;
   if indp<n+2
   if xg(ig)<= xi(indp+1)
    auxv(ig)=aux(ig);
    indp=indp+1;
   end;
   end;
  end;
  axis([a,b,-E,E]);plot(xg,aux,xg,ee*ones(1,np+1),xg,-ee*ones(1,np+1), ...
     xg,auxv);
  titre=sprintf( ' %3.0f   %8.4f   %8.4f ',it,abs(ee),E);title(titre);
  for i=1:n+2;text(xi(i),ee*signes(i)-0.1,'*');end;
% echange
  if indext == 0
   if signext ~= sign(ee)
    for i=n+2:-1:2;xi(i)=xi(i-1);end;
   end;
   xi(1)=xgi;
  elseif indext == n+2
   if signext ~= signes(n+2)* sign(ee)
    for i=1:n+1;xi(i)=xi(i+1);end;
   end;
   xi(n+2)=xgi;
  else
   if signext == signes(indext)*sign(ee)
                       xi(indext)=xgi
   end;
   if signext == signes(indext+1)*sign(ee)
                       xi(indext+1)=xgi
   end;
  end;


% meta remezpg ; % faire gpp remezpg.met -dps
print -dps remezpg
end ;  % iterations

% polynome dans la base x^n   x^(n-1)  ...  1

format long;
disp('             p                 alpha                      beta2 ');
disp([interpf(1:n+1) alpha' beta2']);
polm0=zeros(1,n+1);polm0(n+1)=interpf(n+1);
if n>0
polm1=zeros(1,n+1);polm1(n)=polm0(n+1);
polm1(n+1)=-alpha(n)*polm0(n+1)+interpf(n);polm2=zeros(1,n+1);
for ip=n-1:-1:1
 polm2(ip:n)=polm1(ip+1:n+1);polm2(n+1)=0;
 polm2=polm2-alpha(ip)*polm1-beta2(ip+1)*polm0;polm2(n+1)=polm2(n+1)+interpf(ip);
%disp([polm0' polm1' polm2']);
 polm0=polm1;polm1=polm2;
end;
polm0=polm1;
end;
disp(polm0');

coeff=polm0;

return;
