1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
comment
Remez algorithm
endcomment
.. ### Remez algorithm ###
function remezhelp (x,y,n)
## {t,d,h}=remezhelp(x,y,deg) is the case, where x are exactly deg+2
## points.
d1=interp(x,y); d2=interp(x,mod(1:n+2,2)*2-1);
h=d1[n+2]/d2[n+2];
d=d1-d2*h;
return {x[1:n+1],d[1:n+1],h};
endfunction
remezeps=0.00001;
function remez(x,y,deg,tracing=0)
## {t,d,h,r}=remez(x,y,deg) computes the divided difference form
## of the polynomial best approximation to y at x of degree deg.
## remez(x,y,deg,1) does the same thing with tracing.
## h is the discrete error (with sign), and r the rightmost
## alternation point, which is missing in t. To evaluate the
## polynomial in v use interpval(t,d,v).
global remezeps;
n=cols(x);
ind=round(linspace(1,n,deg+1),0); h=0;
repeat
{t,d,hnew}=remezhelp(x[ind],y[ind],deg);
r=y-interpval(t,d,x); hh=hnew;
if tracing;
plot(x,r); xgrid(x[ind],ticks=0,color=9);
ygrid([-hnew,0,hnew],ticks=0);
color(1);
title("Weiter mit Taste"); wait(180);
endif;
indnew=ind;
e=extrema(r[1:ind[2]]);
if (hh<0);
indnew[1]=e[2];
else
indnew[1]=e[4];
endif;
hh1=hh;
loop 2 to deg+1;
e=extrema(r[ind[#-1]:ind[#+1]]);
if (hh>0);
indnew[#]=e[2]+ind[#-1]-1;
else
indnew[#]=e[4]+ind[#-1]-1;
endif;
hh=-hh;
end;
e=extrema(r[ind[deg+1]:cols(r)]);
if (hh>0);
indnew[deg+2]=e[2]+ind[deg+1]-1;
else
indnew[deg+2]=e[4]+ind[deg+1]-1;
endif;
hh=-hh;
mm=max(abs(r[indnew]));
if (indnew[deg+2]<cols(r));
e=extrema(r[ind[deg+2]:cols(r)]);
if (hh>0)
if (abs(r[e[1]])>mm);
indnew=indnew[2:deg+2]|(e[2]+ind[deg+2]-1);
endif;
else
if (abs(r[e[3]])>mm);
indnew=indnew[2:deg+2]|(e[4]+ind[deg+2]-1);
endif;
endif;
endif;
mm=max(abs(r[indnew]));
if (indnew[1]>1);
e=extrema(r[1:ind[1]]);
if (hh1>0)
if (abs(r[e[1]])>mm); indnew=e[2]|indnew[1:deg+1]; endif;
else
if (abs(r[e[3]])>mm); indnew=e[4]|indnew[1:deg+1]; endif;
endif;
endif;
ind=indnew;
if (abs(hnew)<(1+remezeps)*abs(h)) break; endif;
h=hnew;
end;
{t,d,h}=remezhelp(x[ind],y[ind],deg);
if tracing;
xplot(x,y-interpval(t,d,x));
color(15);
xgrid(x[ind],ticks=0,color=9);
ygrid([-hnew,0,hnew],ticks=0);
color(1);
title("Weiter mit Taste"); wait(180);
endif;
return {t,d,h,x[ind[deg+2]]};
endfunction
|
