1336 lines · plain
1#!/usr/bin/env julia2# -*- julia -*-3 4# remez.jl - implementation of the Remez algorithm for polynomial approximation5#6# Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.7# See https://llvm.org/LICENSE.txt for license information.8# SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception9 10import Base.\11 12# ----------------------------------------------------------------------13# Helper functions to cope with different Julia versions.14if VERSION >= v"0.7.0"15 array1d(T, d) = Array{T, 1}(undef, d)16 array2d(T, d1, d2) = Array{T, 2}(undef, d1, d2)17else18 array1d(T, d) = Array(T, d)19 array2d(T, d1, d2) = Array(T, d1, d2)20end21if VERSION < v"0.5.0"22 String = ASCIIString23end24if VERSION >= v"0.6.0"25 # Use Base.invokelatest to run functions made using eval(), to26 # avoid "world age" error27 run(f, x...) = Base.invokelatest(f, x...)28else29 # Prior to 0.6.0, invokelatest doesn't exist (but fortunately the30 # world age problem also doesn't seem to exist)31 run(f, x...) = f(x...)32end33 34# ----------------------------------------------------------------------35# Global variables configured by command-line options.36floatsuffix = "" # adjusted by --floatsuffix37xvarname = "x" # adjusted by --variable38epsbits = 256 # adjusted by --bits39debug_facilities = Set() # adjusted by --debug40full_output = false # adjusted by --full41array_format = false # adjusted by --array42preliminary_commands = array1d(String, 0) # adjusted by --pre43 44# ----------------------------------------------------------------------45# Diagnostic and utility functions.46 47# Enable debugging printouts from a particular subpart of this48# program.49#50# Arguments:51# facility Name of the facility to debug. For a list of facility names,52# look through the code for calls to debug().53#54# Return value is a BigFloat.55function enable_debug(facility)56 push!(debug_facilities, facility)57end58 59# Print a diagnostic.60#61# Arguments:62# facility Name of the facility for which this is a debug message.63# printargs Arguments to println() if debugging of that facility is64# enabled.65macro debug(facility, printargs...)66 printit = quote67 print("[", $facility, "] ")68 end69 for arg in printargs70 printit = quote71 $printit72 print($(esc(arg)))73 end74 end75 return quote76 if $facility in debug_facilities77 $printit78 println()79 end80 end81end82 83# Evaluate a polynomial.84 85# Arguments:86# coeffs Array of BigFloats giving the coefficients of the polynomial.87# Starts with the constant term, i.e. coeffs[i] is the88# coefficient of x^(i-1) (because Julia arrays are 1-based).89# x Point at which to evaluate the polynomial.90#91# Return value is a BigFloat.92function poly_eval(coeffs::Array{BigFloat}, x::BigFloat)93 n = length(coeffs)94 if n == 095 return BigFloat(0)96 elseif n == 197 return coeffs[1]98 else99 return coeffs[1] + x * poly_eval(coeffs[2:n], x)100 end101end102 103# Evaluate a rational function.104 105# Arguments:106# ncoeffs Array of BigFloats giving the coefficients of the numerator.107# Starts with the constant term, and 1-based, as above.108# dcoeffs Array of BigFloats giving the coefficients of the denominator.109# Starts with the constant term, and 1-based, as above.110# x Point at which to evaluate the function.111#112# Return value is a BigFloat.113function ratfn_eval(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat},114 x::BigFloat)115 return poly_eval(ncoeffs, x) / poly_eval(dcoeffs, x)116end117 118# Format a BigFloat into an appropriate output format.119# Arguments:120# x BigFloat to format.121#122# Return value is a string.123function float_to_str(x)124 return string(x) * floatsuffix125end126 127# Format a polynomial into an arithmetic expression, for pasting into128# other tools such as gnuplot.129 130# Arguments:131# coeffs Array of BigFloats giving the coefficients of the polynomial.132# Starts with the constant term, and 1-based, as above.133#134# Return value is a string.135function poly_to_string(coeffs::Array{BigFloat})136 n = length(coeffs)137 if n == 0138 return "0"139 elseif n == 1140 return float_to_str(coeffs[1])141 else142 return string(float_to_str(coeffs[1]), "+", xvarname, "*(",143 poly_to_string(coeffs[2:n]), ")")144 end145end146 147# Format a rational function into a string.148 149# Arguments:150# ncoeffs Array of BigFloats giving the coefficients of the numerator.151# Starts with the constant term, and 1-based, as above.152# dcoeffs Array of BigFloats giving the coefficients of the denominator.153# Starts with the constant term, and 1-based, as above.154#155# Return value is a string.156function ratfn_to_string(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat})157 if length(dcoeffs) == 1 && dcoeffs[1] == 1158 # Special case: if the denominator is just 1, leave it out.159 return poly_to_string(ncoeffs)160 else161 return string("(", poly_to_string(ncoeffs), ")/(",162 poly_to_string(dcoeffs), ")")163 end164end165 166# Format a list of x,y pairs into a string.167 168# Arguments:169# xys Array of (x,y) pairs of BigFloats.170#171# Return value is a string.172function format_xylist(xys::Array{Tuple{BigFloat,BigFloat}})173 return ("[\n" *174 join([" "*string(x)*" -> "*string(y) for (x,y) in xys], "\n") *175 "\n]")176end177 178# ----------------------------------------------------------------------179# Matrix-equation solver for matrices of BigFloat.180#181# I had hoped that Julia's type-genericity would allow me to solve the182# matrix equation Mx=V by just writing 'M \ V'. Unfortunately, that183# works by translating the inputs into double precision and handing184# off to an optimised library, which misses the point when I have a185# matrix and vector of BigFloat and want my result in _better_ than186# double precision. So I have to implement my own specialisation of187# the \ operator for that case.188#189# Fortunately, the point of using BigFloats is that we have precision190# to burn, so I can do completely naïve Gaussian elimination without191# worrying about instability.192 193# Arguments:194# matrix_in 2-dimensional array of BigFloats, representing a matrix M195# in row-first order, i.e. matrix_in[r,c] represents the196# entry in row r col c.197# vector_in 1-dimensional array of BigFloats, representing a vector V.198#199# Return value: a 1-dimensional array X of BigFloats, satisfying M X = V.200#201# Expects the input to be an invertible square matrix and a vector of202# the corresponding size, on pain of failing an assertion.203function \(matrix_in :: Array{BigFloat,2},204 vector_in :: Array{BigFloat,1})205 # Copy the inputs, because we'll be mutating them as we go.206 M = copy(matrix_in)207 V = copy(vector_in)208 209 # Input consistency criteria: matrix is square, and vector has210 # length to match.211 n = length(V)212 @assert(n > 0)213 @assert(size(M) == (n,n))214 215 @debug("gausselim", "starting, n=", n)216 217 for i = 1:1:n218 # Straightforward Gaussian elimination: find the largest219 # non-zero entry in column i (and in a row we haven't sorted220 # out already), swap it into row i, scale that row to221 # normalise it to 1, then zero out the rest of the column by222 # subtracting a multiple of that row from each other row.223 224 @debug("gausselim", "matrix=", repr(M))225 @debug("gausselim", "vector=", repr(V))226 227 # Find the best pivot.228 bestrow = 0229 bestval = 0230 for j = i:1:n231 if abs(M[j,i]) > bestval232 bestrow = j233 bestval = M[j,i]234 end235 end236 @assert(bestrow > 0) # make sure we did actually find one237 238 @debug("gausselim", "bestrow=", bestrow)239 240 # Swap it into row i.241 if bestrow != i242 for k = 1:1:n243 M[bestrow,k],M[i,k] = M[i,k],M[bestrow,k]244 end245 V[bestrow],V[i] = V[i],V[bestrow]246 end247 248 # Scale that row so that M[i,i] becomes 1.249 divisor = M[i,i]250 for k = 1:1:n251 M[i,k] = M[i,k] / divisor252 end253 V[i] = V[i] / divisor254 @assert(M[i,i] == 1)255 256 # Zero out all other entries in column i, by subtracting257 # multiples of this row.258 for j = 1:1:n259 if j != i260 factor = M[j,i]261 for k = 1:1:n262 M[j,k] = M[j,k] - M[i,k] * factor263 end264 V[j] = V[j] - V[i] * factor265 @assert(M[j,i] == 0)266 end267 end268 end269 270 @debug("gausselim", "matrix=", repr(M))271 @debug("gausselim", "vector=", repr(V))272 @debug("gausselim", "done!")273 274 # Now we're done: M is the identity matrix, so the equation Mx=V275 # becomes just x=V, i.e. V is already exactly the vector we want276 # to return.277 return V278end279 280# ----------------------------------------------------------------------281# Least-squares fitting of a rational function to a set of (x,y)282# points.283#284# We use this to get an initial starting point for the Remez285# iteration. Therefore, it doesn't really need to be particularly286# accurate; it only needs to be good enough to wiggle back and forth287# across the target function the right number of times (so as to give288# enough error extrema to start optimising from) and not have any289# poles in the target interval.290#291# Least-squares fitting of a _polynomial_ is actually a sensible thing292# to do, and minimises the rms error. Doing the following trick with a293# rational function P/Q is less sensible, because it cannot be made to294# minimise the error function (P/Q-f)^2 that you actually wanted;295# instead it minimises (P-fQ)^2. But that should be good enough to296# have the properties described above.297#298# Some theory: suppose you're trying to choose a set of parameters a_i299# so as to minimise the sum of squares of some error function E_i.300# Basic calculus says, if you do this in one variable, just301# differentiate and solve for zero. In this case, that works fine even302# with multiple variables, because you _partially_ differentiate with303# respect to each a_i, giving a system of equations, and that system304# turns out to be linear so we just solve it as a matrix.305#306# In this case, our parameters are the coefficients of P and Q; to307# avoid underdetermining the system we'll fix Q's constant term at 1,308# so that our error function (as described above) is309#310# E = \sum (p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d)^2311#312# where the sum is over all (x,y) coordinate pairs. Setting dE/dp_j=0313# (for each j) gives an equation of the form314#315# 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) x^j316#317# and setting dE/dq_j=0 gives one of the form318#319# 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) y x^j320#321# And both of those row types, treated as multivariate linear322# equations in the p,q values, have each coefficient being a value of323# the form \sum x^i, \sum y x^i or \sum y^2 x^i, for various i. (Times324# a factor of 2, but we can throw that away.) So we can go through the325# list of input coordinates summing all of those things, and then we326# have enough information to construct our matrix and solve it327# straight off for the rational function coefficients.328 329# Arguments:330# f The function to be approximated. Maps BigFloat -> BigFloat.331# xvals Array of BigFloats, giving the list of x-coordinates at which332# to evaluate f.333# n Degree of the numerator polynomial of the desired rational334# function.335# d Degree of the denominator polynomial of the desired rational336# function.337# w Error-weighting function. Takes two BigFloat arguments x,y338# and returns a scaling factor for the error at that location.339# A larger value indicates that the error should be given340# greater weight in the square sum we try to minimise.341# If unspecified, defaults to giving everything the same weight.342#343# Return values: a pair of arrays of BigFloats (N,D) giving the344# coefficients of the returned rational function. N has size n+1; D345# has size d+1. Both start with the constant term, i.e. N[i] is the346# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will347# be 1.348function ratfn_leastsquares(f::Function, xvals::Array{BigFloat}, n, d,349 w = (x,y)->BigFloat(1))350 # Accumulate sums of x^i y^j, for j={0,1,2} and a range of x.351 # Again because Julia arrays are 1-based, we'll have sums[i,j]352 # being the sum of x^(i-1) y^(j-1).353 maxpow = max(n,d) * 2 + 1354 sums = zeros(BigFloat, maxpow, 3)355 for x = xvals356 y = f(x)357 weight = w(x,y)358 for i = 1:1:maxpow359 for j = 1:1:3360 sums[i,j] += x^(i-1) * y^(j-1) * weight361 end362 end363 end364 365 @debug("leastsquares", "sums=", repr(sums))366 367 # Build the matrix. We're solving n+d+1 equations in n+d+1368 # unknowns. (We actually have to return n+d+2 coefficients, but369 # one of them is hardwired to 1.)370 matrix = array2d(BigFloat, n+d+1, n+d+1)371 vector = array1d(BigFloat, n+d+1)372 for i = 0:1:n373 # Equation obtained by differentiating with respect to p_i,374 # i.e. the numerator coefficient of x^i.375 row = 1+i376 for j = 0:1:n377 matrix[row, 1+j] = sums[1+i+j, 1]378 end379 for j = 1:1:d380 matrix[row, 1+n+j] = -sums[1+i+j, 2]381 end382 vector[row] = sums[1+i, 2]383 end384 for i = 1:1:d385 # Equation obtained by differentiating with respect to q_i,386 # i.e. the denominator coefficient of x^i.387 row = 1+n+i388 for j = 0:1:n389 matrix[row, 1+j] = sums[1+i+j, 2]390 end391 for j = 1:1:d392 matrix[row, 1+n+j] = -sums[1+i+j, 3]393 end394 vector[row] = sums[1+i, 3]395 end396 397 @debug("leastsquares", "matrix=", repr(matrix))398 @debug("leastsquares", "vector=", repr(vector))399 400 # Solve the matrix equation.401 all_coeffs = matrix \ vector402 403 @debug("leastsquares", "all_coeffs=", repr(all_coeffs))404 405 # And marshal the results into two separate polynomial vectors to406 # return.407 ncoeffs = all_coeffs[1:n+1]408 dcoeffs = vcat([1], all_coeffs[n+2:n+d+1])409 return (ncoeffs, dcoeffs)410end411 412# ----------------------------------------------------------------------413# Golden-section search to find a maximum of a function.414 415# Arguments:416# f Function to be maximised/minimised. Maps BigFloat -> BigFloat.417# a,b,c BigFloats bracketing a maximum of the function.418#419# Expects:420# a,b,c are in order (either a<=b<=c or c<=b<=a)421# a != c (but b can equal one or the other if it wants to)422# f(a) <= f(b) >= f(c)423#424# Return value is an (x,y) pair of BigFloats giving the extremal input425# and output. (That is, y=f(x).)426function goldensection(f::Function, a::BigFloat, b::BigFloat, c::BigFloat)427 # Decide on a 'good enough' threshold.428 threshold = abs(c-a) * 2^(-epsbits/2)429 430 # We'll need the golden ratio phi, of course. Or rather, in this431 # case, we need 1/phi = 0.618...432 one_over_phi = 2 / (1 + sqrt(BigFloat(5)))433 434 # Flip round the interval endpoints so that the interval [a,b] is435 # at least as large as [b,c]. (Then we can always pick our new436 # point in [a,b] without having to handle lots of special cases.)437 if abs(b-a) < abs(c-a)438 a, c = c, a439 end440 441 # Evaluate the function at the initial points.442 fa = f(a)443 fb = f(b)444 fc = f(c)445 446 @debug("goldensection", "starting")447 448 while abs(c-a) > threshold449 @debug("goldensection", "a: ", a, " -> ", fa)450 @debug("goldensection", "b: ", b, " -> ", fb)451 @debug("goldensection", "c: ", c, " -> ", fc)452 453 # Check invariants.454 @assert(a <= b <= c || c <= b <= a)455 @assert(fa <= fb >= fc)456 457 # Subdivide the larger of the intervals [a,b] and [b,c]. We've458 # arranged that this is always [a,b], for simplicity.459 d = a + (b-a) * one_over_phi460 461 # Now we have an interval looking like this (possibly462 # reversed):463 #464 # a d b c465 #466 # and we know f(b) is bigger than either f(a) or f(c). We have467 # two cases: either f(d) > f(b), or vice versa. In either468 # case, we can narrow to an interval of 1/phi the size, and469 # still satisfy all our invariants (three ordered points,470 # [a,b] at least the width of [b,c], f(a)<=f(b)>=f(c)).471 fd = f(d)472 @debug("goldensection", "d: ", d, " -> ", fd)473 if fd > fb474 a, b, c = a, d, b475 fa, fb, fc = fa, fd, fb476 @debug("goldensection", "adb case")477 else478 a, b, c = c, b, d479 fa, fb, fc = fc, fb, fd480 @debug("goldensection", "cbd case")481 end482 end483 484 @debug("goldensection", "done: ", b, " -> ", fb)485 return (b, fb)486end487 488# ----------------------------------------------------------------------489# Find the extrema of a function within a given interval.490 491# Arguments:492# f The function to be approximated. Maps BigFloat -> BigFloat.493# grid A set of points at which to evaluate f. Must be high enough494# resolution to make extrema obvious.495#496# Returns an array of (x,y) pairs of BigFloats, with each x,y giving497# the extremum location and its value (i.e. y=f(x)).498function find_extrema(f::Function, grid::Array{BigFloat})499 len = length(grid)500 extrema = array1d(Tuple{BigFloat, BigFloat}, 0)501 for i = 1:1:len502 # We have to provide goldensection() with three points503 # bracketing the extremum. If the extremum is at one end of504 # the interval, then the only way we can do that is to set two505 # of the points equal (which goldensection() will cope with).506 prev = max(1, i-1)507 next = min(i+1, len)508 509 # Find our three pairs of (x,y) coordinates.510 xp, xi, xn = grid[prev], grid[i], grid[next]511 yp, yi, yn = f(xp), f(xi), f(xn)512 513 # See if they look like an extremum, and if so, ask514 # goldensection() to give a more exact location for it.515 if yp <= yi >= yn516 push!(extrema, goldensection(f, xp, xi, xn))517 elseif yp >= yi <= yn518 x, y = goldensection(x->-f(x), xp, xi, xn)519 push!(extrema, (x, -y))520 end521 end522 return extrema523end524 525# ----------------------------------------------------------------------526# Winnow a list of a function's extrema to give a subsequence of a527# specified length, with the extrema in the subsequence alternating528# signs, and with the smallest absolute value of an extremum in the529# subsequence as large as possible.530#531# We do this using a dynamic-programming approach. We work along the532# provided array of extrema, and at all times, we track the best set533# of extrema we have so far seen for each possible (length, sign of534# last extremum) pair. Each new extremum is evaluated to see whether535# it can be added to any previously seen best subsequence to make a536# new subsequence that beats the previous record holder in its slot.537 538# Arguments:539# extrema An array of (x,y) pairs of BigFloats giving the input extrema.540# n Number of extrema required as output.541#542# Returns a new array of (x,y) pairs which is a subsequence of the543# original sequence. (So, in particular, if the input was sorted by x544# then so will the output be.)545function winnow_extrema(extrema::Array{Tuple{BigFloat,BigFloat}}, n)546 # best[i,j] gives the best sequence so far of length i and with547 # sign j (where signs are coded as 1=positive, 2=negative), in the548 # form of a tuple (cost, actual array of x,y pairs).549 best = fill((BigFloat(0), array1d(Tuple{BigFloat,BigFloat}, 0)), n, 2)550 551 for (x,y) = extrema552 if y > 0553 sign = 1554 elseif y < 0555 sign = 2556 else557 # A zero-valued extremum cannot possibly contribute to any558 # optimal sequence, so we simply ignore it!559 continue560 end561 562 for i = 1:1:n563 # See if we can create a new entry for best[i,sign] by564 # appending our current (x,y) to some previous thing.565 if i == 1566 # Special case: we don't store a best zero-length567 # sequence :-)568 candidate = (abs(y), [(x,y)])569 else570 othersign = 3-sign # map 1->2 and 2->1571 oldscore, oldlist = best[i-1, othersign]572 newscore = min(abs(y), oldscore)573 newlist = vcat(oldlist, [(x,y)])574 candidate = (newscore, newlist)575 end576 # If our new candidate improves on the previous value of577 # best[i,sign], then replace it.578 if candidate[1] > best[i,sign][1]579 best[i,sign] = candidate580 end581 end582 end583 584 # Our ultimate return value has to be either best[n,1] or585 # best[n,2], but it could be either. See which one has the higher586 # score.587 if best[n,1][1] > best[n,2][1]588 ret = best[n,1][2]589 else590 ret = best[n,2][2]591 end592 # Make sure we did actually _find_ a good answer.593 @assert(length(ret) == n)594 return ret595end596 597# ----------------------------------------------------------------------598# Construct a rational-function approximation with equal and599# alternating weighted deviation at a specific set of x-coordinates.600 601# Arguments:602# f The function to be approximated. Maps BigFloat -> BigFloat.603# coords An array of BigFloats giving the x-coordinates. There should604# be n+d+2 of them.605# n, d The degrees of the numerator and denominator of the desired606# approximation.607# prev_err A plausible value for the alternating weighted deviation.608# (Required to kickstart a binary search in the nonlinear case;609# see comments below.)610# w Error-weighting function. Takes two BigFloat arguments x,y611# and returns a scaling factor for the error at that location.612# The returned approximation R should have the minimum possible613# maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional614# parameter, defaulting to the always-return-1 function.615#616# Return values: a pair of arrays of BigFloats (N,D) giving the617# coefficients of the returned rational function. N has size n+1; D618# has size d+1. Both start with the constant term, i.e. N[i] is the619# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will620# be 1.621function ratfn_equal_deviation(f::Function, coords::Array{BigFloat},622 n, d, prev_err::BigFloat,623 w = (x,y)->BigFloat(1))624 @debug("equaldev", "n=", n, " d=", d, " coords=", repr(coords))625 @assert(length(coords) == n+d+2)626 627 if d == 0628 # Special case: we're after a polynomial. In this case, we629 # have the particularly easy job of just constructing and630 # solving a system of n+2 linear equations, to find the n+1631 # coefficients of the polynomial and also the amount of632 # deviation at the specified coordinates. Each equation is of633 # the form634 #635 # p_0 x^0 + p_1 x^1 + ... + p_n x^n ± e/w(x) = f(x)636 #637 # in which the p_i and e are the variables, and the powers of638 # x and calls to w and f are the coefficients.639 640 matrix = array2d(BigFloat, n+2, n+2)641 vector = array1d(BigFloat, n+2)642 currsign = +1643 for i = 1:1:n+2644 x = coords[i]645 for j = 0:1:n646 matrix[i,1+j] = x^j647 end648 y = f(x)649 vector[i] = y650 matrix[i, n+2] = currsign / w(x,y)651 currsign = -currsign652 end653 654 @debug("equaldev", "matrix=", repr(matrix))655 @debug("equaldev", "vector=", repr(vector))656 657 outvector = matrix \ vector658 659 @debug("equaldev", "outvector=", repr(outvector))660 661 ncoeffs = outvector[1:n+1]662 dcoeffs = [BigFloat(1)]663 return ncoeffs, dcoeffs664 else665 # For a nontrivial rational function, the system of equations666 # we need to solve becomes nonlinear, because each equation667 # now takes the form668 #669 # p_0 x^0 + p_1 x^1 + ... + p_n x^n670 # --------------------------------- ± e/w(x) = f(x)671 # x^0 + q_1 x^1 + ... + q_d x^d672 #673 # and multiplying up by the denominator gives you a lot of674 # terms containing e × q_i. So we can't do this the really675 # easy way using a matrix equation as above.676 #677 # Fortunately, this is a fairly easy kind of nonlinear system.678 # The equations all become linear if you switch to treating e679 # as a constant, so a reasonably sensible approach is to pick680 # a candidate value of e, solve all but one of the equations681 # for the remaining unknowns, and then see what the error682 # turns out to be in the final equation. The Chebyshev683 # alternation theorem guarantees that that error in the last684 # equation will be anti-monotonic in the input e, so we can685 # just binary-search until we get the two as close to equal as686 # we need them.687 688 function try_e(e)689 # Try a given value of e, derive the coefficients of the690 # resulting rational function by setting up equations691 # based on the first n+d+1 of the n+d+2 coordinates, and692 # see what the error turns out to be at the final693 # coordinate.694 matrix = array2d(BigFloat, n+d+1, n+d+1)695 vector = array1d(BigFloat, n+d+1)696 currsign = +1697 for i = 1:1:n+d+1698 x = coords[i]699 y = f(x)700 y_adj = y - currsign * e / w(x,y)701 for j = 0:1:n702 matrix[i,1+j] = x^j703 end704 for j = 1:1:d705 matrix[i,1+n+j] = -x^j * y_adj706 end707 vector[i] = y_adj708 currsign = -currsign709 end710 711 @debug("equaldev", "trying e=", e)712 @debug("equaldev", "matrix=", repr(matrix))713 @debug("equaldev", "vector=", repr(vector))714 715 outvector = matrix \ vector716 717 @debug("equaldev", "outvector=", repr(outvector))718 719 ncoeffs = outvector[1:n+1]720 dcoeffs = vcat([BigFloat(1)], outvector[n+2:n+d+1])721 722 x = coords[n+d+2]723 y = f(x)724 last_e = (ratfn_eval(ncoeffs, dcoeffs, x) - y) * w(x,y) * -currsign725 726 @debug("equaldev", "last e=", last_e)727 728 return ncoeffs, dcoeffs, last_e729 end730 731 threshold = 2^(-epsbits/2) # convergence threshold732 733 # Start by trying our previous iteration's error value. This734 # value (e0) will be one end of our binary-search interval,735 # and whatever it caused the last point's error to be, that736 # (e1) will be the other end.737 e0 = prev_err738 @debug("equaldev", "e0 = ", e0)739 nc, dc, e1 = try_e(e0)740 @debug("equaldev", "e1 = ", e1)741 if abs(e1-e0) <= threshold742 # If we're _really_ lucky, we hit the error right on the743 # nose just by doing that!744 return nc, dc745 end746 s = sign(e1-e0)747 @debug("equaldev", "s = ", s)748 749 # Verify by assertion that trying our other interval endpoint750 # e1 gives a value that's wrong in the other direction.751 # (Otherwise our binary search won't get a sensible answer at752 # all.)753 nc, dc, e2 = try_e(e1)754 @debug("equaldev", "e2 = ", e2)755 @assert(sign(e2-e1) == -s)756 757 # Now binary-search until our two endpoints narrow enough.758 local emid759 while abs(e1-e0) > threshold760 emid = (e1+e0)/2761 nc, dc, enew = try_e(emid)762 if sign(enew-emid) == s763 e0 = emid764 else765 e1 = emid766 end767 end768 769 @debug("equaldev", "final e=", emid)770 return nc, dc771 end772end773 774# ----------------------------------------------------------------------775# Top-level function to find a minimax rational-function approximation.776 777# Arguments:778# f The function to be approximated. Maps BigFloat -> BigFloat.779# interval A pair of BigFloats giving the endpoints of the interval780# (in either order) on which to approximate f.781# n, d The degrees of the numerator and denominator of the desired782# approximation.783# w Error-weighting function. Takes two BigFloat arguments x,y784# and returns a scaling factor for the error at that location.785# The returned approximation R should have the minimum possible786# maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional787# parameter, defaulting to the always-return-1 function.788#789# Return values: a tuple (N,D,E,X), where790 791# N,D A pair of arrays of BigFloats giving the coefficients792# of the returned rational function. N has size n+1; D793# has size d+1. Both start with the constant term, i.e.794# N[i] is the coefficient of x^(i-1) (because Julia795# arrays are 1-based). D[1] will be 1.796# E The maximum weighted error (BigFloat).797# X An array of pairs of BigFloats giving the locations of n+2798# points and the weighted error at each of those points. The799# weighted error values will have alternating signs, which800# means that the Chebyshev alternation theorem guarantees801# that any other function of the same degree must exceed802# the error of this one at at least one of those points.803function ratfn_minimax(f::Function, interval::Tuple{BigFloat,BigFloat}, n, d,804 w = (x,y)->BigFloat(1))805 # We start off by finding a least-squares approximation. This806 # doesn't need to be perfect, but if we can get it reasonably good807 # then it'll save iterations in the refining stage.808 #809 # Least-squares approximations tend to look nicer in a minimax810 # sense if you evaluate the function at a big pile of Chebyshev811 # nodes rather than uniformly spaced points. These values will812 # also make a good grid to use for the initial search for error813 # extrema, so we'll keep them around for that reason too.814 815 # Construct the grid.816 lo, hi = minimum(interval), maximum(interval)817 local grid818 let819 mid = (hi+lo)/2820 halfwid = (hi-lo)/2821 nnodes = 16 * (n+d+1)822 pi = 2*asin(BigFloat(1))823 grid = [ mid - halfwid * cos(pi*i/nnodes) for i=0:1:nnodes ]824 end825 826 # Find the initial least-squares approximation.827 (nc, dc) = ratfn_leastsquares(f, grid, n, d, w)828 @debug("minimax", "initial leastsquares approx = ",829 ratfn_to_string(nc, dc))830 831 # Threshold of convergence. We stop when the relative difference832 # between the min and max (winnowed) error extrema is less than833 # this.834 #835 # This is set to the cube root of machine epsilon on a more or836 # less empirical basis, because the rational-function case will837 # not converge reliably if you set it to only the square root.838 # (Repeatable by using the --test mode.) On the assumption that839 # input and output error in each iteration can be expected to be840 # related by a simple power law (because it'll just be down to how841 # many leading terms of a Taylor series are zero), the cube root842 # was the next thing to try.843 threshold = 2^(-epsbits/3)844 845 # Main loop.846 while true847 # Find all the error extrema we can.848 function compute_error(x)849 real_y = f(x)850 approx_y = ratfn_eval(nc, dc, x)851 return (approx_y - real_y) * w(x, real_y)852 end853 extrema = find_extrema(compute_error, grid)854 @debug("minimax", "all extrema = ", format_xylist(extrema))855 856 # Winnow the extrema down to the right number, and ensure they857 # have alternating sign.858 extrema = winnow_extrema(extrema, n+d+2)859 @debug("minimax", "winnowed extrema = ", format_xylist(extrema))860 861 # See if we've finished.862 min_err = minimum([abs(y) for (x,y) = extrema])863 max_err = maximum([abs(y) for (x,y) = extrema])864 variation = (max_err - min_err) / max_err865 @debug("minimax", "extremum variation = ", variation)866 if variation < threshold867 @debug("minimax", "done!")868 return nc, dc, max_err, extrema869 end870 871 # If not, refine our function by equalising the error at the872 # extrema points, and go round again.873 (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema),874 n, d, max_err, w)875 @debug("minimax", "refined approx = ", ratfn_to_string(nc, dc))876 end877end878 879# ----------------------------------------------------------------------880# Check if a polynomial is well-conditioned for accurate evaluation in881# a given interval by Horner's rule.882#883# This is true if at every step where Horner's rule computes884# (coefficient + x*value_so_far), the constant coefficient you're885# adding on is of larger magnitude than the x*value_so_far operand.886# And this has to be true for every x in the interval.887#888# Arguments:889# coeffs The coefficients of the polynomial under test. Starts with890# the constant term, i.e. coeffs[i] is the coefficient of891# x^(i-1) (because Julia arrays are 1-based).892# lo, hi The bounds of the interval.893#894# Return value: the largest ratio (x*value_so_far / coefficient), at895# any step of evaluation, for any x in the interval. If this is less896# than 1, the polynomial is at least somewhat well-conditioned;897# ideally you want it to be more like 1/8 or 1/16 or so, so that the898# relative rounding error accumulated at each step are reduced by899# several factors of 2 when the next coefficient is added on.900 901function wellcond(coeffs, lo, hi)902 x = max(abs(lo), abs(hi))903 worst = 0904 so_far = 0905 for i = length(coeffs):-1:1906 coeff = abs(coeffs[i])907 so_far *= x908 if coeff != 0909 thisval = so_far / coeff910 worst = max(worst, thisval)911 so_far += coeff912 end913 end914 return worst915end916 917# ----------------------------------------------------------------------918# Small set of unit tests.919 920function test()921 passes = 0922 fails = 0923 924 function approx_eq(x, y, limit=1e-6)925 return abs(x - y) < limit926 end927 928 function test(condition)929 if condition930 passes += 1931 else932 println("fail")933 fails += 1934 end935 end936 937 # Test Gaussian elimination.938 println("Gaussian test 1:")939 m = BigFloat[1 1 2; 3 5 8; 13 34 21]940 v = BigFloat[1, -1, 2]941 ret = m \ v942 println(" ",repr(ret))943 test(approx_eq(ret[1], 109/26))944 test(approx_eq(ret[2], -105/130))945 test(approx_eq(ret[3], -31/26))946 947 # Test leastsquares rational functions.948 println("Leastsquares test 1:")949 n = 10000950 a = array1d(BigFloat, n+1)951 for i = 0:1:n952 a[1+i] = i/BigFloat(n)953 end954 (nc, dc) = ratfn_leastsquares(x->exp(x), a, 2, 2)955 println(" ",ratfn_to_string(nc, dc))956 for x = a957 test(approx_eq(exp(x), ratfn_eval(nc, dc, x), 1e-4))958 end959 960 # Test golden section search.961 println("Golden section test 1:")962 x, y = goldensection(x->sin(x),963 BigFloat(0), BigFloat(1)/10, BigFloat(4))964 println(" ", x, " -> ", y)965 test(approx_eq(x, asin(BigFloat(1))))966 test(approx_eq(y, 1))967 968 # Test extrema-winnowing algorithm.969 println("Winnow test 1:")970 extrema = [(x, sin(20*x)*sin(197*x))971 for x in BigFloat(0):BigFloat(1)/1000:BigFloat(1)]972 winnowed = winnow_extrema(extrema, 12)973 println(" ret = ", format_xylist(winnowed))974 prevx, prevy = -1, 0975 for (x,y) = winnowed976 test(x > prevx)977 test(y != 0)978 test(prevy * y <= 0) # tolerates initial prevx having no sign979 test(abs(y) > 0.9)980 prevx, prevy = x, y981 end982 983 # Test actual minimax approximation.984 println("Minimax test 1 (polynomial):")985 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0)986 println(" ",e)987 println(" ",ratfn_to_string(nc, dc))988 test(0 < e < 1e-3)989 for x = 0:BigFloat(1)/1000:1990 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)991 end992 993 println("Minimax test 2 (rational):")994 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2)995 println(" ",e)996 println(" ",ratfn_to_string(nc, dc))997 test(0 < e < 1e-3)998 for x = 0:BigFloat(1)/1000:1999 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)1000 end1001 1002 println("Minimax test 3 (polynomial, weighted):")1003 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0,1004 (x,y)->1/y)1005 println(" ",e)1006 println(" ",ratfn_to_string(nc, dc))1007 test(0 < e < 1e-3)1008 for x = 0:BigFloat(1)/1000:11009 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)1010 end1011 1012 println("Minimax test 4 (rational, weighted):")1013 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2,1014 (x,y)->1/y)1015 println(" ",e)1016 println(" ",ratfn_to_string(nc, dc))1017 test(0 < e < 1e-3)1018 for x = 0:BigFloat(1)/1000:11019 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)1020 end1021 1022 println("Minimax test 5 (rational, weighted, odd degree):")1023 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 1,1024 (x,y)->1/y)1025 println(" ",e)1026 println(" ",ratfn_to_string(nc, dc))1027 test(0 < e < 1e-3)1028 for x = 0:BigFloat(1)/1000:11029 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)1030 end1031 1032 total = passes + fails1033 println(passes, " passes ", fails, " fails ", total, " total")1034end1035 1036# ----------------------------------------------------------------------1037# Online help.1038function help()1039 print("""1040Usage:1041 1042 remez.jl [options] <lo> <hi> <n> <d> <expr> [<weight>]1043 1044Arguments:1045 1046 <lo>, <hi>1047 1048 Bounds of the interval on which to approximate the target1049 function. These are parsed and evaluated as Julia expressions,1050 so you can write things like '1/BigFloat(6)' to get an1051 accurate representation of 1/6, or '4*atan(BigFloat(1))' to1052 get pi. (Unfortunately, the obvious 'BigFloat(pi)' doesn't1053 work in Julia.)1054 1055 <n>, <d>1056 1057 The desired degree of polynomial(s) you want for your1058 approximation. These should be non-negative integers. If you1059 want a rational function as output, set <n> to the degree of1060 the numerator, and <d> the denominator. If you just want an1061 ordinary polynomial, set <d> to 0, and <n> to the degree of1062 the polynomial you want.1063 1064 <expr>1065 1066 A Julia expression giving the function to be approximated on1067 the interval. The input value is predefined as 'x' when this1068 expression is evaluated, so you should write something along1069 the lines of 'sin(x)' or 'sqrt(1+tan(x)^2)' etc.1070 1071 <weight>1072 1073 If provided, a Julia expression giving the weighting factor1074 for the approximation error. The output polynomial will1075 minimise the largest absolute value of (P-f) * w at any point1076 in the interval, where P is the value of the polynomial, f is1077 the value of the target function given by <expr>, and w is the1078 weight given by this function.1079 1080 When this expression is evaluated, the input value to P and f1081 is predefined as 'x', and also the true output value f(x) is1082 predefined as 'y'. So you can minimise the relative error by1083 simply writing '1/y'.1084 1085 If the <weight> argument is not provided, the default1086 weighting function always returns 1, so that the polynomial1087 will minimise the maximum absolute error |P-f|.1088 1089Computation options:1090 1091 --pre=<predef_expr>1092 1093 Evaluate the Julia expression <predef_expr> before starting1094 the computation. This permits you to pre-define variables or1095 functions which the Julia expressions in your main arguments1096 can refer to. All of <lo>, <hi>, <expr> and <weight> can make1097 use of things defined by <predef_expr>.1098 1099 One internal remez.jl function that you might sometimes find1100 useful in this expression is 'goldensection', which finds the1101 location and value of a maximum of a function. For example,1102 one implementation strategy for the gamma function involves1103 translating it to put its unique local minimum at the origin,1104 in which case you can write something like this1105 1106 --pre='(m,my) = goldensection(x -> -gamma(x),1107 BigFloat(1), BigFloat(1.5), BigFloat(2))'1108 1109 to predefine 'm' as the location of gamma's minimum, and 'my'1110 as the (negated) value that gamma actually takes at that1111 point, i.e. -gamma(m).1112 1113 (Since 'goldensection' always finds a maximum, we had to1114 negate gamma in the input function to make it find a minimum1115 instead. Consult the comments in the source for more details1116 on the use of this function.)1117 1118 If you use this option more than once, all the expressions you1119 provide will be run in sequence.1120 1121 --bits=<bits>1122 1123 Specify the accuracy to which you want the output polynomial,1124 in bits. Default 256, which should be more than enough.1125 1126 --bigfloatbits=<bits>1127 1128 Turn up the precision used by Julia for its BigFloat1129 evaluation. Default is Julia's default (also 256). You might1130 want to try setting this higher than the --bits value if the1131 algorithm is failing to converge for some reason.1132 1133Output options:1134 1135 --full1136 1137 Instead of just printing the approximation function itself,1138 also print auxiliary information:1139 - the locations of the error extrema, and the actual1140 (weighted) error at each of those locations1141 - the overall maximum error of the function1142 - a 'well-conditioning quotient', giving the worst-case ratio1143 between any polynomial coefficient and the largest possible1144 value of the higher-order terms it will be added to.1145 1146 The well-conditioning quotient should be less than 1, ideally1147 by several factors of two, for accurate evaluation in the1148 target precision. If you request a rational function, a1149 separate well-conditioning quotient will be printed for the1150 numerator and denominator.1151 1152 Use this option when deciding how wide an interval to1153 approximate your function on, and what degree of polynomial1154 you need.1155 1156 --variable=<identifier>1157 1158 When writing the output polynomial or rational function in its1159 usual form as an arithmetic expression, use <identifier> as1160 the name of the input variable. Default is 'x'.1161 1162 --suffix=<suffix>1163 1164 When writing the output polynomial or rational function in its1165 usual form as an arithmetic expression, write <suffix> after1166 every floating-point literal. For example, '--suffix=F' will1167 generate a C expression in which the coefficients are literals1168 of type 'float' rather than 'double'.1169 1170 --array1171 1172 Instead of writing the output polynomial as an arithmetic1173 expression in Horner's rule form, write out just its1174 coefficients, one per line, each with a trailing comma.1175 Suitable for pasting into a C array declaration.1176 1177 This option is not currently supported if the output is a1178 rational function, because you'd need two separate arrays for1179 the numerator and denominator coefficients and there's no1180 obviously right way to provide both of those together.1181 1182Debug and test options:1183 1184 --debug=<facility>1185 1186 Enable debugging output from various parts of the Remez1187 calculation. <facility> should be the name of one of the1188 classes of diagnostic output implemented in the program.1189 Useful values include 'gausselim', 'leastsquares',1190 'goldensection', 'equaldev', 'minimax'. This is probably1191 mostly useful to people debugging problems with the script, so1192 consult the source code for more information about what the1193 diagnostic output for each of those facilities will be.1194 1195 If you want diagnostics from more than one facility, specify1196 this option multiple times with different arguments.1197 1198 --test1199 1200 Run remez.jl's internal test suite. No arguments needed.1201 1202Miscellaneous options:1203 1204 --help1205 1206 Display this text and exit. No arguments needed.1207 1208""")1209end1210 1211# ----------------------------------------------------------------------1212# Main program.1213 1214function main()1215 nargs = length(argwords)1216 if nargs != 5 && nargs != 61217 error("usage: remez.jl <lo> <hi> <n> <d> <expr> [<weight>]\n" *1218 " run 'remez.jl --help' for more help")1219 end1220 1221 for preliminary_command in preliminary_commands1222 eval(Meta.parse(preliminary_command))1223 end1224 1225 lo = BigFloat(eval(Meta.parse(argwords[1])))1226 hi = BigFloat(eval(Meta.parse(argwords[2])))1227 n = parse(Int,argwords[3])1228 d = parse(Int,argwords[4])1229 f = eval(Meta.parse("x -> " * argwords[5]))1230 1231 # Wrap the user-provided function with a function of our own. This1232 # arranges to detect silly FP values (inf,nan) early and diagnose1233 # them sensibly, and also lets us log all evaluations of the1234 # function in case you suspect it's doing the wrong thing at some1235 # special-case point.1236 function func(x)1237 y = run(f,x)1238 @debug("f", x, " -> ", y)1239 if !isfinite(y)1240 error("f(" * string(x) * ") returned non-finite value " * string(y))1241 end1242 return y1243 end1244 1245 if nargs == 61246 # Wrap the user-provided weight function similarly.1247 w = eval(Meta.parse("(x,y) -> " * argwords[6]))1248 function wrapped_weight(x,y)1249 ww = run(w,x,y)1250 if !isfinite(ww)1251 error("w(" * string(x) * "," * string(y) *1252 ") returned non-finite value " * string(ww))1253 end1254 return ww1255 end1256 weight = wrapped_weight1257 else1258 weight = (x,y)->BigFloat(1)1259 end1260 1261 (nc, dc, e, extrema) = ratfn_minimax(func, (lo, hi), n, d, weight)1262 if array_format1263 if d == 01264 functext = join([string(x)*",\n" for x=nc],"")1265 else1266 # It's unclear how you should best format an array of1267 # coefficients for a rational function, so I'll leave1268 # implementing this option until I have a use case.1269 error("--array unsupported for rational functions")1270 end1271 else1272 functext = ratfn_to_string(nc, dc) * "\n"1273 end1274 if full_output1275 # Print everything you might want to know about the function1276 println("extrema = ", format_xylist(extrema))1277 println("maxerror = ", string(e))1278 if length(dc) > 11279 println("wellconditioning_numerator = ",1280 string(wellcond(nc, lo, hi)))1281 println("wellconditioning_denominator = ",1282 string(wellcond(dc, lo, hi)))1283 else1284 println("wellconditioning = ", string(wellcond(nc, lo, hi)))1285 end1286 print("function = ", functext)1287 else1288 # Just print the text people will want to paste into their code1289 print(functext)1290 end1291end1292 1293# ----------------------------------------------------------------------1294# Top-level code: parse the argument list and decide what to do.1295 1296what_to_do = main1297 1298doing_opts = true1299argwords = array1d(String, 0)1300for arg = ARGS1301 global doing_opts, what_to_do, argwords1302 global full_output, array_format, xvarname, floatsuffix, epsbits1303 if doing_opts && startswith(arg, "-")1304 if arg == "--"1305 doing_opts = false1306 elseif arg == "--help"1307 what_to_do = help1308 elseif arg == "--test"1309 what_to_do = test1310 elseif arg == "--full"1311 full_output = true1312 elseif arg == "--array"1313 array_format = true1314 elseif startswith(arg, "--debug=")1315 enable_debug(arg[length("--debug=")+1:end])1316 elseif startswith(arg, "--variable=")1317 xvarname = arg[length("--variable=")+1:end]1318 elseif startswith(arg, "--suffix=")1319 floatsuffix = arg[length("--suffix=")+1:end]1320 elseif startswith(arg, "--bits=")1321 epsbits = parse(Int,arg[length("--bits=")+1:end])1322 elseif startswith(arg, "--bigfloatbits=")1323 set_bigfloat_precision(1324 parse(Int,arg[length("--bigfloatbits=")+1:end]))1325 elseif startswith(arg, "--pre=")1326 push!(preliminary_commands, arg[length("--pre=")+1:end])1327 else1328 error("unrecognised option: ", arg)1329 end1330 else1331 push!(argwords, arg)1332 end1333end1334 1335what_to_do()1336