simonbyrne / Remez.jl

Remez algorithm for computing minimax polynomial approximations
Other
38 stars 9 forks source link

Feature Request: ability to specify even (or odd) polynomials. #7

Open oscardssmith opened 4 years ago

oscardssmith commented 4 years ago

For even/odd functions using the respective type of polynomial typically allows for much higher degree approximations. Would it be relatively easy to add the ability to generate a minimax polynomial with certain terms fixed? Thanks!

simonbyrne commented 4 years ago

I'm not sure if there is a way to fix coefficients, but for you should be able to transform the function, intervals and weights to get what you want, e.g. consider the following:

julia> ratfn_minimax(x -> cos(x), (-1/512,1/512), 6, 0,(x,y) -> x)[1]
7-element Array{BigFloat,1}:
  1.000000000000000000000000021025191453910635870262113046244880090250052048297687
  2.132329974586648507759707957826779622191991364617559271619507547428078552377105e-74
 -0.5000000000000000000437265904769655641659729882711213035055769055850619707633793
 -1.122545207285226850711688473959848539555991331752233164686532858130512495042655e-67
  0.04166666666668928739987928126653566886960053981443506915647204043580989982175183
  2.796152308103702412184446266691273904168636496802224849332706841668390414291359e-62
 -0.001388892152376249734890636889354352572279127984505145883593751672434589118164789

to constrain to an even polynomial, you can transform it to:

julia> ratfn_minimax(x -> cos(sqrt(x)), (0,(1/512)^2), 3, 0,(x,y) -> sqrt(x))[1]
4-element Array{BigFloat,1}:
  0.9999999999999999999999999998153614910642167623930369077428126620089776605273766
 -0.4999999999999999999993546416298739241880404263787103007106799039667498668265747
  0.04166666666666605763009027306279404917620420148441538830588253416884656588659239
 -0.001388888676015171583619608294278269866265742478420274017161303267093260299410331

Odd functions are a bit harder, but it should be possible.

simonbyrne commented 4 years ago

Ah, for an odd function you can convert:

julia> ratfn_minimax(x -> sin(x), (-1/512,1/512), 7, 0,(x,y) -> x)[1]
8-element Array{BigFloat,1}:
 -3.695361030594939689728566991198485091100479970705568016827899105540576263301655e-72
  0.9999999999999999999999999999592745870846129634839722846553382529810933278600539
  5.219441535098615558452287887243420528828018205259353548678162935558488268091111e-66
 -0.1666666666666666666665629408475538969981944769165884652287143574206390546711979
 -2.226090510734476363345472695821139622950093077348705584070966871199051379258093e-60
  0.008333333333333251215355947763285648320671917326106380601207847239202092029266208
  2.914486859514740250310194236135653306346099728098587389472686231335601990767873e-55
 -0.0001984126727912604329760708576819678483891849401418629735426438096911363155412933

to

julia> ratfn_minimax(x -> sin(sqrt(x))/sqrt(x), (1e-20,(1/512)^2), 3, 0,(x,y) -> sqrt(x))[1]
4-element Array{BigFloat,1}:
  0.9999999999999999999999999999794846097010352881605084357836376676395804324813592
 -0.1666666666666666666665949601797662501932435609968704815390181041970388588285878
  0.008333333333333265662601526813530838660751790339238225095391057055136981595896283
 -0.0001984126747600628730137219626995982625405010136389965654192310584210894971590516
nsajko commented 1 year ago

It may be interesting to note that the Sollya tool, which finds polynomial (not rational, sadly) approximations using (among other things) the Remez algorithm, does support this feature (in fact it's even more general than that).

https://sollya.org

https://sollya.org/sollya-weekly/help.php?name=fpminimax

A list of integers indicates the list of desired monomials. For instance, the list [|0,2,4,6|] indicates that the polynomial must be even and of degree at most 6. Giving an integer n as second argument is equivalent as giving [|0,...,n|]. Finally, a list of function g_k indicates that the desired approximation must be a linear combination of the g_k.

https://sollya.org/sollya-weekly/help.php?name=remez

If L is given and is a list of integers, p is searched as a linear combination of monomials X^k where k belongs to L. In the case when L is a list of integers, it may contain ellipses but cannot be end-elliptic. If L is given and is a list of functions g_k, p is searched as a linear combination of the g_k.