-
-
Save Observatorio-de-Matematica/2b5c5afda8be19e4f1fa6045c2fd587d to your computer and use it in GitHub Desktop.
Example of main stream linear program solvers failing on small problems
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| The attached files contain a piecwise degree 3 polynomial L-inf norm approximation for the sine and cosine function | |
| in the interval 0...1 with 4 and 64 pieces. | |
| The functions are sampled at the extrema of a Chebyshev polynomial of degree 4 since the residuals of an L-inf polynomial | |
| approximation are by construction typically close to Chebyshev polynomials. The small deviations in the extrema positions usually | |
| don't lead to substantial deviations in the extrema values. | |
| This leads to rather small problems with only 5 variables (4 coefficients and one maximum error variable) and 10 inequalities | |
| (5 extrema points, each positive and negative bounded). | |
| The solution is implemented in: | |
| - Maxima / wxMaxima using an exact rational solution | |
| - Python using the GLPK wrapper using the the simplex, interior and exact method | |
| - Matlab 2020b | |
| The obersvations for 4 pieces are as follows: | |
| - All 5 implementations (couting the 3 Python variants as 3) give the same worst case error: | |
| - Maxima: | |
| Max Error Sin = 9.752e-7 (-19.968 bits) | |
| Max Error Cos = 1.261e-6 (-19.597 bits) | |
| - Python: | |
| Max Error Sin Simplex = 9.752e-07 (-19.968 bits) | |
| Max Error Cos Simplex = 1.261e-06 (-19.597 bits) | |
| Max Error Sin Interior = 9.753e-07 (-19.968 bits) | |
| Max Error Cos Interior = 1.261e-06 (-19.597 bits) | |
| Max Error Sin Exact = 9.753e-07 (-19.968 bits) | |
| Max Error Cos Exact = 1.261e-06 (-19.597 bits) | |
| - Matlab: | |
| Max Error Sin = 9.752e-07 (19.968 bits) | |
| Max Error Cos = 1.261e-06 (19.597 bits) | |
| - The first segment of the sine approximation is off in Matlab though: | |
| - Maxima: 1.584083112995679·10^-7 (see the top left matrix element in section 1.10.1 of LInfApprox_Sin_3deg_4seg.wxmx) | |
| - Matlab: 8.99674244103643e-07 (see top element in column 5 of variable ResultSin after running LinfApprox_Sin_3deg_4seg.m) | |
| - For the other 3 segments of sine and all 4 segments of cosine, Matlab and Maxima results are the same. | |
| (in the Maxima file in section 1.10.1, the first column is for sine, the second column is for cosine). | |
| - Note the "WARNING: deviation is not well distributed " in Matlab - which is not from the solver but in my code | |
| and checks if the deviation is all 5 points is similar. | |
| - Note that the solver does not give any error or warning in this case in Matlab. | |
| The obersavtions for 64 pieces are as follows: | |
| - The 5 implementations give different results | |
| - Maxima: | |
| Max Error Sin = 1.625e-11 (-35.841 bits) | |
| Max Error Cos = 1.939e-11 (-35.586 bits) | |
| - Python: | |
| Max Error Sin Simplex = 2.810e-08 (-25.085 bits) | |
| Max Error Cos Simplex = 4.984e-08 (-24.258 bits) | |
| Max Error Sin Interior = 1.335e-10 (-32.803 bits) | |
| Max Error Cos Interior = 1.212e-10 (-32.942 bits) | |
| Max Error Sin Exact = 1.335e-10 (-32.803 bits) | |
| Max Error Cos Exact = 1.212e-10 (-32.942 bits) | |
| - Matlab | |
| Max Error Sin = 9.558e-07 (19.997 bits) | |
| Max Error Cos = 7.199e-07 (20.406 bits) | |
| The Matlab solution is obviosuly bogus since it doesn't achieve a better error with 64 segments | |
| than with 4 segments. Almost all of the segments are worse then the worst case error found by Maxima. | |
| The Python solution - at least for the interior point and exact algorithm - is less of but the error | |
| is still 3 bits larger than the error found with Maxima. | |
| Also note that also the exact algorithm in Python produces "WARNING: deviation is not well distributed" | |
| warnings for almost all segments. | |
| The Python exact algorithm is especially disappointing, because it is supposed to be excat. | |
| Also the Python error is more dangerous than the Matlab error, since the Matlab result is quite | |
| obviosuly wrong. The Python result requires more experience. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ | |
| /* [ Created with wxMaxima version 21.02.0 ] */ | |
| /* [wxMaxima: title start ] | |
| Least L∞ approximation of | |
| Sine and Cosine in 0...1 | |
| 64 pieces with degree 3 | |
| [wxMaxima: title end ] */ | |
| /* [wxMaxima: section start ] | |
| Intro | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| funcs:[sin,cos]; | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Main parameter | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| degree:3$ | |
| n_segments:4$ | |
| x_min:0$ | |
| x_max:1$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Scalings | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| Scaling denominators for x^0 ... x^n | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: comment start ] | |
| Note: giving 6 extra bits for c0 increases precision after rounding to to 37.2 bits | |
| Note: this might be a case where ILP is worthwhile | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| nbits:31$ | |
| max_deriv:1$ | |
| denoms:makelist(2^nbits,i,0,degree); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Libraries | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| load(simplex)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Settings | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| Some corner cases involve 0^0 terms. | |
| 0^0 is usually undefined, but we tell the maxima simplifier to replace it with 1 here. | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| simp : false $ | |
| tellsimp (0^0, 1) $ | |
| tellsimp (0.0^0, 1) $ | |
| simp : true $ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| The simplex algorithm epsilon can be set to 0 for rational input | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| epsilon_lp:0$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| Setting for float to rat | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| ratepsilon:1.0e-18$ | |
| ratprint:false$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Linf approximation of a function sampled at the extrema of the Chebyshev polynomial with given degree | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| We put the samples at the extrema of the Chebeshev poylnomial, because the residual is usually close to a chebychev polynomial | |
| Chebyshev extrama: http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| Chebyshev_Extrema(n,xmin,xmax):=makelist(float((cos(%pi*i/n)-1)/-2*(xmax-xmin)+xmin),i,0,n); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approximation_poly_linf_cheb(f,degree,xmin,xmax,eps_y,denom):=block( | |
| [ratepsilon,xp,yp,ydiff,gep,gem,ineq,sol,coeff,coeff_rnd,err_sol,err_pnt,err_rnd], | |
| ratepsilon:1e-4, | |
| xp:rat(Chebyshev_Extrema(degree+1,xmin,xmax)), | |
| ratepsilon:eps_y, | |
| yp:makelist(rat(f(float(xp[i]))),i,1,degree+2), | |
| ydiff:makelist(yp[i]-sum((a[k])*xp[i]^k,k,0,degree),i,1,degree+2), | |
| gep:makelist(umax>=ydiff[i],i,1,degree+2), | |
| gem:makelist(umax>=-ydiff[i],i,1,degree+2), | |
| ineq:append(gep,gem), | |
| sol:second(minimize_lp(umax,ineq,[umax])), | |
| coeff:at(makelist(a[i],i,0,degree),sol), | |
| coeff_rnd:round(coeff*denom)/denom, | |
| err_sol:at(umax,sol), | |
| err_pnt:lmax(abs(at(ydiff,sol))), | |
| err_rnd:lmax(abs(at(ydiff,makelist(a[i]=coeff_rnd[i+1],i,0,degree)))), | |
| /* | |
| print("err_sol = ", err_sol, float(err_sol)), | |
| print("err_pnt = ", err_pnt, float(err_pnt)), | |
| print("err_rnd = ", err_rnd, float(err_rnd)), | |
| print("rnd/pnt-1 = ", float(err_rnd/err_pnt-1)), | |
| */ | |
| [coeff,coeff_rnd, err_pnt,err_rnd] | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approxomate a function rescaled to x in -1/2..1/2 | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approximation_poly_linf_rescale(f_full,degree,x_min,x_max,eps_y,denom,max_deriv):=block( | |
| [x_center, x_scale, y_scale, f_scaled], | |
| x_center:(x_max+x_min)/2, | |
| x_scale:(x_max-x_min), | |
| y_scale:1/(max_deriv*x_scale), | |
| f_scaled(x):=y_scale*f_full(x_scale*x+x_center), | |
| [coeff,coeff_rnd, err_pnt,err_rnd]:approximation_poly_linf_cheb(f_scaled,degree,-1/2,1/2,eps_y*y_scale,denom), | |
| [coeff,coeff_rnd, err_pnt/y_scale,err_rnd/y_scale] | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Utility | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| log2(x):=if is(float(x)=0.0) then -1000 else log(float(x))/log(2.0)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approximate all functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approx_all(xmin,xmax):= makelist( | |
| approximation_poly_linf_rescale(funcs[i],degree,xmin,xmax,2^-(50),denoms,max_deriv), | |
| i,1,length(funcs))$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| pos:makelist((x_max-x_min)*i/n_segments+x_min,i,0,n_segments)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| res:makelist(approx_all(pos[i],pos[i+1]),i,1,length(pos)-1)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Evaluate precision | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: subsubsect start ] | |
| Coefficients not rounded | |
| [wxMaxima: subsubsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],float(third(res[r][c]))),length(pos)-1,length(funcs)); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| In bits | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],log2(third(res[r][c]))),length(pos)-1,length(funcs)); | |
| maxerr:lmax(flatten(args(%))); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsubsect start ] | |
| Coefficients rounded | |
| [wxMaxima: subsubsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],float(fourth(res[r][c]))),length(pos)-1,length(funcs)); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| In bits | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],log2(fourth(res[r][c]))),length(pos)-1,length(funcs)); | |
| maxerr:lmax(flatten(args(%))); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Plot | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: input start ] */ | |
| nbits_plot:31$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| colors:[blue,red]$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| makepoly(coeff,x_min,x_max,max_deriv):=block([center,x_scale], | |
| x_center:(x_max+x_min)/2, | |
| x_scale:1/(x_max-x_min), | |
| y_scale:max_deriv/x_scale, | |
| expand(y_scale*sum(float(coeff[i])*(x_scale*(x-x_center))^(i-1),i,1,length(coeff))) | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approximated functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| wxdraw2d( | |
| makelist( | |
| makelist( | |
| [ | |
| color=colors[ifunc], | |
| explicit(makepoly(second(res[ipos][ifunc]),pos[ipos],pos[ipos+1],max_deriv),x,pos[ipos],pos[ipos+1]) | |
| ], | |
| ipos,1,length(pos)-1), | |
| ifunc,1,length(funcs) | |
| ) | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Deviation in 31 bit LSBs (rounded coefficients, but assuming precise real arithmetic) | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| makelist( | |
| wxdraw2d( | |
| makelist( | |
| [ | |
| color=colors[ifunc], | |
| explicit(2^nbits_plot*(funcs[ifunc](x)-makepoly(first(res[ipos][ifunc]),pos[ipos],pos[ipos+1],max_deriv)),x,pos[ipos],pos[ipos+1]) | |
| ], | |
| ipos,1,length(pos)-1), | |
| color=gray,explicit(2^(maxerr+nbits_plot),x,x_min,x_max),explicit(-2^(maxerr+nbits_plot),x,x_min,x_max) | |
| ), | |
| ifunc,1,length(funcs) | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Dump coefficients | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coq format | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| for ifunc from 1 thru length(funcs) do ( | |
| printf(true,"Definition coeff_~a : list (list Z):= [~%", funcs[ifunc]), | |
| for ipos from 1 thru length(pos)-1 do ( | |
| printf(true, " ["), | |
| vals:second(res[ipos][ifunc])*denoms, | |
| for ival from 1 thru length(vals) do (printf(true, "~d", vals[ival]), if ival<length(vals) then printf(true, "; ")), | |
| printf(true, "]~a~%",if ipos<length(pos)-1 then ";" else "") | |
| ), | |
| printf(true, "].~%") | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Final results for comparison | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coefficients rounded | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| err_sin:float(lmax(makelist(fourth(res[r][1]),r,1,length(pos)-1)))$ | |
| err_cos:float(lmax(makelist(fourth(res[r][2]),r,1,length(pos)-1)))$ | |
| printf(false,"Max Error Sin = ~,3e (~,3f bits)~%", err_sin, log2(err_sin)); | |
| printf(false,"Max Error Cos = ~,3e (~,3f bits)~%", err_cos, log2(err_cos)); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coefficients not rounded | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| err_sin:float(lmax(makelist(third(res[r][1]),r,1,length(pos)-1)))$ | |
| err_cos:float(lmax(makelist(third(res[r][2]),r,1,length(pos)-1)))$ | |
| printf(false,"Max Error Sin = ~,3e (~,3f bits)~%", err_sin, log2(err_sin)); | |
| printf(false,"Max Error Cos = ~,3e (~,3f bits)~%", err_cos, log2(err_cos)); | |
| /* [wxMaxima: input end ] */ | |
| /* Old versions of Maxima abort on loading files that end in a comment. */ | |
| "Created with wxMaxima 21.02.0"$ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ | |
| /* [ Created with wxMaxima version 21.02.0 ] */ | |
| /* [wxMaxima: title start ] | |
| Least L∞ approximation of | |
| Sine and Cosine in 0...1 | |
| 64 pieces with degree 3 | |
| [wxMaxima: title end ] */ | |
| /* [wxMaxima: section start ] | |
| Intro | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| funcs:[sin,cos]; | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Main parameter | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| degree:3$ | |
| n_segments:64$ | |
| x_min:0$ | |
| x_max:1$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Scalings | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| Scaling denominators for x^0 ... x^n | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: comment start ] | |
| Note: giving 6 extra bits for c0 increases precision after rounding to to 37.2 bits | |
| Note: this might be a case where ILP is worthwhile | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| nbits:31$ | |
| max_deriv:1$ | |
| denoms:makelist(2^nbits,i,0,degree); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Libraries | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| load(simplex)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Settings | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| Some corner cases involve 0^0 terms. | |
| 0^0 is usually undefined, but we tell the maxima simplifier to replace it with 1 here. | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| simp : false $ | |
| tellsimp (0^0, 1) $ | |
| tellsimp (0.0^0, 1) $ | |
| simp : true $ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| The simplex algorithm epsilon can be set to 0 for rational input | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| epsilon_lp:0$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: comment start ] | |
| Setting for float to rat | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| ratepsilon:1.0e-18$ | |
| ratprint:false$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Linf approximation of a function sampled at the extrema of the Chebyshev polynomial with given degree | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: comment start ] | |
| We put the samples at the extrema of the Chebeshev poylnomial, because the residual is usually close to a chebychev polynomial | |
| Chebyshev extrama: http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html | |
| [wxMaxima: comment end ] */ | |
| /* [wxMaxima: input start ] */ | |
| Chebyshev_Extrema(n,xmin,xmax):=makelist(float((cos(%pi*i/n)-1)/-2*(xmax-xmin)+xmin),i,0,n); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approximation_poly_linf_cheb(f,degree,xmin,xmax,eps_y,denom):=block( | |
| [ratepsilon,xp,yp,ydiff,gep,gem,ineq,sol,coeff,coeff_rnd,err_sol,err_pnt,err_rnd], | |
| ratepsilon:1e-4, | |
| xp:rat(Chebyshev_Extrema(degree+1,xmin,xmax)), | |
| ratepsilon:eps_y, | |
| yp:makelist(rat(f(float(xp[i]))),i,1,degree+2), | |
| ydiff:makelist(yp[i]-sum((a[k])*xp[i]^k,k,0,degree),i,1,degree+2), | |
| gep:makelist(umax>=ydiff[i],i,1,degree+2), | |
| gem:makelist(umax>=-ydiff[i],i,1,degree+2), | |
| ineq:append(gep,gem), | |
| sol:second(minimize_lp(umax,ineq,[umax])), | |
| coeff:at(makelist(a[i],i,0,degree),sol), | |
| coeff_rnd:round(coeff*denom)/denom, | |
| err_sol:at(umax,sol), | |
| err_pnt:lmax(abs(at(ydiff,sol))), | |
| err_rnd:lmax(abs(at(ydiff,makelist(a[i]=coeff_rnd[i+1],i,0,degree)))), | |
| /* | |
| print("err_sol = ", err_sol, float(err_sol)), | |
| print("err_pnt = ", err_pnt, float(err_pnt)), | |
| print("err_rnd = ", err_rnd, float(err_rnd)), | |
| print("rnd/pnt-1 = ", float(err_rnd/err_pnt-1)), | |
| */ | |
| [coeff,coeff_rnd, err_pnt,err_rnd] | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approxomate a function rescaled to x in -1/2..1/2 | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approximation_poly_linf_rescale(f_full,degree,x_min,x_max,eps_y,denom,max_deriv):=block( | |
| [x_center, x_scale, y_scale, f_scaled], | |
| x_center:(x_max+x_min)/2, | |
| x_scale:(x_max-x_min), | |
| y_scale:1/(max_deriv*x_scale), | |
| f_scaled(x):=y_scale*f_full(x_scale*x+x_center), | |
| [coeff,coeff_rnd, err_pnt,err_rnd]:approximation_poly_linf_cheb(f_scaled,degree,-1/2,1/2,eps_y*y_scale,denom), | |
| [coeff,coeff_rnd, err_pnt/y_scale,err_rnd/y_scale] | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Utility | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| log2(x):=if is(float(x)=0.0) then -1000 else log(float(x))/log(2.0)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approximate all functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| approx_all(xmin,xmax):= makelist( | |
| approximation_poly_linf_rescale(funcs[i],degree,xmin,xmax,2^-(50),denoms,max_deriv), | |
| i,1,length(funcs))$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| pos:makelist((x_max-x_min)*i/n_segments+x_min,i,0,n_segments)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| res:makelist(approx_all(pos[i],pos[i+1]),i,1,length(pos)-1)$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Evaluate precision | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: subsubsect start ] | |
| Coefficients not rounded | |
| [wxMaxima: subsubsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],log2(third(res[r][c]))),length(pos)-1,length(funcs)); | |
| maxerr:lmax(flatten(args(%))); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsubsect start ] | |
| Coefficients rounded | |
| [wxMaxima: subsubsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| genmatrix(lambda([r,c],log2(fourth(res[r][c]))),length(pos)-1,length(funcs)); | |
| maxerr:lmax(flatten(args(%))); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Plot | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: input start ] */ | |
| nbits_plot:31$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| colors:[blue,red]$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: input start ] */ | |
| makepoly(coeff,x_min,x_max,max_deriv):=block([center,x_scale], | |
| x_center:(x_max+x_min)/2, | |
| x_scale:1/(x_max-x_min), | |
| y_scale:max_deriv/x_scale, | |
| expand(y_scale*sum(float(coeff[i])*(x_scale*(x-x_center))^(i-1),i,1,length(coeff))) | |
| )$ | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Approximated functions | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| wxdraw2d( | |
| makelist( | |
| makelist( | |
| [ | |
| color=colors[ifunc], | |
| explicit(makepoly(second(res[ipos][ifunc]),pos[ipos],pos[ipos+1],max_deriv),x,pos[ipos],pos[ipos+1]) | |
| ], | |
| ipos,1,length(pos)-1), | |
| ifunc,1,length(funcs) | |
| ) | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Deviation in 31 bit LSBs (rounded coefficients, but assuming precise real arithmetic) | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| makelist( | |
| wxdraw2d( | |
| makelist( | |
| [ | |
| color=colors[ifunc], | |
| explicit(2^nbits_plot*(funcs[ifunc](x)-makepoly(first(res[ipos][ifunc]),pos[ipos],pos[ipos+1],max_deriv)),x,pos[ipos],pos[ipos+1]) | |
| ], | |
| ipos,1,length(pos)-1), | |
| color=gray,explicit(2^(maxerr+nbits_plot),x,x_min,x_max),explicit(-2^(maxerr+nbits_plot),x,x_min,x_max) | |
| ), | |
| ifunc,1,length(funcs) | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Dump coefficients | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coq format | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| for ifunc from 1 thru length(funcs) do ( | |
| printf(true,"Definition coeff_~a : list (list Z):= [~%", funcs[ifunc]), | |
| for ipos from 1 thru length(pos)-1 do ( | |
| printf(true, " ["), | |
| vals:second(res[ipos][ifunc])*denoms, | |
| for ival from 1 thru length(vals) do (printf(true, "~d", vals[ival]), if ival<length(vals) then printf(true, "; ")), | |
| printf(true, "]~a~%",if ipos<length(pos)-1 then ";" else "") | |
| ), | |
| printf(true, "].~%") | |
| ); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: section start ] | |
| Final results for comparison | |
| [wxMaxima: section end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coefficients rounded | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| err_sin:float(lmax(makelist(fourth(res[r][1]),r,1,length(pos)-1)))$ | |
| err_cos:float(lmax(makelist(fourth(res[r][2]),r,1,length(pos)-1)))$ | |
| printf(false,"Max Error Sin = ~,3e (~,3f bits)~%", err_sin, log2(err_sin)); | |
| printf(false,"Max Error Cos = ~,3e (~,3f bits)~%", err_cos, log2(err_cos)); | |
| /* [wxMaxima: input end ] */ | |
| /* [wxMaxima: subsect start ] | |
| Coefficients not rounded | |
| [wxMaxima: subsect end ] */ | |
| /* [wxMaxima: input start ] */ | |
| err_sin:float(lmax(makelist(third(res[r][1]),r,1,length(pos)-1)))$ | |
| err_cos:float(lmax(makelist(third(res[r][2]),r,1,length(pos)-1)))$ | |
| printf(false,"Max Error Sin = ~,3e (~,3f bits)~%", err_sin, log2(err_sin)); | |
| printf(false,"Max Error Cos = ~,3e (~,3f bits)~%", err_cos, log2(err_cos)); | |
| /* [wxMaxima: input end ] */ | |
| /* Old versions of Maxima abort on loading files that end in a comment. */ | |
| "Created with wxMaxima 21.02.0"$ |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment