ABSTRACT
B. sndremez Exchanges at most n+1 numbers with numbers out of a reference set (see [Me64]). Given a sequence of real numbers gi, i=0,...,m, which contains a subsequence gs(j), j=0,1,...,n<m; s(j) > s(j-1), for which a) gs(j) = G, j=0,...,n, and b) sign gs(j) = -sign gs(j-1), j=1,...,n, a sequence of integers S(j), j=0,1,...,n; S(j) > S(j-1); S(0) ≥ 0, S(n) ≤ m, is determined such that A) gS(j) ≥ G, j=0,...,n, B) sign gS(j) = -sign gS(j-1), j=1,...,n, and C) gS(j) > G for at least one j in the range 0 ≤ j ≤ n if gi > G for at least one i in the range 0 ≤ i ≤ n. The process used is as follows: let H(1) = maxgi for s(0) < i < s(1) and gI(1) = H(1) (H(1) is zero if s(0) = s(1) and I(1) is then undefined; a similar convention is adopted below). If H(1) ≤ G then S(0) = s(0), S(1) = s(1); otherwise if sign gI(1) = sign gs(0) then S(0) = I(1), S(1)=s(1), whereas if sign gI(1) ≠ sign gs(0), then S(0) = s(0), S(1) = I(1). Let H(2) = maxgi for s(1) < i < s(2) and gI(2) = H(2). If H(2) ≤ G then S(2) = s(2); otherwise if sign gI(2) = sign gS(1) then S(1) is replaced by I(2) and S(2) = s(2), whereas if sign gI(2) ≠ sign gS(1) then S(2) = I(2). The numbers gi, s(2)<i<s(3), are inspected, S(2) is possibly replaced, and S(3) is determined as before. In this way S(0),...,S(n) are determined. Let H' = maxgi for 0 ≤ i < s(0) with gI'=H', and H" = maxgi for s(n) < i ≤ m, with gI" = H". If H', H" < G, the numbers S(j) are left as they stand. If not, then if H' > H" then a) if sign gI' = sign gS(0) then a') if H'>gs(0), S(0) is replaced by I' and a") if sign gI" = gs(n) and H" > gs(n), S(n) is replaced by I"; b) if sign gI' ≠ sign gS(0) then b') if H' > gs(n), S(k) is replaced by S(k-1), k=n,...,1, and
S(0) is replaced by I' and b") (with the new value of S(n)) if sign gI" = sign gS(n) and H" > gs(n), S(n) is replaced by I". If H" ≥ H' then similar modifications take place. sndremez is an auxiliary procedure used in minmaxpol. Procedure parameters: void sndremez (n,m,s,g,em) n,m: int; entry: the number of points to be exchanged is smaller than or equal to
n+1; the reference set contains the numbers 0,1,...,m, (m ≥ n); s: int s[0:n]; entry: in s one must give n+1 (strictly) monotone increasing numbers out
of 0,1,...,m; exit: n+1 (strictly) monotone increasing numbers out of the numbers
0,1,...,m; g: double g[0:m]; entry: in array g[0:m] one must give the function values; em: double em[0:1]; entry: 0 < em[0] ≤ g[i], i=0,...,m; exit: em[1] = infinity norm of array g[0:m]. Procedure used: infnrmvec. public static void sndremez(int n, int m, int s[], double g[], double em[]) { int s0,sn,sjp1,i,j,k,up,low,nm1; double max,msjp1,hi,hj,he,abse,h,temp1,temp2; int itmp[] = new int[1]; int jtmp[] = new int[1];
s0=sjp1=s[0]; he=em[0]; low=s0+1; max=msjp1=abse=Math.abs(he); nm1=n-1; for (j=0; j<=nm1; j++) { up=s[j+1]-1; h=Basic.infnrmvec(low,up,itmp,g); i=itmp[0]; if (h > max) max=h; if (h > abse) if (he*g[i] > 0.0) { s[j] = (msjp1 < h) ? i : sjp1; sjp1=s[j+1]; msjp1=abse; } else { s[j]=sjp1; sjp1=i; msjp1=h; }
else { s[j]=sjp1; sjp1=s[j+1]; msjp1=abse; } he = -he; low=up+2; } sn=s[n]; s[n]=sjp1; hi=Basic.infnrmvec(0,s0-1,itmp,g); i=itmp[0]; hj=Basic.infnrmvec(sn+1,m,jtmp,g); j=jtmp[0]; if (j > m) j=m; if (hi > hj) { if (hi > max) max=hi; temp1 = (g[i] == 0.0) ? 0.0 : ((g[i] > 0.0) ? 1.0 : -1.0); temp2 = (g[s[0]]==0.0) ? 0.0 : ((g[s[0]]>0.0) ? 1.0 : -1.0); if (temp1 == temp2) { if (hi > Math.abs(g[s[0]])) { s[0]=i; if (g[j]/g[s[n]] > 1.0) s[n]=j; } } else { if (hi > Math.abs(g[s[n]])) { s[n] = (g[j]/g[s[nm1]] > 1.0) ? j : s[nm1]; for (k=nm1; k>=1; k--) s[k]=s[k-1]; s[0]=i; } } } else { if (hj > max) max=hj; temp1 = (g[j] == 0.0) ? 0.0 : ((g[j] > 0.0) ? 1.0 : -1.0); temp2 = (g[s[n]]==0.0) ? 0.0 : ((g[s[n]]>0.0) ? 1.0 : -1.0); if (temp1 == temp2) { if (hj > Math.abs(g[s[n]])) { s[n]=j; if (g[i]/g[s[0]] > 1.0) s[0]=i; } } else if (hj > Math.abs(g[s[0]])) { s[0] = (g[i]/g[s[1]] > 1.0) ? i : s[1]; for (k=1; k<=nm1; k++) s[k]=s[k+1]; s[n]=j; } } em[1]=max; }
C. minmaxpol Calculates the coefficients of the polynomial (as a sum of powers) which approximates a function, given for discrete arguments, in such a way that the infinity norm of the error vector is minimized. With yj (j=0,...,m) being a sequence of real arguments (yj-1 < yj; j=1,...,m) and f(yj) corresponding function values, minmaxpol determines the coefficients {ci} of that polynomial
∑ −
=
n
i
i i yc
for which
∑ −
=
− 1
0 )(max
n
i
i jij ycyf (1)
is a minimum as j ranges over the values 0,...,m (m > n). The method used [Me64] involves the iterative construction of polynomials
∑ −
=
0 ,
n
i
i ik yc
for which
( )njycyf kj n
i
i jksikjks ,,0)1()(
0 ),(,),( K=−=−∑
−
=
λ (2)
where s(k,j), j=0,...,n, is a strictly increasing sequence selected from 0,1,...,m. The discrepancies
∑ −
=
0 ,, )(
n
i
i jikjjk ycyfg
at all points j=0,...,m are then inspected, and from these a sequence gk+1,j = gk,s(k+1,j)', j=0,...,n, is constructed which a) possesses the alternating property sign gk+1,j = -sign gk+1,j-1, j=1,...,n, and b) for which gk+1,j ≥ gk,j, j=0,...,n, with (unless the process is to be terminated) gk+1,j>gk,j for at least one j (for the details of this construction, see the documentation to sndremez). The coefficients {ck+1,i} in the succeeding polynomial are determined (by use of Newton's interpolation formula, and subsequent reduction to the given polynomial form) from condition (2) with k replaced by k+1. Initially s(0,j), j=0,...,n, is a sequence of integer approximations to the points at which Tn((2x-m)/m), where Tn(y) = cos{n arc cos(y)} being a Chebyshev polynomial (see the documentation of ini) assumes its maximal absolute values in the range 0 ≤ x ≤ m. The procedure ini is used to set the s(0,j). Procedure parameters: void minmaxpol (n,m,y,fy,co,em) n: int; entry: the degree of the approximating polynomial; n ≥ 0; m: int;
entry: the number of reference function values is m+1; y,fy: double y[0:m], fy[0:m]; entry: fy[i] is the function values at y[i], i=0,...,m; co: double co[0:n]; exit: the coefficients of the approximating polynomial (co[i] is the
coefficient of yi); em: double em[0:3]; entry: em[2]: maximum allowed number of iterations, say 10*n+5; exit: em[0]: the difference of the given function and the polynomial in
the first approximation point; em[1]: the infinity norm of the error of approximation over the
discrete interval; em[3]: the number of iterations performed. Procedures used: elmvec, dupvec, newton, pol, newgrn, ini, sndremez. public static void minmaxpol(int n, int m, double y[], double fy[], double co[], double em[]) { int np1,k,pomk,count,cnt,j,mi,sjm1,sj,s0,up; double e,abse,abseh; int s[] = new int[n+2]; double x[] = new double[n+2]; double b[] = new double[n+2]; double coef[] = new double[n+2]; double g[] = new double[m+1];
sj=0; np1=n+1; ini(np1,m,s); mi=(int)em[2]; abse=0.0; count=1; do { pomk=1; for (k=0; k<=np1; k++) { x[k]=y[s[k]]; coef[k]=fy[s[k]]; b[k]=pomk; pomk = -pomk; } newton(np1,x,coef); newton(np1,x,b); em[0]=e=coef[np1]/b[np1]; Basic.elmvec(0,n,0,coef,b,-e); Algebraic_eval.newgrn(n,x,coef); s0=sjm1=s[0]; g[s0]=e; for (j=1; j<=np1; j++) { sj=s[j];
up=sj-1; for (k=sjm1+1; k<=up; k++) g[k]=fy[k]-Algebraic_eval.pol(n,y[k],coef); g[sj] = e = -e; sjm1=sj; } for (k=s0-1; k>=0; k--) g[k]=fy[k]-Algebraic_eval.pol(n,y[k],coef); for (k=sj+1; k<=m; k++) g[k]=fy[k]-Algebraic_eval.pol(n,y[k],coef); sndremez(np1,m,s,g,em); abseh=abse; abse=Math.abs(e); cnt=count;
count++;
} while (count <= mi && abse > abseh); em[2]=mi; em[3]=cnt; Basic.dupvec(0,n,0,co,coef); }
