Example 1:
> P = fpminimax(cos(x),6,[|DD, DD, D...|],[-1b-5;1b-5]);
> printexpansion(P);
(0x3ff0000000000000 + 0xbc09fda15e029b00) + x * ((0x3af9eb57163024a8 + 0x37942c2f3f3e3839) + x * (0xbfdfffffffffff98 + x * (0xbbd1693f9c028849 + x * (0x3fa5555555145337 + x * (0x3c7a25f610ad9ebc + x * 0xbf56c138142da5b0)))))
Example 2:
> P = fpminimax(sin(x),6,[|32...|],[-1b-5;1b-5], fixed, absolute);
> display = powers!;
> P;
x * (1 + x^2 * (-357913941 * 2^(-31) + x^2 * (35789873 * 2^(-32))))
Example 3:
> P = fpminimax(exp(x), [|3,4|], [|D,24|], [-1/256; 1/246], 1+x+x^2/2);
> display = powers!;
> P;
1 + x * (1 + x * (1 * 2^(-1) + x * (375300225001191 * 2^(-51) + x * (5592621 * 2^(-27)))))
Example 4:
> print("Weighted approximation: look for p such that (p-f)*w is small");
Weighted approximation: look for p such that (p-f)*w is small
> n = 3;
> f = exp(x);
> w = sqrt(1-x^2);
> monomials = [||];
> for i from 0 to n do { monomials = monomials :. x^i*sqrt(1-x^2); };
> q = fpminimax(f*w, monomials, [|D...|], [-99/100,99/100],absolute);
> display = dyadic!;
> print(q);
3194834406077049b-54 * x^3 * sqrt(1 - x^2) + 4792048811386295b-53 * x^2 * sqrt(1 - x^2) + 8982938415377717b-53 * x * sqrt(1 - x^2) + 561436817847349b-49 * sqrt(1 - x^2)
> leadingcoeff = proc(expr) {
var res;
match(expr) with
a*b+c : {res = [|a,c|]; }
a*b : {res = [|a,0|]; }
default: {res = error; };
return res;
};
> p = 0;
> for i from n to 0 by -1 do {
L = leadingcoeff(q);
p = p + L[0]*x^i;
q = L[1];
};
> print(p);
3194834406077049b-54 * x^3 + 4792048811386295b-53 * x^2 + 8982938415377717b-53 * x + 561436817847349b-49 * x^0
> display=decimal!;
> dirtyinfnorm((p-f)*w, [-1,1]);
2.7416001717508988167160223812127007458632536307517e-3
> r = remez(f*w,n,[-99/100,99/100],w);
> dirtyinfnorm((r-f)*w, [-1,1]);
2.7416001717508734413506680717675258974378965304354e-3
Example 5:
> I = [-1b-20, 1b-22];
> P0 = fpminimax(exp(x), 4, [|DD,DD,DD,D,SG|], I);
> P1 = fpminimax(exp(x), 4, [|D,D,D,D,SG|], I, 1+x+x^2/2);
> for k from 0 to 2 do { DD(coeff(P0,k)) == coeff(P0,k); };
true
true
true
> for k from 0 to 2 do { DD(coeff(P1,k)) == coeff(P1,k); };
true
true
true
> dirtyinfnorm(P0/exp(x)-1, I);
6.2744565211592155928616387674345117702132844508175e-35
> dirtyinfnorm(P1/exp(x)-1, I);
5.8284777124834959383794298527146965060524845952998e-35
Example 6:
> f = cos(exp(x));
> pstar = remez(f, 5, [-1b-7;1b-7]);
> listpoints = dirtyfindzeros(f-pstar, [-1b-7; 1b-7]);
> P1 = fpminimax(f, 5, [|DD...|], listpoints, absolute, default, default, pstar);
> P2 = fpminimax(f, 5, [|D...|], listpoints, absolute, default, default, pstar);
> P3 = fpminimax(f, 5, [|D, D, D, 24...|], listpoints, absolute, default, default, pstar);
> print("Error of pstar: ", dirtyinfnorm(f-pstar, [-1b-7; 1b-7]));
Error of pstar: 7.9048441259903026332577436001060063099817726177425e-16
> print("Error of P1: ", dirtyinfnorm(f-P1, [-1b-7; 1b-7]));
Error of P1: 7.9048441259903026580081299123420463921479618202064e-16
> print("Error of P2: ", dirtyinfnorm(f-P2, [-1b-7; 1b-7]));
Error of P2: 8.2477144579950871737109573536791331686347620955985e-16
> print("Error of P3: ", dirtyinfnorm(f-P3, [-1b-7; 1b-7]));
Error of P3: 1.08454277156993282593701156841863009789063333951055e-15
Example 7:
> L = [|exp(x), sin(x), cos(x)-1, sin(x^3)|];
> g = (2^x-1)/x;
> p = fpminimax(g, L, [|D...|], [-1/16;1/16],absolute);
> display = powers!;
> p;
-3267884797436153 * 2^(-54) * sin(x^3) + 5247089102535885 * 2^(-53) * (cos(x) - 1) + -8159095033730771 * 2^(-54) * sin(x) + 6243315658446641 * 2^(-53) * exp(x)
Example 8:
> n = 9;
> T = [|1, x|];
> for i from 2 to n do T[i] = canonical(2*x*T[i-1]-T[i-2]);
> g = (2^x-1)/x;
> PCheb = fpminimax(g, T, [|DD,DE...|], [-1/16;1/16],absolute);
> display = dyadic!;
> print(PCheb);
112469560253864905423137861008261b-107 + x * (76130818934358736339243350051b-98 + x * (34355379230514621195759363b-89 + x * (381013348011182609334519067b-95 + x * (825307274780189286345649b-89 + x * (3050983523250669954220137b-94 + x * (1180123116639845252291b-86 + x * (6543992088039485657053b-92 + x * (15750497046710770365b-87 + x * 8733930098894247371b-90))))))))
Example 9:
> d = [46768052394588893382517914646921056628989841b-165;1b-7];
> n = 4; f = 1; w = sqrt(x)/acos(1-x);
> p = remez(f, n, d, w);
> monomials = [||]; Lcoeffs = [||];
> for i from 0 to n do {
monomials = monomials :. x^i*sqrt(x);
Lcoeffs = Lcoeffs :. coeff(p, i);
};
> pf = fpminimax(acos(1-x), monomials, [|D...|], d, relative, floating, 0, Lcoeffs);
> display = dyadic!;
> print(pf);
6237591828287231b-61 * x^4 * sqrt(x) + 1137310727877715b-57 * x^3 * sqrt(x) + 7642862123083659b-58 * x^2 * sqrt(x) + 1061508612083747b-53 * x * sqrt(x) + 6369051672525773b-52 * sqrt(x)