Great work, thanks for sharing!
Could you maybe write some short blog post/Jupyter tutorial about the process on coding function in julia (how did you found the polynomial, decided its length, tested for both precisions, optimised for performance etc.)?
I know it’s a lot to ask for, but it could be very useful in order to learn to e.g. how program custom or special function over some range in Julia. I haven’t seen such tutorial anywhere (if anyone knows about one I’ll be thankful).
It doesn’t answer all your questions, but I wrote up an introduction to minimax polynomials for myself as I was teaching myself the concept.[*] Hopefully it’s a good introduction to the range reduction and minimax polynomial approximation, though the example I gives is a much lower accuracy version since it’s purpose was only to clearly show the concepts. (Throwing my own derived coefficients into BetterExp’s implementation, it’s doing only barely worse — <1% less accurate on average, and a peak error 3% larger but that is still less than 0.6ulps inaccurate in absolute terms.)
[*] (In particular, I thought most descriptions I was finding skipped the simple connections between properties of N-th degree polynomials and solving a system of equations to explain why there are N+2 nodes in the fit — it’s there implicitly if you think in terms of degrees of freedom, but I think it was more helpful to build up to minimax polynomials from the perspective of Taylor series.)
That’s a great blogpost! The range reduction for expm1 isn’t one I’ve seen before. When I have some time, I’ll update BetterExp with an expm1 implementation.
@dpsanders I just realized that a slightly modified version of this would give you a version of pow suitible for interval arithmetic.
I’m interested in knowing more of course!
I’m looking into sollya. It seems really powerful, but P = fpminimax(expm1(x),[|1,2,3|],[|double...|],[-1/512;1/512],absolute); gives Error: fpminimax: singular matrix. Any idea why that would happen? I can get lots of other polynomials, but not the one with a fixed 0 as the first coefficient.
omit absolute, add fixed:
dom = [-1/512;1/512];
coefftype = [|double...|];
terms = [|1,2,3|];
P = fpminimax(expm1(x), terms, coefftype, dom, fixed);
hint: when you start up Sollya, the first line should be > x;, this establishes x as the default var name (e.g. in polynomials);
here is information on the settings available:
/*
Global variables controlling numerical computations
`prec`
The precision in significand bits of the float format used in numerical computations.
Commands try to adapt their working precision to have approximately `prec` correct bits.
If, in assignment to `diam`, the value assigned is followed by an exclamation mark, no message indicating the new state is displayed.
`diam`
Parameter used in safe algorithms to control the maximal length of involved intervals.
More precisely, `diam` is relative to the width of the input interval of the command.
[`infnorm`, `checkinfnorm`, `accurateinfnorm`, `integral`, `findzeros`, `supnorm`].
If, in assignment to `diam`, the value assigned is followed by an exclamation mark, no message indicating the new state is displayed.
`supnorm` falls back to bisecting the working interval until safe supremum norm bounds can be computed with the required accuracy or until the width of the subintervals becomes less than `diam` times the original interval I, in which case supnorm fails.
`points`
Its value represents the number of points used in numerical algorithms.
[`dirtyinfnorm`, `dirtyintegral`, `dirtyfindzeros`, `plot`]
Global variables controlling recursive refinement
`taylorrecursions`
The number of recursive steps used when applying Taylor's rule.
This rule is applied by the interval evaluator (used in `infnorm`, `evaluate`).
If, in assignment to `diam`, the value assigned is followed by an exclamation mark, no message indicating the new state is displayed.
`hopitalrecursions`
The number of recursive steps used when applying L'Hopital's rule.
This rule is applied by the interval evaluator (used in `infnorm`, `findzeros`, `evaluate`).
If, in assignment to `diam`, the value assigned is followed by an exclamation mark, no message indicating the new state is displayed.
Global variables controlling Sollya system behavior
`dieonerrormode`
Controls if Sollya is exited on an error or not.
`rationalmode`
Controls whether or not rational arithmetic is used.
When its value is off, which is the default, Sollya will not use rational arithmetic to simplify expressions. All computations, including the evaluation of constant expressions given on the Sollya prompt, will be performed using floating-point and interval arithmetic. Constant expressions will be approximated by floating-point numbers, which are in most cases faithful roundings of the expressions, when shown at the prompt.
When the value of the global variable rationalmode is on, Sollya will use rational arithmetic when simplifying expressions. Constant expressions, given at the Sollya prompt, will be simplified to rational numbers and displayed as such when they are in the set of the rational numbers. Otherwise, flaoting-point and interval arithmetic will be used to compute a floating-point approximation, which is in most cases a faithful rounding of the constant expression.
`display`
sets or inspects the global variable specifying number notation
one of decimal|binary|dyadic|powers|hexadecimal
If the global notation variable display is decimal, all numbers will be output in scientific decimal notation. If the global notation variable display is dyadic, all numbers will be output as dyadic numbers with Gappa notation. If the global notation variable display is powers, all numbers will be output as dyadic numbers with a notation compatible with Maple and PARI/GP. If the global notation variable display is binary, all numbers will be output in binary notation. If the global notation variable display is hexadecimal, all numbers will be output in C99/ IEEE754-2008 notation. All output notations can be parsed back by Sollya, inducing no error if the input and output precisions are the same (see prec).
`midpointmode`
Controls the way intervals are displayed.
When its value is off, intervals are displayed as usual (in the form [a;b]). When its value is on, and if a and b have the same first significant digits, the interval in displayed in a way that lets one immediately see the common digits of the two bounds.
This mode is supported only with display set to decimal. In other modes of display, midpointmode value is simply ignored.
`fullparentheses`
An assignment fullparentheses = activation value, where activation value is one of on or off, activates respectively deactivates the output of expressions with full parenthesising. In full parenthesising mode, Sollya commands like print, write and the implicit command when an expression is given at the prompt will output expressions with parenthesising at all places where it is necessary for expressions containing infix operators to be parsed back with the same result. Otherwise parentheses around associative operators are omitted.
`autosimplify`
An assignment autosimplify = activation value, where activation value is one of on or off, activates respectively deactivates the automatic safe simplification of expressions of functions generated by the evaluation of commands or in argument of other commands.
Sollya commands like remez, taylor or rationalapprox sometimes produce expressions that can be simplified. Constant subexpressions can be evaluated to dyadic floating-point numbers, monomials with coefficients 0 can be eliminated. Further, expressions indicated by the user perform better in many commands when simplified before being passed in argument to a command. When the automatic simplification of expressions is activated, Sollya automatically performs a safe (not value changing) simplification process on such expressions.
The automatic generation of subexpressions can be annoying, in particular if it takes too much time for not enough benefit. Further the user might want to inspect the structure of the expression tree returned by a command. In this case, the automatic simplification should be deactivated.
`verbosity`
Controls the amount of information displayed by commands.
At level 0, commands do not display anything on standard output. Note that very critical information may however be displayed on standard error.
Default level is 1. It displays important information such as warnings when roundings happen. For higher levels more information is displayed depending on the command.
`roundingwarnings`
Controls whether or not a warning is displayed when roundings occur.
`timing`
Controls whether or not timing is measured.
When its value is on, the time spent in each command is measured and displayed (for verbosity levels higher than 1).
Related facilities
`on` `off`
settings available for:
[`autosimplify`, `canonical`, `timing`, `fullparentheses`, `midpointmode`, `rationalmode`, `roundingwarnings`, `timing`, `dieonerrormode`]
hint: after generating a polynomial, setting > display = dyadic; makes it easy to see the actual (machine) values of the coeffs you determine. restore the display mode with > display = default;.
These may be useful.
procedure poly2coeffs(ply) {
var dg, cfs;
var i;
dg = degree(ply);
cfs = [| |];
for i from 0 to dg do {
cfs = cfs :. coeff(ply, i);
};
return cfs;
};
procedure coeffs2poly(cfs) {
var n, ply, i;
n = length(cfs) - 1;
ply = cfs[0];
for i from 1 to n do {
ply = ply + cfs[i]*(x^i);
};
return ply;
};
procedure coeffs2dyadic(cfs) {
var cfsexp, cfsfrac, res;
cfsexp = [||];
cfsfrac = [||];
for v in cfs do {
cfsexp = cfsexp :. exponent(v);
cfsfrac = cfsfrac :. mantissa(v);
};
res = [| cfsfrac, cfsexp |];
return res;
};
Here are some other non-numeric procedures that can be handy.
/* control the display mode with less typing */
procedure dec(){ display=decimal!; };
procedure hex(){ display=hexadecimal!; };
procedure bin(){ display=binary!; };
procedure dyad(){ display=dyadic!; };
procedure pows(){ display=powers!; };
/* type testing */
/*
isnumber(x)
isinterval(x) -- zero-width intervals are numbers (false here)
islist(x)
ispoly(x) -- degree 0 polynomials are numbers (false here)
*/
isnumber = proc(a) { inf(a) == sup(a); };
isinterval = proc(a) { inf(a) != sup(a); };
islist = proc(a) { var n; n=length(a); return n == n; };
ispoly = proc(a) { var n; n=degree(a); return n == max(n,1); };
/* essential constructs */
/*
ivl(a,b) -> interval(a,b)
lst(...) -> list(...)
*/
ivl = proc(a, b){ return [min(a,b); max(a,b)]; };
lst = proc(args = ...) {
var res, i;
res = [| |];
for i in args do res = res :. i;
return res;
};
/* basic list handling */
/*
reverse(list)
prepend(item, list)
postpend(list, item)
append(list, list)
dropfirst(list)
droplast(list)
dropidx(list, idx) -- idx==0 drops first, idx==length(list)-1 drops last
*/
procedure reverse(lst) {
var res, lastidx, i;
res = [| |];
lastidx = length(lst) - 1;
if lastidx >= 0 then {
for i from 0 to lastidx do { res = res :. lst[lastidx-i]; };
};
return res;
};
prepend = proc(itm, lst) { return [|itm|] @ lst; };
postpend = proc(lst, itm) { return lst :. itm; };
append = proc(lst1, lst2) { return lst1 @ lst2; };
dropfirst = proc(lst) { return tail(lst); };
droplast = proc(lst) { return reverse(tail(reverse(lst))); };
procedure dropidx(lst, idx) {
var res, lastidx, i;
res = [| |];
lastidx = length(lst) - 1;
if (idx==0) then res = dropfirst(lst) else {
if (idx==lastidx) then res = droplast(lst) else {
for i from 0 to lastidx do { if (i!=idx) then res = res :. lst[i]; };
};
};
return res;
};
I’ve just pushed a new commit that fixes Inf handling and greatly simplifies the constant definitions by using generated functions.
Edit: committing and pushing are not the same thing. I’ve now pushed the changes.
This Package is now officially in the general registry. Hopefully it will be redundant soon.