4 # remez.jl - implementation of the Remez algorithm for polynomial approximation
6 # Copyright (c) 2015-2019, Arm Limited.
7 # SPDX-License-Identifier: MIT
11 # ----------------------------------------------------------------------
12 # Helper functions to cope with different Julia versions.
13 if VERSION >= v"0.7.0"
14 array1d(T, d) = Array{T, 1}(undef, d)
15 array2d(T, d1, d2) = Array{T, 2}(undef, d1, d2)
17 array1d(T, d) = Array(T, d)
18 array2d(T, d1, d2) = Array(T, d1, d2)
23 if VERSION >= v"0.6.0"
24 # Use Base.invokelatest to run functions made using eval(), to
25 # avoid "world age" error
26 run(f, x...) = Base.invokelatest(f, x...)
28 # Prior to 0.6.0, invokelatest doesn't exist (but fortunately the
29 # world age problem also doesn't seem to exist)
30 run(f, x...) = f(x...)
33 # ----------------------------------------------------------------------
34 # Global variables configured by command-line options.
35 floatsuffix = "" # adjusted by --floatsuffix
36 xvarname = "x" # adjusted by --variable
37 epsbits = 256 # adjusted by --bits
38 debug_facilities = Set() # adjusted by --debug
39 full_output = false # adjusted by --full
40 array_format = false # adjusted by --array
41 preliminary_commands = array1d(String, 0) # adjusted by --pre
43 # ----------------------------------------------------------------------
44 # Diagnostic and utility functions.
46 # Enable debugging printouts from a particular subpart of this
50 # facility Name of the facility to debug. For a list of facility names,
51 # look through the code for calls to debug().
53 # Return value is a BigFloat.
54 function enable_debug(facility)
55 push!(debug_facilities, facility)
61 # facility Name of the facility for which this is a debug message.
62 # printargs Arguments to println() if debugging of that facility is
64 macro debug(facility, printargs...)
66 print("[", $facility, "] ")
75 if $facility in debug_facilities
82 # Evaluate a polynomial.
85 # coeffs Array of BigFloats giving the coefficients of the polynomial.
86 # Starts with the constant term, i.e. coeffs[i] is the
87 # coefficient of x^(i-1) (because Julia arrays are 1-based).
88 # x Point at which to evaluate the polynomial.
90 # Return value is a BigFloat.
91 function poly_eval(coeffs::Array{BigFloat}, x::BigFloat)
98 return coeffs[1] + x * poly_eval(coeffs[2:n], x)
102 # Evaluate a rational function.
105 # ncoeffs Array of BigFloats giving the coefficients of the numerator.
106 # Starts with the constant term, and 1-based, as above.
107 # dcoeffs Array of BigFloats giving the coefficients of the denominator.
108 # Starts with the constant term, and 1-based, as above.
109 # x Point at which to evaluate the function.
111 # Return value is a BigFloat.
112 function ratfn_eval(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat},
114 return poly_eval(ncoeffs, x) / poly_eval(dcoeffs, x)
117 # Format a BigFloat into an appropriate output format.
119 # x BigFloat to format.
121 # Return value is a string.
122 function float_to_str(x)
123 return string(x) * floatsuffix
126 # Format a polynomial into an arithmetic expression, for pasting into
127 # other tools such as gnuplot.
130 # coeffs Array of BigFloats giving the coefficients of the polynomial.
131 # Starts with the constant term, and 1-based, as above.
133 # Return value is a string.
134 function poly_to_string(coeffs::Array{BigFloat})
139 return float_to_str(coeffs[1])
141 return string(float_to_str(coeffs[1]), "+", xvarname, "*(",
142 poly_to_string(coeffs[2:n]), ")")
146 # Format a rational function into a string.
149 # ncoeffs Array of BigFloats giving the coefficients of the numerator.
150 # Starts with the constant term, and 1-based, as above.
151 # dcoeffs Array of BigFloats giving the coefficients of the denominator.
152 # Starts with the constant term, and 1-based, as above.
154 # Return value is a string.
155 function ratfn_to_string(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat})
156 if length(dcoeffs) == 1 && dcoeffs[1] == 1
157 # Special case: if the denominator is just 1, leave it out.
158 return poly_to_string(ncoeffs)
160 return string("(", poly_to_string(ncoeffs), ")/(",
161 poly_to_string(dcoeffs), ")")
165 # Format a list of x,y pairs into a string.
168 # xys Array of (x,y) pairs of BigFloats.
170 # Return value is a string.
171 function format_xylist(xys::Array{Tuple{BigFloat,BigFloat}})
173 join([" "*string(x)*" -> "*string(y) for (x,y) in xys], "\n") *
177 # ----------------------------------------------------------------------
178 # Matrix-equation solver for matrices of BigFloat.
180 # I had hoped that Julia's type-genericity would allow me to solve the
181 # matrix equation Mx=V by just writing 'M \ V'. Unfortunately, that
182 # works by translating the inputs into double precision and handing
183 # off to an optimised library, which misses the point when I have a
184 # matrix and vector of BigFloat and want my result in _better_ than
185 # double precision. So I have to implement my own specialisation of
186 # the \ operator for that case.
188 # Fortunately, the point of using BigFloats is that we have precision
189 # to burn, so I can do completely naïve Gaussian elimination without
190 # worrying about instability.
193 # matrix_in 2-dimensional array of BigFloats, representing a matrix M
194 # in row-first order, i.e. matrix_in[r,c] represents the
195 # entry in row r col c.
196 # vector_in 1-dimensional array of BigFloats, representing a vector V.
198 # Return value: a 1-dimensional array X of BigFloats, satisfying M X = V.
200 # Expects the input to be an invertible square matrix and a vector of
201 # the corresponding size, on pain of failing an assertion.
202 function \(matrix_in :: Array{BigFloat,2},
203 vector_in :: Array{BigFloat,1})
204 # Copy the inputs, because we'll be mutating them as we go.
208 # Input consistency criteria: matrix is square, and vector has
212 @assert(size(M) == (n,n))
214 @debug("gausselim", "starting, n=", n)
217 # Straightforward Gaussian elimination: find the largest
218 # non-zero entry in column i (and in a row we haven't sorted
219 # out already), swap it into row i, scale that row to
220 # normalise it to 1, then zero out the rest of the column by
221 # subtracting a multiple of that row from each other row.
223 @debug("gausselim", "matrix=", repr(M))
224 @debug("gausselim", "vector=", repr(V))
226 # Find the best pivot.
230 if abs(M[j,i]) > bestval
235 @assert(bestrow > 0) # make sure we did actually find one
237 @debug("gausselim", "bestrow=", bestrow)
239 # Swap it into row i.
242 M[bestrow,k],M[i,k] = M[i,k],M[bestrow,k]
244 V[bestrow],V[i] = V[i],V[bestrow]
247 # Scale that row so that M[i,i] becomes 1.
250 M[i,k] = M[i,k] / divisor
252 V[i] = V[i] / divisor
255 # Zero out all other entries in column i, by subtracting
256 # multiples of this row.
261 M[j,k] = M[j,k] - M[i,k] * factor
263 V[j] = V[j] - V[i] * factor
269 @debug("gausselim", "matrix=", repr(M))
270 @debug("gausselim", "vector=", repr(V))
271 @debug("gausselim", "done!")
273 # Now we're done: M is the identity matrix, so the equation Mx=V
274 # becomes just x=V, i.e. V is already exactly the vector we want
279 # ----------------------------------------------------------------------
280 # Least-squares fitting of a rational function to a set of (x,y)
283 # We use this to get an initial starting point for the Remez
284 # iteration. Therefore, it doesn't really need to be particularly
285 # accurate; it only needs to be good enough to wiggle back and forth
286 # across the target function the right number of times (so as to give
287 # enough error extrema to start optimising from) and not have any
288 # poles in the target interval.
290 # Least-squares fitting of a _polynomial_ is actually a sensible thing
291 # to do, and minimises the rms error. Doing the following trick with a
292 # rational function P/Q is less sensible, because it cannot be made to
293 # minimise the error function (P/Q-f)^2 that you actually wanted;
294 # instead it minimises (P-fQ)^2. But that should be good enough to
295 # have the properties described above.
297 # Some theory: suppose you're trying to choose a set of parameters a_i
298 # so as to minimise the sum of squares of some error function E_i.
299 # Basic calculus says, if you do this in one variable, just
300 # differentiate and solve for zero. In this case, that works fine even
301 # with multiple variables, because you _partially_ differentiate with
302 # respect to each a_i, giving a system of equations, and that system
303 # turns out to be linear so we just solve it as a matrix.
305 # In this case, our parameters are the coefficients of P and Q; to
306 # avoid underdetermining the system we'll fix Q's constant term at 1,
307 # so that our error function (as described above) is
309 # E = \sum (p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d)^2
311 # where the sum is over all (x,y) coordinate pairs. Setting dE/dp_j=0
312 # (for each j) gives an equation of the form
314 # 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) x^j
316 # and setting dE/dq_j=0 gives one of the form
318 # 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) y x^j
320 # And both of those row types, treated as multivariate linear
321 # equations in the p,q values, have each coefficient being a value of
322 # the form \sum x^i, \sum y x^i or \sum y^2 x^i, for various i. (Times
323 # a factor of 2, but we can throw that away.) So we can go through the
324 # list of input coordinates summing all of those things, and then we
325 # have enough information to construct our matrix and solve it
326 # straight off for the rational function coefficients.
329 # f The function to be approximated. Maps BigFloat -> BigFloat.
330 # xvals Array of BigFloats, giving the list of x-coordinates at which
332 # n Degree of the numerator polynomial of the desired rational
334 # d Degree of the denominator polynomial of the desired rational
336 # w Error-weighting function. Takes two BigFloat arguments x,y
337 # and returns a scaling factor for the error at that location.
338 # A larger value indicates that the error should be given
339 # greater weight in the square sum we try to minimise.
340 # If unspecified, defaults to giving everything the same weight.
342 # Return values: a pair of arrays of BigFloats (N,D) giving the
343 # coefficients of the returned rational function. N has size n+1; D
344 # has size d+1. Both start with the constant term, i.e. N[i] is the
345 # coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
347 function ratfn_leastsquares(f::Function, xvals::Array{BigFloat}, n, d,
348 w = (x,y)->BigFloat(1))
349 # Accumulate sums of x^i y^j, for j={0,1,2} and a range of x.
350 # Again because Julia arrays are 1-based, we'll have sums[i,j]
351 # being the sum of x^(i-1) y^(j-1).
352 maxpow = max(n,d) * 2 + 1
353 sums = zeros(BigFloat, maxpow, 3)
359 sums[i,j] += x^(i-1) * y^(j-1) * weight
364 @debug("leastsquares", "sums=", repr(sums))
366 # Build the matrix. We're solving n+d+1 equations in n+d+1
367 # unknowns. (We actually have to return n+d+2 coefficients, but
368 # one of them is hardwired to 1.)
369 matrix = array2d(BigFloat, n+d+1, n+d+1)
370 vector = array1d(BigFloat, n+d+1)
372 # Equation obtained by differentiating with respect to p_i,
373 # i.e. the numerator coefficient of x^i.
376 matrix[row, 1+j] = sums[1+i+j, 1]
379 matrix[row, 1+n+j] = -sums[1+i+j, 2]
381 vector[row] = sums[1+i, 2]
384 # Equation obtained by differentiating with respect to q_i,
385 # i.e. the denominator coefficient of x^i.
388 matrix[row, 1+j] = sums[1+i+j, 2]
391 matrix[row, 1+n+j] = -sums[1+i+j, 3]
393 vector[row] = sums[1+i, 3]
396 @debug("leastsquares", "matrix=", repr(matrix))
397 @debug("leastsquares", "vector=", repr(vector))
399 # Solve the matrix equation.
400 all_coeffs = matrix \ vector
402 @debug("leastsquares", "all_coeffs=", repr(all_coeffs))
404 # And marshal the results into two separate polynomial vectors to
406 ncoeffs = all_coeffs[1:n+1]
407 dcoeffs = vcat([1], all_coeffs[n+2:n+d+1])
408 return (ncoeffs, dcoeffs)
411 # ----------------------------------------------------------------------
412 # Golden-section search to find a maximum of a function.
415 # f Function to be maximised/minimised. Maps BigFloat -> BigFloat.
416 # a,b,c BigFloats bracketing a maximum of the function.
419 # a,b,c are in order (either a<=b<=c or c<=b<=a)
420 # a != c (but b can equal one or the other if it wants to)
421 # f(a) <= f(b) >= f(c)
423 # Return value is an (x,y) pair of BigFloats giving the extremal input
424 # and output. (That is, y=f(x).)
425 function goldensection(f::Function, a::BigFloat, b::BigFloat, c::BigFloat)
426 # Decide on a 'good enough' threshold.
427 threshold = abs(c-a) * 2^(-epsbits/2)
429 # We'll need the golden ratio phi, of course. Or rather, in this
430 # case, we need 1/phi = 0.618...
431 one_over_phi = 2 / (1 + sqrt(BigFloat(5)))
433 # Flip round the interval endpoints so that the interval [a,b] is
434 # at least as large as [b,c]. (Then we can always pick our new
435 # point in [a,b] without having to handle lots of special cases.)
436 if abs(b-a) < abs(c-a)
440 # Evaluate the function at the initial points.
445 @debug("goldensection", "starting")
447 while abs(c-a) > threshold
448 @debug("goldensection", "a: ", a, " -> ", fa)
449 @debug("goldensection", "b: ", b, " -> ", fb)
450 @debug("goldensection", "c: ", c, " -> ", fc)
453 @assert(a <= b <= c || c <= b <= a)
454 @assert(fa <= fb >= fc)
456 # Subdivide the larger of the intervals [a,b] and [b,c]. We've
457 # arranged that this is always [a,b], for simplicity.
458 d = a + (b-a) * one_over_phi
460 # Now we have an interval looking like this (possibly
465 # and we know f(b) is bigger than either f(a) or f(c). We have
466 # two cases: either f(d) > f(b), or vice versa. In either
467 # case, we can narrow to an interval of 1/phi the size, and
468 # still satisfy all our invariants (three ordered points,
469 # [a,b] at least the width of [b,c], f(a)<=f(b)>=f(c)).
471 @debug("goldensection", "d: ", d, " -> ", fd)
474 fa, fb, fc = fa, fd, fb
475 @debug("goldensection", "adb case")
478 fa, fb, fc = fc, fb, fd
479 @debug("goldensection", "cbd case")
483 @debug("goldensection", "done: ", b, " -> ", fb)
487 # ----------------------------------------------------------------------
488 # Find the extrema of a function within a given interval.
491 # f The function to be approximated. Maps BigFloat -> BigFloat.
492 # grid A set of points at which to evaluate f. Must be high enough
493 # resolution to make extrema obvious.
495 # Returns an array of (x,y) pairs of BigFloats, with each x,y giving
496 # the extremum location and its value (i.e. y=f(x)).
497 function find_extrema(f::Function, grid::Array{BigFloat})
499 extrema = array1d(Tuple{BigFloat, BigFloat}, 0)
501 # We have to provide goldensection() with three points
502 # bracketing the extremum. If the extremum is at one end of
503 # the interval, then the only way we can do that is to set two
504 # of the points equal (which goldensection() will cope with).
508 # Find our three pairs of (x,y) coordinates.
509 xp, xi, xn = grid[prev], grid[i], grid[next]
510 yp, yi, yn = f(xp), f(xi), f(xn)
512 # See if they look like an extremum, and if so, ask
513 # goldensection() to give a more exact location for it.
515 push!(extrema, goldensection(f, xp, xi, xn))
516 elseif yp >= yi <= yn
517 x, y = goldensection(x->-f(x), xp, xi, xn)
518 push!(extrema, (x, -y))
524 # ----------------------------------------------------------------------
525 # Winnow a list of a function's extrema to give a subsequence of a
526 # specified length, with the extrema in the subsequence alternating
527 # signs, and with the smallest absolute value of an extremum in the
528 # subsequence as large as possible.
530 # We do this using a dynamic-programming approach. We work along the
531 # provided array of extrema, and at all times, we track the best set
532 # of extrema we have so far seen for each possible (length, sign of
533 # last extremum) pair. Each new extremum is evaluated to see whether
534 # it can be added to any previously seen best subsequence to make a
535 # new subsequence that beats the previous record holder in its slot.
538 # extrema An array of (x,y) pairs of BigFloats giving the input extrema.
539 # n Number of extrema required as output.
541 # Returns a new array of (x,y) pairs which is a subsequence of the
542 # original sequence. (So, in particular, if the input was sorted by x
543 # then so will the output be.)
544 function winnow_extrema(extrema::Array{Tuple{BigFloat,BigFloat}}, n)
545 # best[i,j] gives the best sequence so far of length i and with
546 # sign j (where signs are coded as 1=positive, 2=negative), in the
547 # form of a tuple (cost, actual array of x,y pairs).
548 best = fill((BigFloat(0), array1d(Tuple{BigFloat,BigFloat}, 0)), n, 2)
556 # A zero-valued extremum cannot possibly contribute to any
557 # optimal sequence, so we simply ignore it!
562 # See if we can create a new entry for best[i,sign] by
563 # appending our current (x,y) to some previous thing.
565 # Special case: we don't store a best zero-length
567 candidate = (abs(y), [(x,y)])
569 othersign = 3-sign # map 1->2 and 2->1
570 oldscore, oldlist = best[i-1, othersign]
571 newscore = min(abs(y), oldscore)
572 newlist = vcat(oldlist, [(x,y)])
573 candidate = (newscore, newlist)
575 # If our new candidate improves on the previous value of
576 # best[i,sign], then replace it.
577 if candidate[1] > best[i,sign][1]
578 best[i,sign] = candidate
583 # Our ultimate return value has to be either best[n,1] or
584 # best[n,2], but it could be either. See which one has the higher
586 if best[n,1][1] > best[n,2][1]
591 # Make sure we did actually _find_ a good answer.
592 @assert(length(ret) == n)
596 # ----------------------------------------------------------------------
597 # Construct a rational-function approximation with equal and
598 # alternating weighted deviation at a specific set of x-coordinates.
601 # f The function to be approximated. Maps BigFloat -> BigFloat.
602 # coords An array of BigFloats giving the x-coordinates. There should
604 # n, d The degrees of the numerator and denominator of the desired
606 # prev_err A plausible value for the alternating weighted deviation.
607 # (Required to kickstart a binary search in the nonlinear case;
608 # see comments below.)
609 # w Error-weighting function. Takes two BigFloat arguments x,y
610 # and returns a scaling factor for the error at that location.
611 # The returned approximation R should have the minimum possible
612 # maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
613 # parameter, defaulting to the always-return-1 function.
615 # Return values: a pair of arrays of BigFloats (N,D) giving the
616 # coefficients of the returned rational function. N has size n+1; D
617 # has size d+1. Both start with the constant term, i.e. N[i] is the
618 # coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
620 function ratfn_equal_deviation(f::Function, coords::Array{BigFloat},
621 n, d, prev_err::BigFloat,
622 w = (x,y)->BigFloat(1))
623 @debug("equaldev", "n=", n, " d=", d, " coords=", repr(coords))
624 @assert(length(coords) == n+d+2)
627 # Special case: we're after a polynomial. In this case, we
628 # have the particularly easy job of just constructing and
629 # solving a system of n+2 linear equations, to find the n+1
630 # coefficients of the polynomial and also the amount of
631 # deviation at the specified coordinates. Each equation is of
634 # p_0 x^0 + p_1 x^1 + ... + p_n x^n ± e/w(x) = f(x)
636 # in which the p_i and e are the variables, and the powers of
637 # x and calls to w and f are the coefficients.
639 matrix = array2d(BigFloat, n+2, n+2)
640 vector = array1d(BigFloat, n+2)
649 matrix[i, n+2] = currsign / w(x,y)
653 @debug("equaldev", "matrix=", repr(matrix))
654 @debug("equaldev", "vector=", repr(vector))
656 outvector = matrix \ vector
658 @debug("equaldev", "outvector=", repr(outvector))
660 ncoeffs = outvector[1:n+1]
661 dcoeffs = [BigFloat(1)]
662 return ncoeffs, dcoeffs
664 # For a nontrivial rational function, the system of equations
665 # we need to solve becomes nonlinear, because each equation
668 # p_0 x^0 + p_1 x^1 + ... + p_n x^n
669 # --------------------------------- ± e/w(x) = f(x)
670 # x^0 + q_1 x^1 + ... + q_d x^d
672 # and multiplying up by the denominator gives you a lot of
673 # terms containing e × q_i. So we can't do this the really
674 # easy way using a matrix equation as above.
676 # Fortunately, this is a fairly easy kind of nonlinear system.
677 # The equations all become linear if you switch to treating e
678 # as a constant, so a reasonably sensible approach is to pick
679 # a candidate value of e, solve all but one of the equations
680 # for the remaining unknowns, and then see what the error
681 # turns out to be in the final equation. The Chebyshev
682 # alternation theorem guarantees that that error in the last
683 # equation will be anti-monotonic in the input e, so we can
684 # just binary-search until we get the two as close to equal as
688 # Try a given value of e, derive the coefficients of the
689 # resulting rational function by setting up equations
690 # based on the first n+d+1 of the n+d+2 coordinates, and
691 # see what the error turns out to be at the final
693 matrix = array2d(BigFloat, n+d+1, n+d+1)
694 vector = array1d(BigFloat, n+d+1)
699 y_adj = y - currsign * e / w(x,y)
704 matrix[i,1+n+j] = -x^j * y_adj
710 @debug("equaldev", "trying e=", e)
711 @debug("equaldev", "matrix=", repr(matrix))
712 @debug("equaldev", "vector=", repr(vector))
714 outvector = matrix \ vector
716 @debug("equaldev", "outvector=", repr(outvector))
718 ncoeffs = outvector[1:n+1]
719 dcoeffs = vcat([BigFloat(1)], outvector[n+2:n+d+1])
723 last_e = (ratfn_eval(ncoeffs, dcoeffs, x) - y) * w(x,y) * -currsign
725 @debug("equaldev", "last e=", last_e)
727 return ncoeffs, dcoeffs, last_e
730 threshold = 2^(-epsbits/2) # convergence threshold
732 # Start by trying our previous iteration's error value. This
733 # value (e0) will be one end of our binary-search interval,
734 # and whatever it caused the last point's error to be, that
735 # (e1) will be the other end.
737 @debug("equaldev", "e0 = ", e0)
738 nc, dc, e1 = try_e(e0)
739 @debug("equaldev", "e1 = ", e1)
740 if abs(e1-e0) <= threshold
741 # If we're _really_ lucky, we hit the error right on the
742 # nose just by doing that!
746 @debug("equaldev", "s = ", s)
748 # Verify by assertion that trying our other interval endpoint
749 # e1 gives a value that's wrong in the other direction.
750 # (Otherwise our binary search won't get a sensible answer at
752 nc, dc, e2 = try_e(e1)
753 @debug("equaldev", "e2 = ", e2)
754 @assert(sign(e2-e1) == -s)
756 # Now binary-search until our two endpoints narrow enough.
758 while abs(e1-e0) > threshold
760 nc, dc, enew = try_e(emid)
761 if sign(enew-emid) == s
768 @debug("equaldev", "final e=", emid)
773 # ----------------------------------------------------------------------
774 # Top-level function to find a minimax rational-function approximation.
777 # f The function to be approximated. Maps BigFloat -> BigFloat.
778 # interval A pair of BigFloats giving the endpoints of the interval
779 # (in either order) on which to approximate f.
780 # n, d The degrees of the numerator and denominator of the desired
782 # w Error-weighting function. Takes two BigFloat arguments x,y
783 # and returns a scaling factor for the error at that location.
784 # The returned approximation R should have the minimum possible
785 # maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
786 # parameter, defaulting to the always-return-1 function.
788 # Return values: a tuple (N,D,E,X), where
790 # N,D A pair of arrays of BigFloats giving the coefficients
791 # of the returned rational function. N has size n+1; D
792 # has size d+1. Both start with the constant term, i.e.
793 # N[i] is the coefficient of x^(i-1) (because Julia
794 # arrays are 1-based). D[1] will be 1.
795 # E The maximum weighted error (BigFloat).
796 # X An array of pairs of BigFloats giving the locations of n+2
797 # points and the weighted error at each of those points. The
798 # weighted error values will have alternating signs, which
799 # means that the Chebyshev alternation theorem guarantees
800 # that any other function of the same degree must exceed
801 # the error of this one at at least one of those points.
802 function ratfn_minimax(f::Function, interval::Tuple{BigFloat,BigFloat}, n, d,
803 w = (x,y)->BigFloat(1))
804 # We start off by finding a least-squares approximation. This
805 # doesn't need to be perfect, but if we can get it reasonably good
806 # then it'll save iterations in the refining stage.
808 # Least-squares approximations tend to look nicer in a minimax
809 # sense if you evaluate the function at a big pile of Chebyshev
810 # nodes rather than uniformly spaced points. These values will
811 # also make a good grid to use for the initial search for error
812 # extrema, so we'll keep them around for that reason too.
814 # Construct the grid.
815 lo, hi = minimum(interval), maximum(interval)
820 nnodes = 16 * (n+d+1)
821 pi = 2*asin(BigFloat(1))
822 grid = [ mid - halfwid * cos(pi*i/nnodes) for i=0:1:nnodes ]
825 # Find the initial least-squares approximation.
826 (nc, dc) = ratfn_leastsquares(f, grid, n, d, w)
827 @debug("minimax", "initial leastsquares approx = ",
828 ratfn_to_string(nc, dc))
830 # Threshold of convergence. We stop when the relative difference
831 # between the min and max (winnowed) error extrema is less than
834 # This is set to the cube root of machine epsilon on a more or
835 # less empirical basis, because the rational-function case will
836 # not converge reliably if you set it to only the square root.
837 # (Repeatable by using the --test mode.) On the assumption that
838 # input and output error in each iteration can be expected to be
839 # related by a simple power law (because it'll just be down to how
840 # many leading terms of a Taylor series are zero), the cube root
841 # was the next thing to try.
842 threshold = 2^(-epsbits/3)
846 # Find all the error extrema we can.
847 function compute_error(x)
849 approx_y = ratfn_eval(nc, dc, x)
850 return (approx_y - real_y) * w(x, real_y)
852 extrema = find_extrema(compute_error, grid)
853 @debug("minimax", "all extrema = ", format_xylist(extrema))
855 # Winnow the extrema down to the right number, and ensure they
856 # have alternating sign.
857 extrema = winnow_extrema(extrema, n+d+2)
858 @debug("minimax", "winnowed extrema = ", format_xylist(extrema))
860 # See if we've finished.
861 min_err = minimum([abs(y) for (x,y) = extrema])
862 max_err = maximum([abs(y) for (x,y) = extrema])
863 variation = (max_err - min_err) / max_err
864 @debug("minimax", "extremum variation = ", variation)
865 if variation < threshold
866 @debug("minimax", "done!")
867 return nc, dc, max_err, extrema
870 # If not, refine our function by equalising the error at the
871 # extrema points, and go round again.
872 (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema),
874 @debug("minimax", "refined approx = ", ratfn_to_string(nc, dc))
878 # ----------------------------------------------------------------------
879 # Check if a polynomial is well-conditioned for accurate evaluation in
880 # a given interval by Horner's rule.
882 # This is true if at every step where Horner's rule computes
883 # (coefficient + x*value_so_far), the constant coefficient you're
884 # adding on is of larger magnitude than the x*value_so_far operand.
885 # And this has to be true for every x in the interval.
888 # coeffs The coefficients of the polynomial under test. Starts with
889 # the constant term, i.e. coeffs[i] is the coefficient of
890 # x^(i-1) (because Julia arrays are 1-based).
891 # lo, hi The bounds of the interval.
893 # Return value: the largest ratio (x*value_so_far / coefficient), at
894 # any step of evaluation, for any x in the interval. If this is less
895 # than 1, the polynomial is at least somewhat well-conditioned;
896 # ideally you want it to be more like 1/8 or 1/16 or so, so that the
897 # relative rounding error accumulated at each step are reduced by
898 # several factors of 2 when the next coefficient is added on.
900 function wellcond(coeffs, lo, hi)
901 x = max(abs(lo), abs(hi))
904 for i = length(coeffs):-1:1
905 coeff = abs(coeffs[i])
908 thisval = so_far / coeff
909 worst = max(worst, thisval)
916 # ----------------------------------------------------------------------
917 # Small set of unit tests.
923 function approx_eq(x, y, limit=1e-6)
924 return abs(x - y) < limit
927 function test(condition)
936 # Test Gaussian elimination.
937 println("Gaussian test 1:")
938 m = BigFloat[1 1 2; 3 5 8; 13 34 21]
939 v = BigFloat[1, -1, 2]
941 println(" ",repr(ret))
942 test(approx_eq(ret[1], 109/26))
943 test(approx_eq(ret[2], -105/130))
944 test(approx_eq(ret[3], -31/26))
946 # Test leastsquares rational functions.
947 println("Leastsquares test 1:")
949 a = array1d(BigFloat, n+1)
951 a[1+i] = i/BigFloat(n)
953 (nc, dc) = ratfn_leastsquares(x->exp(x), a, 2, 2)
954 println(" ",ratfn_to_string(nc, dc))
956 test(approx_eq(exp(x), ratfn_eval(nc, dc, x), 1e-4))
959 # Test golden section search.
960 println("Golden section test 1:")
961 x, y = goldensection(x->sin(x),
962 BigFloat(0), BigFloat(1)/10, BigFloat(4))
963 println(" ", x, " -> ", y)
964 test(approx_eq(x, asin(BigFloat(1))))
965 test(approx_eq(y, 1))
967 # Test extrema-winnowing algorithm.
968 println("Winnow test 1:")
969 extrema = [(x, sin(20*x)*sin(197*x))
970 for x in BigFloat(0):BigFloat(1)/1000:BigFloat(1)]
971 winnowed = winnow_extrema(extrema, 12)
972 println(" ret = ", format_xylist(winnowed))
977 test(prevy * y <= 0) # tolerates initial prevx having no sign
982 # Test actual minimax approximation.
983 println("Minimax test 1 (polynomial):")
984 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0)
986 println(" ",ratfn_to_string(nc, dc))
988 for x = 0:BigFloat(1)/1000:1
989 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)
992 println("Minimax test 2 (rational):")
993 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2)
995 println(" ",ratfn_to_string(nc, dc))
997 for x = 0:BigFloat(1)/1000:1
998 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)
1001 println("Minimax test 3 (polynomial, weighted):")
1002 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0,
1005 println(" ",ratfn_to_string(nc, dc))
1007 for x = 0:BigFloat(1)/1000:1
1008 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1011 println("Minimax test 4 (rational, weighted):")
1012 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2,
1015 println(" ",ratfn_to_string(nc, dc))
1017 for x = 0:BigFloat(1)/1000:1
1018 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1021 println("Minimax test 5 (rational, weighted, odd degree):")
1022 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 1,
1025 println(" ",ratfn_to_string(nc, dc))
1027 for x = 0:BigFloat(1)/1000:1
1028 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1031 total = passes + fails
1032 println(passes, " passes ", fails, " fails ", total, " total")
1035 # ----------------------------------------------------------------------
1041 remez.jl [options] <lo> <hi> <n> <d> <expr> [<weight>]
1047 Bounds of the interval on which to approximate the target
1048 function. These are parsed and evaluated as Julia expressions,
1049 so you can write things like '1/BigFloat(6)' to get an
1050 accurate representation of 1/6, or '4*atan(BigFloat(1))' to
1051 get pi. (Unfortunately, the obvious 'BigFloat(pi)' doesn't
1056 The desired degree of polynomial(s) you want for your
1057 approximation. These should be non-negative integers. If you
1058 want a rational function as output, set <n> to the degree of
1059 the numerator, and <d> the denominator. If you just want an
1060 ordinary polynomial, set <d> to 0, and <n> to the degree of
1061 the polynomial you want.
1065 A Julia expression giving the function to be approximated on
1066 the interval. The input value is predefined as 'x' when this
1067 expression is evaluated, so you should write something along
1068 the lines of 'sin(x)' or 'sqrt(1+tan(x)^2)' etc.
1072 If provided, a Julia expression giving the weighting factor
1073 for the approximation error. The output polynomial will
1074 minimise the largest absolute value of (P-f) * w at any point
1075 in the interval, where P is the value of the polynomial, f is
1076 the value of the target function given by <expr>, and w is the
1077 weight given by this function.
1079 When this expression is evaluated, the input value to P and f
1080 is predefined as 'x', and also the true output value f(x) is
1081 predefined as 'y'. So you can minimise the relative error by
1082 simply writing '1/y'.
1084 If the <weight> argument is not provided, the default
1085 weighting function always returns 1, so that the polynomial
1086 will minimise the maximum absolute error |P-f|.
1088 Computation options:
1092 Evaluate the Julia expression <predef_expr> before starting
1093 the computation. This permits you to pre-define variables or
1094 functions which the Julia expressions in your main arguments
1095 can refer to. All of <lo>, <hi>, <expr> and <weight> can make
1096 use of things defined by <predef_expr>.
1098 One internal remez.jl function that you might sometimes find
1099 useful in this expression is 'goldensection', which finds the
1100 location and value of a maximum of a function. For example,
1101 one implementation strategy for the gamma function involves
1102 translating it to put its unique local minimum at the origin,
1103 in which case you can write something like this
1105 --pre='(m,my) = goldensection(x -> -gamma(x),
1106 BigFloat(1), BigFloat(1.5), BigFloat(2))'
1108 to predefine 'm' as the location of gamma's minimum, and 'my'
1109 as the (negated) value that gamma actually takes at that
1110 point, i.e. -gamma(m).
1112 (Since 'goldensection' always finds a maximum, we had to
1113 negate gamma in the input function to make it find a minimum
1114 instead. Consult the comments in the source for more details
1115 on the use of this function.)
1117 If you use this option more than once, all the expressions you
1118 provide will be run in sequence.
1122 Specify the accuracy to which you want the output polynomial,
1123 in bits. Default 256, which should be more than enough.
1125 --bigfloatbits=<bits>
1127 Turn up the precision used by Julia for its BigFloat
1128 evaluation. Default is Julia's default (also 256). You might
1129 want to try setting this higher than the --bits value if the
1130 algorithm is failing to converge for some reason.
1136 Instead of just printing the approximation function itself,
1137 also print auxiliary information:
1138 - the locations of the error extrema, and the actual
1139 (weighted) error at each of those locations
1140 - the overall maximum error of the function
1141 - a 'well-conditioning quotient', giving the worst-case ratio
1142 between any polynomial coefficient and the largest possible
1143 value of the higher-order terms it will be added to.
1145 The well-conditioning quotient should be less than 1, ideally
1146 by several factors of two, for accurate evaluation in the
1147 target precision. If you request a rational function, a
1148 separate well-conditioning quotient will be printed for the
1149 numerator and denominator.
1151 Use this option when deciding how wide an interval to
1152 approximate your function on, and what degree of polynomial
1155 --variable=<identifier>
1157 When writing the output polynomial or rational function in its
1158 usual form as an arithmetic expression, use <identifier> as
1159 the name of the input variable. Default is 'x'.
1163 When writing the output polynomial or rational function in its
1164 usual form as an arithmetic expression, write <suffix> after
1165 every floating-point literal. For example, '--suffix=F' will
1166 generate a C expression in which the coefficients are literals
1167 of type 'float' rather than 'double'.
1171 Instead of writing the output polynomial as an arithmetic
1172 expression in Horner's rule form, write out just its
1173 coefficients, one per line, each with a trailing comma.
1174 Suitable for pasting into a C array declaration.
1176 This option is not currently supported if the output is a
1177 rational function, because you'd need two separate arrays for
1178 the numerator and denominator coefficients and there's no
1179 obviously right way to provide both of those together.
1181 Debug and test options:
1185 Enable debugging output from various parts of the Remez
1186 calculation. <facility> should be the name of one of the
1187 classes of diagnostic output implemented in the program.
1188 Useful values include 'gausselim', 'leastsquares',
1189 'goldensection', 'equaldev', 'minimax'. This is probably
1190 mostly useful to people debugging problems with the script, so
1191 consult the source code for more information about what the
1192 diagnostic output for each of those facilities will be.
1194 If you want diagnostics from more than one facility, specify
1195 this option multiple times with different arguments.
1199 Run remez.jl's internal test suite. No arguments needed.
1201 Miscellaneous options:
1205 Display this text and exit. No arguments needed.
1210 # ----------------------------------------------------------------------
1214 nargs = length(argwords)
1215 if nargs != 5 && nargs != 6
1216 error("usage: remez.jl <lo> <hi> <n> <d> <expr> [<weight>]\n" *
1217 " run 'remez.jl --help' for more help")
1220 for preliminary_command in preliminary_commands
1221 eval(Meta.parse(preliminary_command))
1224 lo = BigFloat(eval(Meta.parse(argwords[1])))
1225 hi = BigFloat(eval(Meta.parse(argwords[2])))
1226 n = parse(Int,argwords[3])
1227 d = parse(Int,argwords[4])
1228 f = eval(Meta.parse("x -> " * argwords[5]))
1230 # Wrap the user-provided function with a function of our own. This
1231 # arranges to detect silly FP values (inf,nan) early and diagnose
1232 # them sensibly, and also lets us log all evaluations of the
1233 # function in case you suspect it's doing the wrong thing at some
1234 # special-case point.
1237 @debug("f", x, " -> ", y)
1239 error("f(" * string(x) * ") returned non-finite value " * string(y))
1245 # Wrap the user-provided weight function similarly.
1246 w = eval(Meta.parse("(x,y) -> " * argwords[6]))
1247 function wrapped_weight(x,y)
1250 error("w(" * string(x) * "," * string(y) *
1251 ") returned non-finite value " * string(ww))
1255 weight = wrapped_weight
1257 weight = (x,y)->BigFloat(1)
1260 (nc, dc, e, extrema) = ratfn_minimax(func, (lo, hi), n, d, weight)
1263 functext = join([string(x)*",\n" for x=nc],"")
1265 # It's unclear how you should best format an array of
1266 # coefficients for a rational function, so I'll leave
1267 # implementing this option until I have a use case.
1268 error("--array unsupported for rational functions")
1271 functext = ratfn_to_string(nc, dc) * "\n"
1274 # Print everything you might want to know about the function
1275 println("extrema = ", format_xylist(extrema))
1276 println("maxerror = ", string(e))
1278 println("wellconditioning_numerator = ",
1279 string(wellcond(nc, lo, hi)))
1280 println("wellconditioning_denominator = ",
1281 string(wellcond(dc, lo, hi)))
1283 println("wellconditioning = ", string(wellcond(nc, lo, hi)))
1285 print("function = ", functext)
1287 # Just print the text people will want to paste into their code
1292 # ----------------------------------------------------------------------
1293 # Top-level code: parse the argument list and decide what to do.
1298 argwords = array1d(String, 0)
1300 global doing_opts, what_to_do, argwords
1301 global full_output, array_format, xvarname, floatsuffix, epsbits
1302 if doing_opts && startswith(arg, "-")
1305 elseif arg == "--help"
1307 elseif arg == "--test"
1309 elseif arg == "--full"
1311 elseif arg == "--array"
1313 elseif startswith(arg, "--debug=")
1314 enable_debug(arg[length("--debug=")+1:end])
1315 elseif startswith(arg, "--variable=")
1316 xvarname = arg[length("--variable=")+1:end]
1317 elseif startswith(arg, "--suffix=")
1318 floatsuffix = arg[length("--suffix=")+1:end]
1319 elseif startswith(arg, "--bits=")
1320 epsbits = parse(Int,arg[length("--bits=")+1:end])
1321 elseif startswith(arg, "--bigfloatbits=")
1322 set_bigfloat_precision(
1323 parse(Int,arg[length("--bigfloatbits=")+1:end]))
1324 elseif startswith(arg, "--pre=")
1325 push!(preliminary_commands, arg[length("--pre=")+1:end])
1327 error("unrecognised option: ", arg)
1330 push!(argwords, arg)