> restart;
> with(linalg): Digits:=20:
Warning, new definition for norm
Warning, new definition for trace
> remez:=proc(a,b,N,z,x,f,m)
> # a, b Intervallgrenzen
> # N Grad des gesuchten Polynoms bester Approximation
> # z ist die Folge der Stutzstellen, Anzahl M>>N+2
> # x ist die Startreferenz mit N+2 Stuetzstellen aus z
> # f ist die zu approximierende Funktion vorgegeben an den Stuetzstellen z[k] 
> # m Anzahl der Austauschschritte
> global p;
> local  Newtoninterpol,i,k,j,l,ind,f1,h1,h, S, hi, d, delta, M,sig, indi;
> # option trace;
> # Prozedur zur Berechnung des Interpolationspolynomes
> Newtoninterpol:=proc(x,f,z)
> # f ist der Vektor der Stuetzwerte
> # x ist der Vektor der Stuetzstellen
> # z ist der Vektor der Stuetzstellen, fuer den ein Interpolationswert gesucht ist
> local Newtonkoeff,y,p,m,n,j,k;
> Newtonkoeff:=proc(x,f)
> #x ist der Vektor der Stuetzstellen
> #f ist der Vektor der Stuetzwerte
> local y,n,i,k;
> n:=vectdim(x)-1;
> y:=vector(n);
> for i from 1 to n do y[i]:=f[i] od;
> for k from 1 to n-1 do 
>    for i from n by -1 to k+1 do
>       y[i]:=evalf((y[i]-y[i-1])/(x[i]-x[i-k]))
>    od
> od;
> evalm(y)
> end;
> # Berechnung der Funktionswerte des Interpolationspolynoms
> y:=Newtonkoeff(x,f);
> n:=vectdim(y);
> m:=vectdim(z);
> p:=vector(m);
> for j from 1 to m do
>     p[j]:=y[n];
>     for k from n-1 by (-1) to 1 do p[j]:=evalf(y[k]+(z[j]-x[k])*p[j]) od
> od;
> evalm(p)
> end:
> #Berechnung der f(x[k]) und der Indizes der Knoten x[l] bezueglich z
> M:=vectdim(z);
> ind:=vector(N+2);
> for l from 1 to N+2 do 
> for j from 1 to M do if abs(x[l]-z[j])<10^(-5) then ind[l]:=j; f1[l]:=f[j] fi od od;
> for i from 1 to m do  # Zaehler der Austauschschritte
>  # Berechnung von h
>   h:=0; S:=0; hi:=-1;  # hi Hilfsgroesse, Vorzeichen
>     for k from 1 to N+2 do
>       # Berechnung der d_k
>       d:=1; for j from 1 to k-1 do d:=d*(x[k]-x[j]) od;
>       for j from k+1 to N+2 do d:=d*(x[j]-x[k]) od;
>       d:=1/abs(d);
>       h:=h+hi*d*f1[k]; S:=S+d; hi:=-hi
>     od;
>   h:=h/S;
>   for j from 1 to N+1 do f1[j]:=f1[j]+(-1)^(j+1)*h od;
>   p:=Newtoninterpol(x,f1,z);
>   # Berechnung der Fehlerfunktion delta
>   delta:=vector(M);
>   for j from 1 to M do delta[j]:=(p[j]-f[j]) od;
>   # Berechnung der neuen Referenz x
>   if h=0 then 
>     for j from 1 to ind[2]-1 do
>       h1:=0;
>       if (abs(delta[j])> h1) then sig:=signum(delta[j]); x[1]:=z[j]; h1:=abs(delta[j]); ind[1]:=j fi
>     od  
>   else
>     sig:=signum(delta[ind[1]]); indi:=ind[1]; h1:=abs(h); 
>     for j from 1 to  ind[1]-1 do
>       if (abs(delta[j])>h1) then indi:=j; h1:=abs(delta[j]) fi;
>     od;
>     if signum(delta[indi]) <> sig  then sig:=-sig; ind[1]:=indi
>     else
>      h1:=abs(h);
>      for j from 1 to ind[2]-1 do
>        if (abs(delta[j])> h1) then 
>          if (signum(delta[j])=sig) then x[1]:=z[j]; h1:=abs(delta[j]); ind[1]:=j fi
>        fi;
>      od
>     fi
>   fi;
>  for k from 2 to N+1 do
>   sig:=-sig; h1:=abs(h);
>   for j from ind[k-1]+1 to ind[k+1]-1 do 
>     if (abs(delta[j])> h1) then 
>      if (signum(delta[j])=sig) then x[k]:=z[j]; h1:=abs(delta[j]); ind[k]:=j fi
>     fi;
>   od;
>  od;
>  sig:=-sig; h1:=abs(h);
>  for j from ind[N+1]+1 to M do
>     if (abs(delta[j])> h1) then 
>      if (signum(delta[j])=sig) then x[N+2]:=z[j]; h1:=abs(delta[j]); ind[N+2]:=j fi
>     fi;
>  od;
>  # Neuberechnung von f1 fuer neue Referenz
>  for l from 1 to N+2 do
>   f1[l]:=f[ind[l]] 
>  od 
> od;
> # Berechnung von h fuer neue Referenz
>   h:=0; S:=0; hi:=-1;  # hi Hilfsgroesse, Vorzeichen
>   for k from 1 to N+2 do
>     # Berechnung der d_k
>     d:=1; for j from 1 to k-1 do d:=d/(x[k]-x[j]) od;
>     for j from k+1 to N+2 do d:=d/(x[j]-x[k]) od;
>     h:=h+hi*d*f1[k]; S:=S+d; hi:=-hi
>   od;
>  h:=h/S;
> for j from 1 to N+1 do f1[j]:=f1[j]+(-1)^(j+1)*h od;
> p:=Newtoninterpol(x,f1,z);
> evalm(p)
> end:
> a:=-5: b:=5: N:=7: M:=159: # Beachte dass M+1 durch N+1 teilbar ist!
> r:=round((M+1)/(N+1)):
> z:=vector(M+2): x:=vector(N+2):
> for k from 0 to M+1 do z[k+1]:=evalf((a+b)/2 + (a-b)*cos(k*Pi/(M+1))/2) od:
> x[1]:=z[1]:
> for k from 1 to N+1 do x[k+1]:=z[1+r*k] od:
> print(`Startreferenz`,x);
> f:=vector(M+2):
> for k from 1 to M+2 do f[k]:= 1/(25+z[k]^2)od:
> m:=4:
> remez(a,b,N,z,x,f,m):

Startreferenz, [-5., -4.6193976625564337808, -3.5355339059327376220,

    -1.9134171618254488587, 0, 1.9134171618254488587,

    3.5355339059327376220, 4.6193976625564337808, 5.]

> print(`Endreferenz`,x);

Endreferenz, [-5., -4.5809397855856795542, -3.3940037276647087069,

    -1.7305852853874648824, 0, 1.7305852853874648824,

    3.3940037276647087069, 4.5809397855856795542, 5.]

> with(plots):
> F:=plot(1/(25+x^2),x=-5..5):
> G:=listplot([seq([z[t],p[t]],t=1..M+2)],thickness=0):
> display({F,G});
> listplot([seq([z[t],p[t]-1/(25+z[t]^2)],t=1..M+2)],thickness=0);
> 


>