Example of main stream linear program solvers failing on small problems

__init__.py

# For languages which need empty files and servers which don't support this ...

chebyshev_extrema.m

function [xlist] = chebyshev_extrema(n,xmin,xmax)
% Compute an array of x positions of the extrama of a Chebyshev polynomial
% of given degree
xlist = arrayfun(@(i) (cos(pi*i/n)-1)/-2*(xmax-xmin)+xmin, 0:n);
end

Description.txt

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.

linf_min.m

function [coeff,sol_err,max_err,min_err] = linf_min(f,xmin,xmax,deg)
    % Find an approximation polynomial with optimized L-inf norm deviation from
    % a given function
    % f: function to optimize
    % xmin, xmax: range in which the polynomial is optimized
    % deg: degree of polynomial

    xlist = chebyshev_extrema(deg+1,xmin,xmax);

    % Column (variable) names
    % VC: coefficient variable
    % VM: max error variable

    % Row (condition) names
    % CXpos: bound positive error at one x position
    % CXneg: bound negative error at one x position

    A_VC_CXpos = zeros(length(xlist), deg+1);
    A_VC_CXneg = zeros(length(xlist), deg+1);

    A_VM_CXpos = zeros(length(xlist), 1);
    A_VM_CXneg = zeros(length(xlist), 1);

    B_CXpos = zeros(length(xlist), 1);
    B_CXneg = zeros(length(xlist), 1);

    for IX=1:length(xlist)
        for IC=1:deg+1
            % Set matrix elements
            %  + c0 + c1*x + c2*x^3 + c3*x^3 - f(x) <= maxerr
            %  - c0 - c1*x - c2*x^3 - c3*x^3 + f(x) <= maxerr
            % <=>
            % + c0 + c1*x + c2*x^3 + c3*x^3 - maxerr <= f(x)
            % - c0 - c1*x - c2*x^3 - c3*x^3 - maxerr <= -f(x)
            A_VC_CXpos(IX,IC) =   xlist(IX)^(IC-1);
            A_VC_CXneg(IX,IC) = - xlist(IX)^(IC-1);
        end
        A_VM_CXpos(IX,1) = -1;
        A_VM_CXneg(IX,1) = -1;
        B_CXpos(IX,1) =   f(xlist(IX));
        B_CXneg(IX,1) = - f(xlist(IX));
    end

    A = [A_VC_CXpos,A_VM_CXpos; A_VC_CXneg,A_VM_CXneg];
    B = [B_CXpos; B_CXneg];

    cost = [zeros(1,deg+1), 1];

    % options = optimoptions('linprog','Display','iter');
    options = optimoptions('linprog');
    sol = linprog(cost,A,B,[],[],[],[],options);
    coeff = sol(1:deg+1);
    sol_err = sol(deg+2);

    % Compute actual errors at points
    list_devs = arrayfun(@(x) sum(arrayfun(@(i) coeff(i+1) * x^i, 0:deg)) - f(x), xlist);
    max_err = max(abs(list_devs));
    min_err = min(abs(list_devs));
    % Warn if the residual is not close to a Chebyshev Polynomial (or very
    % small)
    if min_err/max_err < 0.99 && max_err > 1e-14 
        disp("WARNING: deviation is not well distributed ");
        disp(list_devs);
    end
    
    % fprintf("Solerr = %.3e, Maxerr = %.3e, Minerr = %.3e\n", sol_err, max_err, min_err);
end

linf_min.py

import glpk
from math import cos,pi,log2

# GLPK C ref see  https://www.gnu.org/software/glpk/
#                 (need to download the .tar.gz to get the PDF)
# PYGLPK see      https://github.com/bradfordboyle/pyglpk
#                 http://tfinley.net/software/pyglpk/

# For a GLPK brief reference see:     http://tfinley.net/software/pyglpk/brief_ref.html
# For a more extensive reference see: https://www.cs.cornell.edu/~tomf/pyglpk/glpk.html#LPX-intopt

def chebyshev_extrema(n,xmin,xmax):
    return [ (cos(pi*i/n)-1)/-2*(xmax-xmin)+xmin for i in range(0,n+1) ]

def linf_min(f,xmin,xmax,deg,method):
    xlist = chebyshev_extrema(deg+1,xmin,xmax);

    lp = glpk.LPX()
    lp.name = 'linf'
    lp.obj.maximize = False

    # Add 2 rows for each x-point, one for yerr<=max and one for -yerr<=max
    lp.rows.add(2*len(xlist))

    for i in range(len(xlist)):
        lp.rows[2*i+0].name = "pos"+str(i)
        lp.rows[2*i+1].name = "neg"+str(i)

    # Add 1 col per coefficient
    ncoeff = deg+1

    # Add 1 maxydiff variable
    imaxerr = ncoeff

    lp.cols.add(ncoeff+1)

    for i in range(ncoeff):
        lp.cols[i].name = 'c%d' % i
        lp.cols[i].kind = int
        lp.cols[i].bounds = (None, None) # -inf <= ci < inf
    
    lp.cols[imaxerr].name = "maxerr"
    lp.cols[imaxerr].bounds = (0, None) # 0 <= dmax < inf

    # Set objective coefficients (minimize dmax)
    lp.obj[:] = [0.0]*(ncoeff)+[1.0]

    # Set matrix elements
    #  + c0 + c1*x + c2*x^3 + c3*x^3 - f(x) <= maxerr
    #  - c0 - c1*x - c2*x^3 - c3*x^3 + f(x) <= maxerr
    # <=>
    # + c0 + c1*x + c2*x^3 + c3*x^3 - maxerr <= f(x)
    # - c0 - c1*x - c2*x^3 - c3*x^3 - maxerr <= -f(x)

    matrix = []
    for i in range(len(xlist)):
        matrix.append( (2*i+0, imaxerr, - 1) )
        matrix.append( (2*i+1, imaxerr, - 1) )
        for j in range(deg+1):
            matrix.append( (2*i+0, j,   xlist[i]**j) )
            matrix.append( (2*i+1, j, - xlist[i]**j) )
        lp.rows[2*i+0].bounds = (None, f(xlist[i]))
        lp.rows[2*i+1].bounds = (None, -f(xlist[i]))
    
    lp.matrix = matrix

    if method==0:
        lp.simplex()         # Solve this LP with the simplex method
    elif method==1:
        lp.interior()        # Solve this LP with the interior method
    else:
        lp.exact()           # Solve this LP with the exact method

    # print('Z = %g;' % lp.obj.value)  # Retrieve and print obj func value
    # print('; '.join('%s = %g' % (c.name, c.primal) for c in lp.cols))

    coeff = [lp.cols[i].primal for i in range(0,deg+1)]
    sol_err = lp.cols[deg+1].primal

    list_devs = [sum(lp.cols[j].primal * x**j for j in range(deg+1)) - f(x) for x in xlist]
    max_err = max([abs(x) for x in list_devs])
    min_err = min([abs(x) for x in list_devs])
    if min_err/max_err < 0.99 and max_err > 1e-14:
        print("WARNING: deviation is not well distributed ")
        print(list_devs)
    return (coeff,sol_err,max_err,min_err)

LinfApprox_Sin_3deg_4seg.m

% Compute linf approximation for Sine and Cosine
% 3rd order in 1/4 intervals

deg = 3;
nseg = 4;

ResultSin = [];
ResultCos = [];

for xpos = 0:(nseg-1)
    xmin = xpos/nseg;
    xmax = (xpos+1)/nseg;

    [coeff,sol_err,max_err,min_err] = linf_min(@sin,xmin,xmax,deg);
    ResultSin = [ResultSin; coeff', max_err, min_err, min_err/max_err];
    
    [coeff,sol_err,max_err,min_err] = linf_min(@cos,xmin,xmax,deg);
    ResultCos = [ResultCos; coeff', max_err, min_err, min_err/max_err];
    
end

fprintf("Max Error Sin = %.3e (%.3f bits)\n", max(ResultSin(:,deg+2)), -log2(max(ResultSin(:,deg+2))));
fprintf("Max Error Cos = %.3e (%.3f bits)\n", max(ResultCos(:,deg+2)), -log2(max(ResultCos(:,deg+2))));

LinfApprox_Sin_3deg_4seg.py

# Compute linf approximation for Sine and Cosine
# 3rd order in 1/4 intervals

from math import sin,cos,log2
from linf_min import *

deg = 3
nseg = 4

ResultSinSi = []
ResultCosSi = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,0)
    ResultSinSi.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,0)
    ResultCosSi.append( (coeff, max_err, min_err, min_err/max_err) )

ResultSinIn = []
ResultCosIn = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,2)
    ResultSinIn.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,2)
    ResultCosIn.append( (coeff, max_err, min_err, min_err/max_err) )

ResultSinEx = []
ResultCosEx = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,2)
    ResultSinEx.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,2)
    ResultCosEx.append( (coeff, max_err, min_err, min_err/max_err) )

err_max_sin_si = max([res[1] for res in ResultSinSi])
err_max_cos_si = max([res[1] for res in ResultCosSi])

err_max_sin_in = max([res[1] for res in ResultSinIn])
err_max_cos_in = max([res[1] for res in ResultCosIn])

err_max_sin_ex = max([res[1] for res in ResultSinEx])
err_max_cos_ex = max([res[1] for res in ResultCosEx])

print( "Max Error Sin Simplex  = {0:.3e} ({1:.3f} bits)".format(err_max_sin_si, log2(err_max_sin_si)))
print( "Max Error Cos Simplex  = {0:.3e} ({1:.3f} bits)".format(err_max_cos_si, log2(err_max_cos_si)))
print( "Max Error Sin Interior = {0:.3e} ({1:.3f} bits)".format(err_max_sin_in, log2(err_max_sin_in)))
print( "Max Error Cos Interior = {0:.3e} ({1:.3f} bits)".format(err_max_cos_in, log2(err_max_cos_in)))
print( "Max Error Sin Exact    = {0:.3e} ({1:.3f} bits)".format(err_max_sin_ex, log2(err_max_sin_ex)))
print( "Max Error Cos Exact    = {0:.3e} ({1:.3f} bits)".format(err_max_cos_ex, log2(err_max_cos_ex)))

LInfApprox_Sin_3deg_4seg.wxm

/* [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"$

LinfApprox_Sin_3deg_64seg.m

% Compute linf approximation for Sine and Cosine
% 3rd order in 1/64 intervals

deg = 3;
nseg = 64;

ResultSin = [];
ResultCos = [];

for xpos = 0:(nseg-1)
    xmin = xpos/nseg;
    xmax = (xpos+1)/nseg;

    [coeff,sol_err,max_err,min_err] = linf_min(@sin,xmin,xmax,deg);
    ResultSin = [ResultSin; coeff', max_err, min_err, min_err/max_err];
    
    [coeff,sol_err,max_err,min_err] = linf_min(@cos,xmin,xmax,deg);
    ResultCos = [ResultCos; coeff', max_err, min_err, min_err/max_err];
    
end

fprintf("Max Error Sin = %.3e (%.3f bits)\n", max(ResultSin(:,deg+2)), -log2(max(ResultSin(:,deg+2))));
fprintf("Max Error Cos = %.3e (%.3f bits)\n", max(ResultCos(:,deg+2)), -log2(max(ResultCos(:,deg+2))));

LinfApprox_Sin_3deg_64seg.py

# Compute linf approximation for Sine and Cosine
# 3rd order in 1/64 intervals

from math import sin,cos,log2
from linf_min import *

deg = 3
nseg = 64

ResultSinSi = []
ResultCosSi = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,0)
    ResultSinSi.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,0)
    ResultCosSi.append( (coeff, max_err, min_err, min_err/max_err) )

ResultSinIn = []
ResultCosIn = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,2)
    ResultSinIn.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,2)
    ResultCosIn.append( (coeff, max_err, min_err, min_err/max_err) )

ResultSinEx = []
ResultCosEx = []

for xpos in range(0,nseg):
    xmin = xpos/nseg
    xmax = (xpos+1)/nseg

    [coeff,sol_err,max_err,min_err] = linf_min(sin,xmin,xmax,deg,2)
    ResultSinEx.append( (coeff, max_err, min_err, min_err/max_err) )

    [coeff,sol_err,max_err,min_err] = linf_min(cos,xmin,xmax,deg,2)
    ResultCosEx.append( (coeff, max_err, min_err, min_err/max_err) )

err_max_sin_si = max([res[1] for res in ResultSinSi])
err_max_cos_si = max([res[1] for res in ResultCosSi])

err_max_sin_in = max([res[1] for res in ResultSinIn])
err_max_cos_in = max([res[1] for res in ResultCosIn])

err_max_sin_ex = max([res[1] for res in ResultSinEx])
err_max_cos_ex = max([res[1] for res in ResultCosEx])

print( "Max Error Sin Simplex  = {0:.3e} ({1:.3f} bits)".format(err_max_sin_si, log2(err_max_sin_si)))
print( "Max Error Cos Simplex  = {0:.3e} ({1:.3f} bits)".format(err_max_cos_si, log2(err_max_cos_si)))
print( "Max Error Sin Interior = {0:.3e} ({1:.3f} bits)".format(err_max_sin_in, log2(err_max_sin_in)))
print( "Max Error Cos Interior = {0:.3e} ({1:.3f} bits)".format(err_max_cos_in, log2(err_max_cos_in)))
print( "Max Error Sin Exact    = {0:.3e} ({1:.3f} bits)".format(err_max_sin_ex, log2(err_max_sin_ex)))
print( "Max Error Cos Exact    = {0:.3e} ({1:.3f} bits)".format(err_max_cos_ex, log2(err_max_cos_ex)))

LInfApprox_Sin_3deg_64seg.wxm

/* [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"$