1bbbf1280Sopenharmony_ci#!/usr/bin/env julia 2bbbf1280Sopenharmony_ci# -*- julia -*- 3bbbf1280Sopenharmony_ci 4bbbf1280Sopenharmony_ci# remez.jl - implementation of the Remez algorithm for polynomial approximation 5bbbf1280Sopenharmony_ci# 6bbbf1280Sopenharmony_ci# Copyright (c) 2015-2019, Arm Limited. 7bbbf1280Sopenharmony_ci# SPDX-License-Identifier: MIT 8bbbf1280Sopenharmony_ci 9bbbf1280Sopenharmony_ciimport Base.\ 10bbbf1280Sopenharmony_ci 11bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 12bbbf1280Sopenharmony_ci# Helper functions to cope with different Julia versions. 13bbbf1280Sopenharmony_ciif VERSION >= v"0.7.0" 14bbbf1280Sopenharmony_ci array1d(T, d) = Array{T, 1}(undef, d) 15bbbf1280Sopenharmony_ci array2d(T, d1, d2) = Array{T, 2}(undef, d1, d2) 16bbbf1280Sopenharmony_cielse 17bbbf1280Sopenharmony_ci array1d(T, d) = Array(T, d) 18bbbf1280Sopenharmony_ci array2d(T, d1, d2) = Array(T, d1, d2) 19bbbf1280Sopenharmony_ciend 20bbbf1280Sopenharmony_ciif VERSION < v"0.5.0" 21bbbf1280Sopenharmony_ci String = ASCIIString 22bbbf1280Sopenharmony_ciend 23bbbf1280Sopenharmony_ciif VERSION >= v"0.6.0" 24bbbf1280Sopenharmony_ci # Use Base.invokelatest to run functions made using eval(), to 25bbbf1280Sopenharmony_ci # avoid "world age" error 26bbbf1280Sopenharmony_ci run(f, x...) = Base.invokelatest(f, x...) 27bbbf1280Sopenharmony_cielse 28bbbf1280Sopenharmony_ci # Prior to 0.6.0, invokelatest doesn't exist (but fortunately the 29bbbf1280Sopenharmony_ci # world age problem also doesn't seem to exist) 30bbbf1280Sopenharmony_ci run(f, x...) = f(x...) 31bbbf1280Sopenharmony_ciend 32bbbf1280Sopenharmony_ci 33bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 34bbbf1280Sopenharmony_ci# Global variables configured by command-line options. 35bbbf1280Sopenharmony_cifloatsuffix = "" # adjusted by --floatsuffix 36bbbf1280Sopenharmony_cixvarname = "x" # adjusted by --variable 37bbbf1280Sopenharmony_ciepsbits = 256 # adjusted by --bits 38bbbf1280Sopenharmony_cidebug_facilities = Set() # adjusted by --debug 39bbbf1280Sopenharmony_cifull_output = false # adjusted by --full 40bbbf1280Sopenharmony_ciarray_format = false # adjusted by --array 41bbbf1280Sopenharmony_cipreliminary_commands = array1d(String, 0) # adjusted by --pre 42bbbf1280Sopenharmony_ci 43bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 44bbbf1280Sopenharmony_ci# Diagnostic and utility functions. 45bbbf1280Sopenharmony_ci 46bbbf1280Sopenharmony_ci# Enable debugging printouts from a particular subpart of this 47bbbf1280Sopenharmony_ci# program. 48bbbf1280Sopenharmony_ci# 49bbbf1280Sopenharmony_ci# Arguments: 50bbbf1280Sopenharmony_ci# facility Name of the facility to debug. For a list of facility names, 51bbbf1280Sopenharmony_ci# look through the code for calls to debug(). 52bbbf1280Sopenharmony_ci# 53bbbf1280Sopenharmony_ci# Return value is a BigFloat. 54bbbf1280Sopenharmony_cifunction enable_debug(facility) 55bbbf1280Sopenharmony_ci push!(debug_facilities, facility) 56bbbf1280Sopenharmony_ciend 57bbbf1280Sopenharmony_ci 58bbbf1280Sopenharmony_ci# Print a diagnostic. 59bbbf1280Sopenharmony_ci# 60bbbf1280Sopenharmony_ci# Arguments: 61bbbf1280Sopenharmony_ci# facility Name of the facility for which this is a debug message. 62bbbf1280Sopenharmony_ci# printargs Arguments to println() if debugging of that facility is 63bbbf1280Sopenharmony_ci# enabled. 64bbbf1280Sopenharmony_cimacro debug(facility, printargs...) 65bbbf1280Sopenharmony_ci printit = quote 66bbbf1280Sopenharmony_ci print("[", $facility, "] ") 67bbbf1280Sopenharmony_ci end 68bbbf1280Sopenharmony_ci for arg in printargs 69bbbf1280Sopenharmony_ci printit = quote 70bbbf1280Sopenharmony_ci $printit 71bbbf1280Sopenharmony_ci print($(esc(arg))) 72bbbf1280Sopenharmony_ci end 73bbbf1280Sopenharmony_ci end 74bbbf1280Sopenharmony_ci return quote 75bbbf1280Sopenharmony_ci if $facility in debug_facilities 76bbbf1280Sopenharmony_ci $printit 77bbbf1280Sopenharmony_ci println() 78bbbf1280Sopenharmony_ci end 79bbbf1280Sopenharmony_ci end 80bbbf1280Sopenharmony_ciend 81bbbf1280Sopenharmony_ci 82bbbf1280Sopenharmony_ci# Evaluate a polynomial. 83bbbf1280Sopenharmony_ci 84bbbf1280Sopenharmony_ci# Arguments: 85bbbf1280Sopenharmony_ci# coeffs Array of BigFloats giving the coefficients of the polynomial. 86bbbf1280Sopenharmony_ci# Starts with the constant term, i.e. coeffs[i] is the 87bbbf1280Sopenharmony_ci# coefficient of x^(i-1) (because Julia arrays are 1-based). 88bbbf1280Sopenharmony_ci# x Point at which to evaluate the polynomial. 89bbbf1280Sopenharmony_ci# 90bbbf1280Sopenharmony_ci# Return value is a BigFloat. 91bbbf1280Sopenharmony_cifunction poly_eval(coeffs::Array{BigFloat}, x::BigFloat) 92bbbf1280Sopenharmony_ci n = length(coeffs) 93bbbf1280Sopenharmony_ci if n == 0 94bbbf1280Sopenharmony_ci return BigFloat(0) 95bbbf1280Sopenharmony_ci elseif n == 1 96bbbf1280Sopenharmony_ci return coeffs[1] 97bbbf1280Sopenharmony_ci else 98bbbf1280Sopenharmony_ci return coeffs[1] + x * poly_eval(coeffs[2:n], x) 99bbbf1280Sopenharmony_ci end 100bbbf1280Sopenharmony_ciend 101bbbf1280Sopenharmony_ci 102bbbf1280Sopenharmony_ci# Evaluate a rational function. 103bbbf1280Sopenharmony_ci 104bbbf1280Sopenharmony_ci# Arguments: 105bbbf1280Sopenharmony_ci# ncoeffs Array of BigFloats giving the coefficients of the numerator. 106bbbf1280Sopenharmony_ci# Starts with the constant term, and 1-based, as above. 107bbbf1280Sopenharmony_ci# dcoeffs Array of BigFloats giving the coefficients of the denominator. 108bbbf1280Sopenharmony_ci# Starts with the constant term, and 1-based, as above. 109bbbf1280Sopenharmony_ci# x Point at which to evaluate the function. 110bbbf1280Sopenharmony_ci# 111bbbf1280Sopenharmony_ci# Return value is a BigFloat. 112bbbf1280Sopenharmony_cifunction ratfn_eval(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat}, 113bbbf1280Sopenharmony_ci x::BigFloat) 114bbbf1280Sopenharmony_ci return poly_eval(ncoeffs, x) / poly_eval(dcoeffs, x) 115bbbf1280Sopenharmony_ciend 116bbbf1280Sopenharmony_ci 117bbbf1280Sopenharmony_ci# Format a BigFloat into an appropriate output format. 118bbbf1280Sopenharmony_ci# Arguments: 119bbbf1280Sopenharmony_ci# x BigFloat to format. 120bbbf1280Sopenharmony_ci# 121bbbf1280Sopenharmony_ci# Return value is a string. 122bbbf1280Sopenharmony_cifunction float_to_str(x) 123bbbf1280Sopenharmony_ci return string(x) * floatsuffix 124bbbf1280Sopenharmony_ciend 125bbbf1280Sopenharmony_ci 126bbbf1280Sopenharmony_ci# Format a polynomial into an arithmetic expression, for pasting into 127bbbf1280Sopenharmony_ci# other tools such as gnuplot. 128bbbf1280Sopenharmony_ci 129bbbf1280Sopenharmony_ci# Arguments: 130bbbf1280Sopenharmony_ci# coeffs Array of BigFloats giving the coefficients of the polynomial. 131bbbf1280Sopenharmony_ci# Starts with the constant term, and 1-based, as above. 132bbbf1280Sopenharmony_ci# 133bbbf1280Sopenharmony_ci# Return value is a string. 134bbbf1280Sopenharmony_cifunction poly_to_string(coeffs::Array{BigFloat}) 135bbbf1280Sopenharmony_ci n = length(coeffs) 136bbbf1280Sopenharmony_ci if n == 0 137bbbf1280Sopenharmony_ci return "0" 138bbbf1280Sopenharmony_ci elseif n == 1 139bbbf1280Sopenharmony_ci return float_to_str(coeffs[1]) 140bbbf1280Sopenharmony_ci else 141bbbf1280Sopenharmony_ci return string(float_to_str(coeffs[1]), "+", xvarname, "*(", 142bbbf1280Sopenharmony_ci poly_to_string(coeffs[2:n]), ")") 143bbbf1280Sopenharmony_ci end 144bbbf1280Sopenharmony_ciend 145bbbf1280Sopenharmony_ci 146bbbf1280Sopenharmony_ci# Format a rational function into a string. 147bbbf1280Sopenharmony_ci 148bbbf1280Sopenharmony_ci# Arguments: 149bbbf1280Sopenharmony_ci# ncoeffs Array of BigFloats giving the coefficients of the numerator. 150bbbf1280Sopenharmony_ci# Starts with the constant term, and 1-based, as above. 151bbbf1280Sopenharmony_ci# dcoeffs Array of BigFloats giving the coefficients of the denominator. 152bbbf1280Sopenharmony_ci# Starts with the constant term, and 1-based, as above. 153bbbf1280Sopenharmony_ci# 154bbbf1280Sopenharmony_ci# Return value is a string. 155bbbf1280Sopenharmony_cifunction ratfn_to_string(ncoeffs::Array{BigFloat}, dcoeffs::Array{BigFloat}) 156bbbf1280Sopenharmony_ci if length(dcoeffs) == 1 && dcoeffs[1] == 1 157bbbf1280Sopenharmony_ci # Special case: if the denominator is just 1, leave it out. 158bbbf1280Sopenharmony_ci return poly_to_string(ncoeffs) 159bbbf1280Sopenharmony_ci else 160bbbf1280Sopenharmony_ci return string("(", poly_to_string(ncoeffs), ")/(", 161bbbf1280Sopenharmony_ci poly_to_string(dcoeffs), ")") 162bbbf1280Sopenharmony_ci end 163bbbf1280Sopenharmony_ciend 164bbbf1280Sopenharmony_ci 165bbbf1280Sopenharmony_ci# Format a list of x,y pairs into a string. 166bbbf1280Sopenharmony_ci 167bbbf1280Sopenharmony_ci# Arguments: 168bbbf1280Sopenharmony_ci# xys Array of (x,y) pairs of BigFloats. 169bbbf1280Sopenharmony_ci# 170bbbf1280Sopenharmony_ci# Return value is a string. 171bbbf1280Sopenharmony_cifunction format_xylist(xys::Array{Tuple{BigFloat,BigFloat}}) 172bbbf1280Sopenharmony_ci return ("[\n" * 173bbbf1280Sopenharmony_ci join([" "*string(x)*" -> "*string(y) for (x,y) in xys], "\n") * 174bbbf1280Sopenharmony_ci "\n]") 175bbbf1280Sopenharmony_ciend 176bbbf1280Sopenharmony_ci 177bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 178bbbf1280Sopenharmony_ci# Matrix-equation solver for matrices of BigFloat. 179bbbf1280Sopenharmony_ci# 180bbbf1280Sopenharmony_ci# I had hoped that Julia's type-genericity would allow me to solve the 181bbbf1280Sopenharmony_ci# matrix equation Mx=V by just writing 'M \ V'. Unfortunately, that 182bbbf1280Sopenharmony_ci# works by translating the inputs into double precision and handing 183bbbf1280Sopenharmony_ci# off to an optimised library, which misses the point when I have a 184bbbf1280Sopenharmony_ci# matrix and vector of BigFloat and want my result in _better_ than 185bbbf1280Sopenharmony_ci# double precision. So I have to implement my own specialisation of 186bbbf1280Sopenharmony_ci# the \ operator for that case. 187bbbf1280Sopenharmony_ci# 188bbbf1280Sopenharmony_ci# Fortunately, the point of using BigFloats is that we have precision 189bbbf1280Sopenharmony_ci# to burn, so I can do completely naïve Gaussian elimination without 190bbbf1280Sopenharmony_ci# worrying about instability. 191bbbf1280Sopenharmony_ci 192bbbf1280Sopenharmony_ci# Arguments: 193bbbf1280Sopenharmony_ci# matrix_in 2-dimensional array of BigFloats, representing a matrix M 194bbbf1280Sopenharmony_ci# in row-first order, i.e. matrix_in[r,c] represents the 195bbbf1280Sopenharmony_ci# entry in row r col c. 196bbbf1280Sopenharmony_ci# vector_in 1-dimensional array of BigFloats, representing a vector V. 197bbbf1280Sopenharmony_ci# 198bbbf1280Sopenharmony_ci# Return value: a 1-dimensional array X of BigFloats, satisfying M X = V. 199bbbf1280Sopenharmony_ci# 200bbbf1280Sopenharmony_ci# Expects the input to be an invertible square matrix and a vector of 201bbbf1280Sopenharmony_ci# the corresponding size, on pain of failing an assertion. 202bbbf1280Sopenharmony_cifunction \(matrix_in :: Array{BigFloat,2}, 203bbbf1280Sopenharmony_ci vector_in :: Array{BigFloat,1}) 204bbbf1280Sopenharmony_ci # Copy the inputs, because we'll be mutating them as we go. 205bbbf1280Sopenharmony_ci M = copy(matrix_in) 206bbbf1280Sopenharmony_ci V = copy(vector_in) 207bbbf1280Sopenharmony_ci 208bbbf1280Sopenharmony_ci # Input consistency criteria: matrix is square, and vector has 209bbbf1280Sopenharmony_ci # length to match. 210bbbf1280Sopenharmony_ci n = length(V) 211bbbf1280Sopenharmony_ci @assert(n > 0) 212bbbf1280Sopenharmony_ci @assert(size(M) == (n,n)) 213bbbf1280Sopenharmony_ci 214bbbf1280Sopenharmony_ci @debug("gausselim", "starting, n=", n) 215bbbf1280Sopenharmony_ci 216bbbf1280Sopenharmony_ci for i = 1:1:n 217bbbf1280Sopenharmony_ci # Straightforward Gaussian elimination: find the largest 218bbbf1280Sopenharmony_ci # non-zero entry in column i (and in a row we haven't sorted 219bbbf1280Sopenharmony_ci # out already), swap it into row i, scale that row to 220bbbf1280Sopenharmony_ci # normalise it to 1, then zero out the rest of the column by 221bbbf1280Sopenharmony_ci # subtracting a multiple of that row from each other row. 222bbbf1280Sopenharmony_ci 223bbbf1280Sopenharmony_ci @debug("gausselim", "matrix=", repr(M)) 224bbbf1280Sopenharmony_ci @debug("gausselim", "vector=", repr(V)) 225bbbf1280Sopenharmony_ci 226bbbf1280Sopenharmony_ci # Find the best pivot. 227bbbf1280Sopenharmony_ci bestrow = 0 228bbbf1280Sopenharmony_ci bestval = 0 229bbbf1280Sopenharmony_ci for j = i:1:n 230bbbf1280Sopenharmony_ci if abs(M[j,i]) > bestval 231bbbf1280Sopenharmony_ci bestrow = j 232bbbf1280Sopenharmony_ci bestval = M[j,i] 233bbbf1280Sopenharmony_ci end 234bbbf1280Sopenharmony_ci end 235bbbf1280Sopenharmony_ci @assert(bestrow > 0) # make sure we did actually find one 236bbbf1280Sopenharmony_ci 237bbbf1280Sopenharmony_ci @debug("gausselim", "bestrow=", bestrow) 238bbbf1280Sopenharmony_ci 239bbbf1280Sopenharmony_ci # Swap it into row i. 240bbbf1280Sopenharmony_ci if bestrow != i 241bbbf1280Sopenharmony_ci for k = 1:1:n 242bbbf1280Sopenharmony_ci M[bestrow,k],M[i,k] = M[i,k],M[bestrow,k] 243bbbf1280Sopenharmony_ci end 244bbbf1280Sopenharmony_ci V[bestrow],V[i] = V[i],V[bestrow] 245bbbf1280Sopenharmony_ci end 246bbbf1280Sopenharmony_ci 247bbbf1280Sopenharmony_ci # Scale that row so that M[i,i] becomes 1. 248bbbf1280Sopenharmony_ci divisor = M[i,i] 249bbbf1280Sopenharmony_ci for k = 1:1:n 250bbbf1280Sopenharmony_ci M[i,k] = M[i,k] / divisor 251bbbf1280Sopenharmony_ci end 252bbbf1280Sopenharmony_ci V[i] = V[i] / divisor 253bbbf1280Sopenharmony_ci @assert(M[i,i] == 1) 254bbbf1280Sopenharmony_ci 255bbbf1280Sopenharmony_ci # Zero out all other entries in column i, by subtracting 256bbbf1280Sopenharmony_ci # multiples of this row. 257bbbf1280Sopenharmony_ci for j = 1:1:n 258bbbf1280Sopenharmony_ci if j != i 259bbbf1280Sopenharmony_ci factor = M[j,i] 260bbbf1280Sopenharmony_ci for k = 1:1:n 261bbbf1280Sopenharmony_ci M[j,k] = M[j,k] - M[i,k] * factor 262bbbf1280Sopenharmony_ci end 263bbbf1280Sopenharmony_ci V[j] = V[j] - V[i] * factor 264bbbf1280Sopenharmony_ci @assert(M[j,i] == 0) 265bbbf1280Sopenharmony_ci end 266bbbf1280Sopenharmony_ci end 267bbbf1280Sopenharmony_ci end 268bbbf1280Sopenharmony_ci 269bbbf1280Sopenharmony_ci @debug("gausselim", "matrix=", repr(M)) 270bbbf1280Sopenharmony_ci @debug("gausselim", "vector=", repr(V)) 271bbbf1280Sopenharmony_ci @debug("gausselim", "done!") 272bbbf1280Sopenharmony_ci 273bbbf1280Sopenharmony_ci # Now we're done: M is the identity matrix, so the equation Mx=V 274bbbf1280Sopenharmony_ci # becomes just x=V, i.e. V is already exactly the vector we want 275bbbf1280Sopenharmony_ci # to return. 276bbbf1280Sopenharmony_ci return V 277bbbf1280Sopenharmony_ciend 278bbbf1280Sopenharmony_ci 279bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 280bbbf1280Sopenharmony_ci# Least-squares fitting of a rational function to a set of (x,y) 281bbbf1280Sopenharmony_ci# points. 282bbbf1280Sopenharmony_ci# 283bbbf1280Sopenharmony_ci# We use this to get an initial starting point for the Remez 284bbbf1280Sopenharmony_ci# iteration. Therefore, it doesn't really need to be particularly 285bbbf1280Sopenharmony_ci# accurate; it only needs to be good enough to wiggle back and forth 286bbbf1280Sopenharmony_ci# across the target function the right number of times (so as to give 287bbbf1280Sopenharmony_ci# enough error extrema to start optimising from) and not have any 288bbbf1280Sopenharmony_ci# poles in the target interval. 289bbbf1280Sopenharmony_ci# 290bbbf1280Sopenharmony_ci# Least-squares fitting of a _polynomial_ is actually a sensible thing 291bbbf1280Sopenharmony_ci# to do, and minimises the rms error. Doing the following trick with a 292bbbf1280Sopenharmony_ci# rational function P/Q is less sensible, because it cannot be made to 293bbbf1280Sopenharmony_ci# minimise the error function (P/Q-f)^2 that you actually wanted; 294bbbf1280Sopenharmony_ci# instead it minimises (P-fQ)^2. But that should be good enough to 295bbbf1280Sopenharmony_ci# have the properties described above. 296bbbf1280Sopenharmony_ci# 297bbbf1280Sopenharmony_ci# Some theory: suppose you're trying to choose a set of parameters a_i 298bbbf1280Sopenharmony_ci# so as to minimise the sum of squares of some error function E_i. 299bbbf1280Sopenharmony_ci# Basic calculus says, if you do this in one variable, just 300bbbf1280Sopenharmony_ci# differentiate and solve for zero. In this case, that works fine even 301bbbf1280Sopenharmony_ci# with multiple variables, because you _partially_ differentiate with 302bbbf1280Sopenharmony_ci# respect to each a_i, giving a system of equations, and that system 303bbbf1280Sopenharmony_ci# turns out to be linear so we just solve it as a matrix. 304bbbf1280Sopenharmony_ci# 305bbbf1280Sopenharmony_ci# In this case, our parameters are the coefficients of P and Q; to 306bbbf1280Sopenharmony_ci# avoid underdetermining the system we'll fix Q's constant term at 1, 307bbbf1280Sopenharmony_ci# so that our error function (as described above) is 308bbbf1280Sopenharmony_ci# 309bbbf1280Sopenharmony_ci# E = \sum (p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d)^2 310bbbf1280Sopenharmony_ci# 311bbbf1280Sopenharmony_ci# where the sum is over all (x,y) coordinate pairs. Setting dE/dp_j=0 312bbbf1280Sopenharmony_ci# (for each j) gives an equation of the form 313bbbf1280Sopenharmony_ci# 314bbbf1280Sopenharmony_ci# 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) x^j 315bbbf1280Sopenharmony_ci# 316bbbf1280Sopenharmony_ci# and setting dE/dq_j=0 gives one of the form 317bbbf1280Sopenharmony_ci# 318bbbf1280Sopenharmony_ci# 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 319bbbf1280Sopenharmony_ci# 320bbbf1280Sopenharmony_ci# And both of those row types, treated as multivariate linear 321bbbf1280Sopenharmony_ci# equations in the p,q values, have each coefficient being a value of 322bbbf1280Sopenharmony_ci# the form \sum x^i, \sum y x^i or \sum y^2 x^i, for various i. (Times 323bbbf1280Sopenharmony_ci# a factor of 2, but we can throw that away.) So we can go through the 324bbbf1280Sopenharmony_ci# list of input coordinates summing all of those things, and then we 325bbbf1280Sopenharmony_ci# have enough information to construct our matrix and solve it 326bbbf1280Sopenharmony_ci# straight off for the rational function coefficients. 327bbbf1280Sopenharmony_ci 328bbbf1280Sopenharmony_ci# Arguments: 329bbbf1280Sopenharmony_ci# f The function to be approximated. Maps BigFloat -> BigFloat. 330bbbf1280Sopenharmony_ci# xvals Array of BigFloats, giving the list of x-coordinates at which 331bbbf1280Sopenharmony_ci# to evaluate f. 332bbbf1280Sopenharmony_ci# n Degree of the numerator polynomial of the desired rational 333bbbf1280Sopenharmony_ci# function. 334bbbf1280Sopenharmony_ci# d Degree of the denominator polynomial of the desired rational 335bbbf1280Sopenharmony_ci# function. 336bbbf1280Sopenharmony_ci# w Error-weighting function. Takes two BigFloat arguments x,y 337bbbf1280Sopenharmony_ci# and returns a scaling factor for the error at that location. 338bbbf1280Sopenharmony_ci# A larger value indicates that the error should be given 339bbbf1280Sopenharmony_ci# greater weight in the square sum we try to minimise. 340bbbf1280Sopenharmony_ci# If unspecified, defaults to giving everything the same weight. 341bbbf1280Sopenharmony_ci# 342bbbf1280Sopenharmony_ci# Return values: a pair of arrays of BigFloats (N,D) giving the 343bbbf1280Sopenharmony_ci# coefficients of the returned rational function. N has size n+1; D 344bbbf1280Sopenharmony_ci# has size d+1. Both start with the constant term, i.e. N[i] is the 345bbbf1280Sopenharmony_ci# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will 346bbbf1280Sopenharmony_ci# be 1. 347bbbf1280Sopenharmony_cifunction ratfn_leastsquares(f::Function, xvals::Array{BigFloat}, n, d, 348bbbf1280Sopenharmony_ci w = (x,y)->BigFloat(1)) 349bbbf1280Sopenharmony_ci # Accumulate sums of x^i y^j, for j={0,1,2} and a range of x. 350bbbf1280Sopenharmony_ci # Again because Julia arrays are 1-based, we'll have sums[i,j] 351bbbf1280Sopenharmony_ci # being the sum of x^(i-1) y^(j-1). 352bbbf1280Sopenharmony_ci maxpow = max(n,d) * 2 + 1 353bbbf1280Sopenharmony_ci sums = zeros(BigFloat, maxpow, 3) 354bbbf1280Sopenharmony_ci for x = xvals 355bbbf1280Sopenharmony_ci y = f(x) 356bbbf1280Sopenharmony_ci weight = w(x,y) 357bbbf1280Sopenharmony_ci for i = 1:1:maxpow 358bbbf1280Sopenharmony_ci for j = 1:1:3 359bbbf1280Sopenharmony_ci sums[i,j] += x^(i-1) * y^(j-1) * weight 360bbbf1280Sopenharmony_ci end 361bbbf1280Sopenharmony_ci end 362bbbf1280Sopenharmony_ci end 363bbbf1280Sopenharmony_ci 364bbbf1280Sopenharmony_ci @debug("leastsquares", "sums=", repr(sums)) 365bbbf1280Sopenharmony_ci 366bbbf1280Sopenharmony_ci # Build the matrix. We're solving n+d+1 equations in n+d+1 367bbbf1280Sopenharmony_ci # unknowns. (We actually have to return n+d+2 coefficients, but 368bbbf1280Sopenharmony_ci # one of them is hardwired to 1.) 369bbbf1280Sopenharmony_ci matrix = array2d(BigFloat, n+d+1, n+d+1) 370bbbf1280Sopenharmony_ci vector = array1d(BigFloat, n+d+1) 371bbbf1280Sopenharmony_ci for i = 0:1:n 372bbbf1280Sopenharmony_ci # Equation obtained by differentiating with respect to p_i, 373bbbf1280Sopenharmony_ci # i.e. the numerator coefficient of x^i. 374bbbf1280Sopenharmony_ci row = 1+i 375bbbf1280Sopenharmony_ci for j = 0:1:n 376bbbf1280Sopenharmony_ci matrix[row, 1+j] = sums[1+i+j, 1] 377bbbf1280Sopenharmony_ci end 378bbbf1280Sopenharmony_ci for j = 1:1:d 379bbbf1280Sopenharmony_ci matrix[row, 1+n+j] = -sums[1+i+j, 2] 380bbbf1280Sopenharmony_ci end 381bbbf1280Sopenharmony_ci vector[row] = sums[1+i, 2] 382bbbf1280Sopenharmony_ci end 383bbbf1280Sopenharmony_ci for i = 1:1:d 384bbbf1280Sopenharmony_ci # Equation obtained by differentiating with respect to q_i, 385bbbf1280Sopenharmony_ci # i.e. the denominator coefficient of x^i. 386bbbf1280Sopenharmony_ci row = 1+n+i 387bbbf1280Sopenharmony_ci for j = 0:1:n 388bbbf1280Sopenharmony_ci matrix[row, 1+j] = sums[1+i+j, 2] 389bbbf1280Sopenharmony_ci end 390bbbf1280Sopenharmony_ci for j = 1:1:d 391bbbf1280Sopenharmony_ci matrix[row, 1+n+j] = -sums[1+i+j, 3] 392bbbf1280Sopenharmony_ci end 393bbbf1280Sopenharmony_ci vector[row] = sums[1+i, 3] 394bbbf1280Sopenharmony_ci end 395bbbf1280Sopenharmony_ci 396bbbf1280Sopenharmony_ci @debug("leastsquares", "matrix=", repr(matrix)) 397bbbf1280Sopenharmony_ci @debug("leastsquares", "vector=", repr(vector)) 398bbbf1280Sopenharmony_ci 399bbbf1280Sopenharmony_ci # Solve the matrix equation. 400bbbf1280Sopenharmony_ci all_coeffs = matrix \ vector 401bbbf1280Sopenharmony_ci 402bbbf1280Sopenharmony_ci @debug("leastsquares", "all_coeffs=", repr(all_coeffs)) 403bbbf1280Sopenharmony_ci 404bbbf1280Sopenharmony_ci # And marshal the results into two separate polynomial vectors to 405bbbf1280Sopenharmony_ci # return. 406bbbf1280Sopenharmony_ci ncoeffs = all_coeffs[1:n+1] 407bbbf1280Sopenharmony_ci dcoeffs = vcat([1], all_coeffs[n+2:n+d+1]) 408bbbf1280Sopenharmony_ci return (ncoeffs, dcoeffs) 409bbbf1280Sopenharmony_ciend 410bbbf1280Sopenharmony_ci 411bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 412bbbf1280Sopenharmony_ci# Golden-section search to find a maximum of a function. 413bbbf1280Sopenharmony_ci 414bbbf1280Sopenharmony_ci# Arguments: 415bbbf1280Sopenharmony_ci# f Function to be maximised/minimised. Maps BigFloat -> BigFloat. 416bbbf1280Sopenharmony_ci# a,b,c BigFloats bracketing a maximum of the function. 417bbbf1280Sopenharmony_ci# 418bbbf1280Sopenharmony_ci# Expects: 419bbbf1280Sopenharmony_ci# a,b,c are in order (either a<=b<=c or c<=b<=a) 420bbbf1280Sopenharmony_ci# a != c (but b can equal one or the other if it wants to) 421bbbf1280Sopenharmony_ci# f(a) <= f(b) >= f(c) 422bbbf1280Sopenharmony_ci# 423bbbf1280Sopenharmony_ci# Return value is an (x,y) pair of BigFloats giving the extremal input 424bbbf1280Sopenharmony_ci# and output. (That is, y=f(x).) 425bbbf1280Sopenharmony_cifunction goldensection(f::Function, a::BigFloat, b::BigFloat, c::BigFloat) 426bbbf1280Sopenharmony_ci # Decide on a 'good enough' threshold. 427bbbf1280Sopenharmony_ci threshold = abs(c-a) * 2^(-epsbits/2) 428bbbf1280Sopenharmony_ci 429bbbf1280Sopenharmony_ci # We'll need the golden ratio phi, of course. Or rather, in this 430bbbf1280Sopenharmony_ci # case, we need 1/phi = 0.618... 431bbbf1280Sopenharmony_ci one_over_phi = 2 / (1 + sqrt(BigFloat(5))) 432bbbf1280Sopenharmony_ci 433bbbf1280Sopenharmony_ci # Flip round the interval endpoints so that the interval [a,b] is 434bbbf1280Sopenharmony_ci # at least as large as [b,c]. (Then we can always pick our new 435bbbf1280Sopenharmony_ci # point in [a,b] without having to handle lots of special cases.) 436bbbf1280Sopenharmony_ci if abs(b-a) < abs(c-a) 437bbbf1280Sopenharmony_ci a, c = c, a 438bbbf1280Sopenharmony_ci end 439bbbf1280Sopenharmony_ci 440bbbf1280Sopenharmony_ci # Evaluate the function at the initial points. 441bbbf1280Sopenharmony_ci fa = f(a) 442bbbf1280Sopenharmony_ci fb = f(b) 443bbbf1280Sopenharmony_ci fc = f(c) 444bbbf1280Sopenharmony_ci 445bbbf1280Sopenharmony_ci @debug("goldensection", "starting") 446bbbf1280Sopenharmony_ci 447bbbf1280Sopenharmony_ci while abs(c-a) > threshold 448bbbf1280Sopenharmony_ci @debug("goldensection", "a: ", a, " -> ", fa) 449bbbf1280Sopenharmony_ci @debug("goldensection", "b: ", b, " -> ", fb) 450bbbf1280Sopenharmony_ci @debug("goldensection", "c: ", c, " -> ", fc) 451bbbf1280Sopenharmony_ci 452bbbf1280Sopenharmony_ci # Check invariants. 453bbbf1280Sopenharmony_ci @assert(a <= b <= c || c <= b <= a) 454bbbf1280Sopenharmony_ci @assert(fa <= fb >= fc) 455bbbf1280Sopenharmony_ci 456bbbf1280Sopenharmony_ci # Subdivide the larger of the intervals [a,b] and [b,c]. We've 457bbbf1280Sopenharmony_ci # arranged that this is always [a,b], for simplicity. 458bbbf1280Sopenharmony_ci d = a + (b-a) * one_over_phi 459bbbf1280Sopenharmony_ci 460bbbf1280Sopenharmony_ci # Now we have an interval looking like this (possibly 461bbbf1280Sopenharmony_ci # reversed): 462bbbf1280Sopenharmony_ci # 463bbbf1280Sopenharmony_ci # a d b c 464bbbf1280Sopenharmony_ci # 465bbbf1280Sopenharmony_ci # and we know f(b) is bigger than either f(a) or f(c). We have 466bbbf1280Sopenharmony_ci # two cases: either f(d) > f(b), or vice versa. In either 467bbbf1280Sopenharmony_ci # case, we can narrow to an interval of 1/phi the size, and 468bbbf1280Sopenharmony_ci # still satisfy all our invariants (three ordered points, 469bbbf1280Sopenharmony_ci # [a,b] at least the width of [b,c], f(a)<=f(b)>=f(c)). 470bbbf1280Sopenharmony_ci fd = f(d) 471bbbf1280Sopenharmony_ci @debug("goldensection", "d: ", d, " -> ", fd) 472bbbf1280Sopenharmony_ci if fd > fb 473bbbf1280Sopenharmony_ci a, b, c = a, d, b 474bbbf1280Sopenharmony_ci fa, fb, fc = fa, fd, fb 475bbbf1280Sopenharmony_ci @debug("goldensection", "adb case") 476bbbf1280Sopenharmony_ci else 477bbbf1280Sopenharmony_ci a, b, c = c, b, d 478bbbf1280Sopenharmony_ci fa, fb, fc = fc, fb, fd 479bbbf1280Sopenharmony_ci @debug("goldensection", "cbd case") 480bbbf1280Sopenharmony_ci end 481bbbf1280Sopenharmony_ci end 482bbbf1280Sopenharmony_ci 483bbbf1280Sopenharmony_ci @debug("goldensection", "done: ", b, " -> ", fb) 484bbbf1280Sopenharmony_ci return (b, fb) 485bbbf1280Sopenharmony_ciend 486bbbf1280Sopenharmony_ci 487bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 488bbbf1280Sopenharmony_ci# Find the extrema of a function within a given interval. 489bbbf1280Sopenharmony_ci 490bbbf1280Sopenharmony_ci# Arguments: 491bbbf1280Sopenharmony_ci# f The function to be approximated. Maps BigFloat -> BigFloat. 492bbbf1280Sopenharmony_ci# grid A set of points at which to evaluate f. Must be high enough 493bbbf1280Sopenharmony_ci# resolution to make extrema obvious. 494bbbf1280Sopenharmony_ci# 495bbbf1280Sopenharmony_ci# Returns an array of (x,y) pairs of BigFloats, with each x,y giving 496bbbf1280Sopenharmony_ci# the extremum location and its value (i.e. y=f(x)). 497bbbf1280Sopenharmony_cifunction find_extrema(f::Function, grid::Array{BigFloat}) 498bbbf1280Sopenharmony_ci len = length(grid) 499bbbf1280Sopenharmony_ci extrema = array1d(Tuple{BigFloat, BigFloat}, 0) 500bbbf1280Sopenharmony_ci for i = 1:1:len 501bbbf1280Sopenharmony_ci # We have to provide goldensection() with three points 502bbbf1280Sopenharmony_ci # bracketing the extremum. If the extremum is at one end of 503bbbf1280Sopenharmony_ci # the interval, then the only way we can do that is to set two 504bbbf1280Sopenharmony_ci # of the points equal (which goldensection() will cope with). 505bbbf1280Sopenharmony_ci prev = max(1, i-1) 506bbbf1280Sopenharmony_ci next = min(i+1, len) 507bbbf1280Sopenharmony_ci 508bbbf1280Sopenharmony_ci # Find our three pairs of (x,y) coordinates. 509bbbf1280Sopenharmony_ci xp, xi, xn = grid[prev], grid[i], grid[next] 510bbbf1280Sopenharmony_ci yp, yi, yn = f(xp), f(xi), f(xn) 511bbbf1280Sopenharmony_ci 512bbbf1280Sopenharmony_ci # See if they look like an extremum, and if so, ask 513bbbf1280Sopenharmony_ci # goldensection() to give a more exact location for it. 514bbbf1280Sopenharmony_ci if yp <= yi >= yn 515bbbf1280Sopenharmony_ci push!(extrema, goldensection(f, xp, xi, xn)) 516bbbf1280Sopenharmony_ci elseif yp >= yi <= yn 517bbbf1280Sopenharmony_ci x, y = goldensection(x->-f(x), xp, xi, xn) 518bbbf1280Sopenharmony_ci push!(extrema, (x, -y)) 519bbbf1280Sopenharmony_ci end 520bbbf1280Sopenharmony_ci end 521bbbf1280Sopenharmony_ci return extrema 522bbbf1280Sopenharmony_ciend 523bbbf1280Sopenharmony_ci 524bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 525bbbf1280Sopenharmony_ci# Winnow a list of a function's extrema to give a subsequence of a 526bbbf1280Sopenharmony_ci# specified length, with the extrema in the subsequence alternating 527bbbf1280Sopenharmony_ci# signs, and with the smallest absolute value of an extremum in the 528bbbf1280Sopenharmony_ci# subsequence as large as possible. 529bbbf1280Sopenharmony_ci# 530bbbf1280Sopenharmony_ci# We do this using a dynamic-programming approach. We work along the 531bbbf1280Sopenharmony_ci# provided array of extrema, and at all times, we track the best set 532bbbf1280Sopenharmony_ci# of extrema we have so far seen for each possible (length, sign of 533bbbf1280Sopenharmony_ci# last extremum) pair. Each new extremum is evaluated to see whether 534bbbf1280Sopenharmony_ci# it can be added to any previously seen best subsequence to make a 535bbbf1280Sopenharmony_ci# new subsequence that beats the previous record holder in its slot. 536bbbf1280Sopenharmony_ci 537bbbf1280Sopenharmony_ci# Arguments: 538bbbf1280Sopenharmony_ci# extrema An array of (x,y) pairs of BigFloats giving the input extrema. 539bbbf1280Sopenharmony_ci# n Number of extrema required as output. 540bbbf1280Sopenharmony_ci# 541bbbf1280Sopenharmony_ci# Returns a new array of (x,y) pairs which is a subsequence of the 542bbbf1280Sopenharmony_ci# original sequence. (So, in particular, if the input was sorted by x 543bbbf1280Sopenharmony_ci# then so will the output be.) 544bbbf1280Sopenharmony_cifunction winnow_extrema(extrema::Array{Tuple{BigFloat,BigFloat}}, n) 545bbbf1280Sopenharmony_ci # best[i,j] gives the best sequence so far of length i and with 546bbbf1280Sopenharmony_ci # sign j (where signs are coded as 1=positive, 2=negative), in the 547bbbf1280Sopenharmony_ci # form of a tuple (cost, actual array of x,y pairs). 548bbbf1280Sopenharmony_ci best = fill((BigFloat(0), array1d(Tuple{BigFloat,BigFloat}, 0)), n, 2) 549bbbf1280Sopenharmony_ci 550bbbf1280Sopenharmony_ci for (x,y) = extrema 551bbbf1280Sopenharmony_ci if y > 0 552bbbf1280Sopenharmony_ci sign = 1 553bbbf1280Sopenharmony_ci elseif y < 0 554bbbf1280Sopenharmony_ci sign = 2 555bbbf1280Sopenharmony_ci else 556bbbf1280Sopenharmony_ci # A zero-valued extremum cannot possibly contribute to any 557bbbf1280Sopenharmony_ci # optimal sequence, so we simply ignore it! 558bbbf1280Sopenharmony_ci continue 559bbbf1280Sopenharmony_ci end 560bbbf1280Sopenharmony_ci 561bbbf1280Sopenharmony_ci for i = 1:1:n 562bbbf1280Sopenharmony_ci # See if we can create a new entry for best[i,sign] by 563bbbf1280Sopenharmony_ci # appending our current (x,y) to some previous thing. 564bbbf1280Sopenharmony_ci if i == 1 565bbbf1280Sopenharmony_ci # Special case: we don't store a best zero-length 566bbbf1280Sopenharmony_ci # sequence :-) 567bbbf1280Sopenharmony_ci candidate = (abs(y), [(x,y)]) 568bbbf1280Sopenharmony_ci else 569bbbf1280Sopenharmony_ci othersign = 3-sign # map 1->2 and 2->1 570bbbf1280Sopenharmony_ci oldscore, oldlist = best[i-1, othersign] 571bbbf1280Sopenharmony_ci newscore = min(abs(y), oldscore) 572bbbf1280Sopenharmony_ci newlist = vcat(oldlist, [(x,y)]) 573bbbf1280Sopenharmony_ci candidate = (newscore, newlist) 574bbbf1280Sopenharmony_ci end 575bbbf1280Sopenharmony_ci # If our new candidate improves on the previous value of 576bbbf1280Sopenharmony_ci # best[i,sign], then replace it. 577bbbf1280Sopenharmony_ci if candidate[1] > best[i,sign][1] 578bbbf1280Sopenharmony_ci best[i,sign] = candidate 579bbbf1280Sopenharmony_ci end 580bbbf1280Sopenharmony_ci end 581bbbf1280Sopenharmony_ci end 582bbbf1280Sopenharmony_ci 583bbbf1280Sopenharmony_ci # Our ultimate return value has to be either best[n,1] or 584bbbf1280Sopenharmony_ci # best[n,2], but it could be either. See which one has the higher 585bbbf1280Sopenharmony_ci # score. 586bbbf1280Sopenharmony_ci if best[n,1][1] > best[n,2][1] 587bbbf1280Sopenharmony_ci ret = best[n,1][2] 588bbbf1280Sopenharmony_ci else 589bbbf1280Sopenharmony_ci ret = best[n,2][2] 590bbbf1280Sopenharmony_ci end 591bbbf1280Sopenharmony_ci # Make sure we did actually _find_ a good answer. 592bbbf1280Sopenharmony_ci @assert(length(ret) == n) 593bbbf1280Sopenharmony_ci return ret 594bbbf1280Sopenharmony_ciend 595bbbf1280Sopenharmony_ci 596bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 597bbbf1280Sopenharmony_ci# Construct a rational-function approximation with equal and 598bbbf1280Sopenharmony_ci# alternating weighted deviation at a specific set of x-coordinates. 599bbbf1280Sopenharmony_ci 600bbbf1280Sopenharmony_ci# Arguments: 601bbbf1280Sopenharmony_ci# f The function to be approximated. Maps BigFloat -> BigFloat. 602bbbf1280Sopenharmony_ci# coords An array of BigFloats giving the x-coordinates. There should 603bbbf1280Sopenharmony_ci# be n+d+2 of them. 604bbbf1280Sopenharmony_ci# n, d The degrees of the numerator and denominator of the desired 605bbbf1280Sopenharmony_ci# approximation. 606bbbf1280Sopenharmony_ci# prev_err A plausible value for the alternating weighted deviation. 607bbbf1280Sopenharmony_ci# (Required to kickstart a binary search in the nonlinear case; 608bbbf1280Sopenharmony_ci# see comments below.) 609bbbf1280Sopenharmony_ci# w Error-weighting function. Takes two BigFloat arguments x,y 610bbbf1280Sopenharmony_ci# and returns a scaling factor for the error at that location. 611bbbf1280Sopenharmony_ci# The returned approximation R should have the minimum possible 612bbbf1280Sopenharmony_ci# maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional 613bbbf1280Sopenharmony_ci# parameter, defaulting to the always-return-1 function. 614bbbf1280Sopenharmony_ci# 615bbbf1280Sopenharmony_ci# Return values: a pair of arrays of BigFloats (N,D) giving the 616bbbf1280Sopenharmony_ci# coefficients of the returned rational function. N has size n+1; D 617bbbf1280Sopenharmony_ci# has size d+1. Both start with the constant term, i.e. N[i] is the 618bbbf1280Sopenharmony_ci# coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will 619bbbf1280Sopenharmony_ci# be 1. 620bbbf1280Sopenharmony_cifunction ratfn_equal_deviation(f::Function, coords::Array{BigFloat}, 621bbbf1280Sopenharmony_ci n, d, prev_err::BigFloat, 622bbbf1280Sopenharmony_ci w = (x,y)->BigFloat(1)) 623bbbf1280Sopenharmony_ci @debug("equaldev", "n=", n, " d=", d, " coords=", repr(coords)) 624bbbf1280Sopenharmony_ci @assert(length(coords) == n+d+2) 625bbbf1280Sopenharmony_ci 626bbbf1280Sopenharmony_ci if d == 0 627bbbf1280Sopenharmony_ci # Special case: we're after a polynomial. In this case, we 628bbbf1280Sopenharmony_ci # have the particularly easy job of just constructing and 629bbbf1280Sopenharmony_ci # solving a system of n+2 linear equations, to find the n+1 630bbbf1280Sopenharmony_ci # coefficients of the polynomial and also the amount of 631bbbf1280Sopenharmony_ci # deviation at the specified coordinates. Each equation is of 632bbbf1280Sopenharmony_ci # the form 633bbbf1280Sopenharmony_ci # 634bbbf1280Sopenharmony_ci # p_0 x^0 + p_1 x^1 + ... + p_n x^n ± e/w(x) = f(x) 635bbbf1280Sopenharmony_ci # 636bbbf1280Sopenharmony_ci # in which the p_i and e are the variables, and the powers of 637bbbf1280Sopenharmony_ci # x and calls to w and f are the coefficients. 638bbbf1280Sopenharmony_ci 639bbbf1280Sopenharmony_ci matrix = array2d(BigFloat, n+2, n+2) 640bbbf1280Sopenharmony_ci vector = array1d(BigFloat, n+2) 641bbbf1280Sopenharmony_ci currsign = +1 642bbbf1280Sopenharmony_ci for i = 1:1:n+2 643bbbf1280Sopenharmony_ci x = coords[i] 644bbbf1280Sopenharmony_ci for j = 0:1:n 645bbbf1280Sopenharmony_ci matrix[i,1+j] = x^j 646bbbf1280Sopenharmony_ci end 647bbbf1280Sopenharmony_ci y = f(x) 648bbbf1280Sopenharmony_ci vector[i] = y 649bbbf1280Sopenharmony_ci matrix[i, n+2] = currsign / w(x,y) 650bbbf1280Sopenharmony_ci currsign = -currsign 651bbbf1280Sopenharmony_ci end 652bbbf1280Sopenharmony_ci 653bbbf1280Sopenharmony_ci @debug("equaldev", "matrix=", repr(matrix)) 654bbbf1280Sopenharmony_ci @debug("equaldev", "vector=", repr(vector)) 655bbbf1280Sopenharmony_ci 656bbbf1280Sopenharmony_ci outvector = matrix \ vector 657bbbf1280Sopenharmony_ci 658bbbf1280Sopenharmony_ci @debug("equaldev", "outvector=", repr(outvector)) 659bbbf1280Sopenharmony_ci 660bbbf1280Sopenharmony_ci ncoeffs = outvector[1:n+1] 661bbbf1280Sopenharmony_ci dcoeffs = [BigFloat(1)] 662bbbf1280Sopenharmony_ci return ncoeffs, dcoeffs 663bbbf1280Sopenharmony_ci else 664bbbf1280Sopenharmony_ci # For a nontrivial rational function, the system of equations 665bbbf1280Sopenharmony_ci # we need to solve becomes nonlinear, because each equation 666bbbf1280Sopenharmony_ci # now takes the form 667bbbf1280Sopenharmony_ci # 668bbbf1280Sopenharmony_ci # p_0 x^0 + p_1 x^1 + ... + p_n x^n 669bbbf1280Sopenharmony_ci # --------------------------------- ± e/w(x) = f(x) 670bbbf1280Sopenharmony_ci # x^0 + q_1 x^1 + ... + q_d x^d 671bbbf1280Sopenharmony_ci # 672bbbf1280Sopenharmony_ci # and multiplying up by the denominator gives you a lot of 673bbbf1280Sopenharmony_ci # terms containing e × q_i. So we can't do this the really 674bbbf1280Sopenharmony_ci # easy way using a matrix equation as above. 675bbbf1280Sopenharmony_ci # 676bbbf1280Sopenharmony_ci # Fortunately, this is a fairly easy kind of nonlinear system. 677bbbf1280Sopenharmony_ci # The equations all become linear if you switch to treating e 678bbbf1280Sopenharmony_ci # as a constant, so a reasonably sensible approach is to pick 679bbbf1280Sopenharmony_ci # a candidate value of e, solve all but one of the equations 680bbbf1280Sopenharmony_ci # for the remaining unknowns, and then see what the error 681bbbf1280Sopenharmony_ci # turns out to be in the final equation. The Chebyshev 682bbbf1280Sopenharmony_ci # alternation theorem guarantees that that error in the last 683bbbf1280Sopenharmony_ci # equation will be anti-monotonic in the input e, so we can 684bbbf1280Sopenharmony_ci # just binary-search until we get the two as close to equal as 685bbbf1280Sopenharmony_ci # we need them. 686bbbf1280Sopenharmony_ci 687bbbf1280Sopenharmony_ci function try_e(e) 688bbbf1280Sopenharmony_ci # Try a given value of e, derive the coefficients of the 689bbbf1280Sopenharmony_ci # resulting rational function by setting up equations 690bbbf1280Sopenharmony_ci # based on the first n+d+1 of the n+d+2 coordinates, and 691bbbf1280Sopenharmony_ci # see what the error turns out to be at the final 692bbbf1280Sopenharmony_ci # coordinate. 693bbbf1280Sopenharmony_ci matrix = array2d(BigFloat, n+d+1, n+d+1) 694bbbf1280Sopenharmony_ci vector = array1d(BigFloat, n+d+1) 695bbbf1280Sopenharmony_ci currsign = +1 696bbbf1280Sopenharmony_ci for i = 1:1:n+d+1 697bbbf1280Sopenharmony_ci x = coords[i] 698bbbf1280Sopenharmony_ci y = f(x) 699bbbf1280Sopenharmony_ci y_adj = y - currsign * e / w(x,y) 700bbbf1280Sopenharmony_ci for j = 0:1:n 701bbbf1280Sopenharmony_ci matrix[i,1+j] = x^j 702bbbf1280Sopenharmony_ci end 703bbbf1280Sopenharmony_ci for j = 1:1:d 704bbbf1280Sopenharmony_ci matrix[i,1+n+j] = -x^j * y_adj 705bbbf1280Sopenharmony_ci end 706bbbf1280Sopenharmony_ci vector[i] = y_adj 707bbbf1280Sopenharmony_ci currsign = -currsign 708bbbf1280Sopenharmony_ci end 709bbbf1280Sopenharmony_ci 710bbbf1280Sopenharmony_ci @debug("equaldev", "trying e=", e) 711bbbf1280Sopenharmony_ci @debug("equaldev", "matrix=", repr(matrix)) 712bbbf1280Sopenharmony_ci @debug("equaldev", "vector=", repr(vector)) 713bbbf1280Sopenharmony_ci 714bbbf1280Sopenharmony_ci outvector = matrix \ vector 715bbbf1280Sopenharmony_ci 716bbbf1280Sopenharmony_ci @debug("equaldev", "outvector=", repr(outvector)) 717bbbf1280Sopenharmony_ci 718bbbf1280Sopenharmony_ci ncoeffs = outvector[1:n+1] 719bbbf1280Sopenharmony_ci dcoeffs = vcat([BigFloat(1)], outvector[n+2:n+d+1]) 720bbbf1280Sopenharmony_ci 721bbbf1280Sopenharmony_ci x = coords[n+d+2] 722bbbf1280Sopenharmony_ci y = f(x) 723bbbf1280Sopenharmony_ci last_e = (ratfn_eval(ncoeffs, dcoeffs, x) - y) * w(x,y) * -currsign 724bbbf1280Sopenharmony_ci 725bbbf1280Sopenharmony_ci @debug("equaldev", "last e=", last_e) 726bbbf1280Sopenharmony_ci 727bbbf1280Sopenharmony_ci return ncoeffs, dcoeffs, last_e 728bbbf1280Sopenharmony_ci end 729bbbf1280Sopenharmony_ci 730bbbf1280Sopenharmony_ci threshold = 2^(-epsbits/2) # convergence threshold 731bbbf1280Sopenharmony_ci 732bbbf1280Sopenharmony_ci # Start by trying our previous iteration's error value. This 733bbbf1280Sopenharmony_ci # value (e0) will be one end of our binary-search interval, 734bbbf1280Sopenharmony_ci # and whatever it caused the last point's error to be, that 735bbbf1280Sopenharmony_ci # (e1) will be the other end. 736bbbf1280Sopenharmony_ci e0 = prev_err 737bbbf1280Sopenharmony_ci @debug("equaldev", "e0 = ", e0) 738bbbf1280Sopenharmony_ci nc, dc, e1 = try_e(e0) 739bbbf1280Sopenharmony_ci @debug("equaldev", "e1 = ", e1) 740bbbf1280Sopenharmony_ci if abs(e1-e0) <= threshold 741bbbf1280Sopenharmony_ci # If we're _really_ lucky, we hit the error right on the 742bbbf1280Sopenharmony_ci # nose just by doing that! 743bbbf1280Sopenharmony_ci return nc, dc 744bbbf1280Sopenharmony_ci end 745bbbf1280Sopenharmony_ci s = sign(e1-e0) 746bbbf1280Sopenharmony_ci @debug("equaldev", "s = ", s) 747bbbf1280Sopenharmony_ci 748bbbf1280Sopenharmony_ci # Verify by assertion that trying our other interval endpoint 749bbbf1280Sopenharmony_ci # e1 gives a value that's wrong in the other direction. 750bbbf1280Sopenharmony_ci # (Otherwise our binary search won't get a sensible answer at 751bbbf1280Sopenharmony_ci # all.) 752bbbf1280Sopenharmony_ci nc, dc, e2 = try_e(e1) 753bbbf1280Sopenharmony_ci @debug("equaldev", "e2 = ", e2) 754bbbf1280Sopenharmony_ci @assert(sign(e2-e1) == -s) 755bbbf1280Sopenharmony_ci 756bbbf1280Sopenharmony_ci # Now binary-search until our two endpoints narrow enough. 757bbbf1280Sopenharmony_ci local emid 758bbbf1280Sopenharmony_ci while abs(e1-e0) > threshold 759bbbf1280Sopenharmony_ci emid = (e1+e0)/2 760bbbf1280Sopenharmony_ci nc, dc, enew = try_e(emid) 761bbbf1280Sopenharmony_ci if sign(enew-emid) == s 762bbbf1280Sopenharmony_ci e0 = emid 763bbbf1280Sopenharmony_ci else 764bbbf1280Sopenharmony_ci e1 = emid 765bbbf1280Sopenharmony_ci end 766bbbf1280Sopenharmony_ci end 767bbbf1280Sopenharmony_ci 768bbbf1280Sopenharmony_ci @debug("equaldev", "final e=", emid) 769bbbf1280Sopenharmony_ci return nc, dc 770bbbf1280Sopenharmony_ci end 771bbbf1280Sopenharmony_ciend 772bbbf1280Sopenharmony_ci 773bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 774bbbf1280Sopenharmony_ci# Top-level function to find a minimax rational-function approximation. 775bbbf1280Sopenharmony_ci 776bbbf1280Sopenharmony_ci# Arguments: 777bbbf1280Sopenharmony_ci# f The function to be approximated. Maps BigFloat -> BigFloat. 778bbbf1280Sopenharmony_ci# interval A pair of BigFloats giving the endpoints of the interval 779bbbf1280Sopenharmony_ci# (in either order) on which to approximate f. 780bbbf1280Sopenharmony_ci# n, d The degrees of the numerator and denominator of the desired 781bbbf1280Sopenharmony_ci# approximation. 782bbbf1280Sopenharmony_ci# w Error-weighting function. Takes two BigFloat arguments x,y 783bbbf1280Sopenharmony_ci# and returns a scaling factor for the error at that location. 784bbbf1280Sopenharmony_ci# The returned approximation R should have the minimum possible 785bbbf1280Sopenharmony_ci# maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional 786bbbf1280Sopenharmony_ci# parameter, defaulting to the always-return-1 function. 787bbbf1280Sopenharmony_ci# 788bbbf1280Sopenharmony_ci# Return values: a tuple (N,D,E,X), where 789bbbf1280Sopenharmony_ci 790bbbf1280Sopenharmony_ci# N,D A pair of arrays of BigFloats giving the coefficients 791bbbf1280Sopenharmony_ci# of the returned rational function. N has size n+1; D 792bbbf1280Sopenharmony_ci# has size d+1. Both start with the constant term, i.e. 793bbbf1280Sopenharmony_ci# N[i] is the coefficient of x^(i-1) (because Julia 794bbbf1280Sopenharmony_ci# arrays are 1-based). D[1] will be 1. 795bbbf1280Sopenharmony_ci# E The maximum weighted error (BigFloat). 796bbbf1280Sopenharmony_ci# X An array of pairs of BigFloats giving the locations of n+2 797bbbf1280Sopenharmony_ci# points and the weighted error at each of those points. The 798bbbf1280Sopenharmony_ci# weighted error values will have alternating signs, which 799bbbf1280Sopenharmony_ci# means that the Chebyshev alternation theorem guarantees 800bbbf1280Sopenharmony_ci# that any other function of the same degree must exceed 801bbbf1280Sopenharmony_ci# the error of this one at at least one of those points. 802bbbf1280Sopenharmony_cifunction ratfn_minimax(f::Function, interval::Tuple{BigFloat,BigFloat}, n, d, 803bbbf1280Sopenharmony_ci w = (x,y)->BigFloat(1)) 804bbbf1280Sopenharmony_ci # We start off by finding a least-squares approximation. This 805bbbf1280Sopenharmony_ci # doesn't need to be perfect, but if we can get it reasonably good 806bbbf1280Sopenharmony_ci # then it'll save iterations in the refining stage. 807bbbf1280Sopenharmony_ci # 808bbbf1280Sopenharmony_ci # Least-squares approximations tend to look nicer in a minimax 809bbbf1280Sopenharmony_ci # sense if you evaluate the function at a big pile of Chebyshev 810bbbf1280Sopenharmony_ci # nodes rather than uniformly spaced points. These values will 811bbbf1280Sopenharmony_ci # also make a good grid to use for the initial search for error 812bbbf1280Sopenharmony_ci # extrema, so we'll keep them around for that reason too. 813bbbf1280Sopenharmony_ci 814bbbf1280Sopenharmony_ci # Construct the grid. 815bbbf1280Sopenharmony_ci lo, hi = minimum(interval), maximum(interval) 816bbbf1280Sopenharmony_ci local grid 817bbbf1280Sopenharmony_ci let 818bbbf1280Sopenharmony_ci mid = (hi+lo)/2 819bbbf1280Sopenharmony_ci halfwid = (hi-lo)/2 820bbbf1280Sopenharmony_ci nnodes = 16 * (n+d+1) 821bbbf1280Sopenharmony_ci pi = 2*asin(BigFloat(1)) 822bbbf1280Sopenharmony_ci grid = [ mid - halfwid * cos(pi*i/nnodes) for i=0:1:nnodes ] 823bbbf1280Sopenharmony_ci end 824bbbf1280Sopenharmony_ci 825bbbf1280Sopenharmony_ci # Find the initial least-squares approximation. 826bbbf1280Sopenharmony_ci (nc, dc) = ratfn_leastsquares(f, grid, n, d, w) 827bbbf1280Sopenharmony_ci @debug("minimax", "initial leastsquares approx = ", 828bbbf1280Sopenharmony_ci ratfn_to_string(nc, dc)) 829bbbf1280Sopenharmony_ci 830bbbf1280Sopenharmony_ci # Threshold of convergence. We stop when the relative difference 831bbbf1280Sopenharmony_ci # between the min and max (winnowed) error extrema is less than 832bbbf1280Sopenharmony_ci # this. 833bbbf1280Sopenharmony_ci # 834bbbf1280Sopenharmony_ci # This is set to the cube root of machine epsilon on a more or 835bbbf1280Sopenharmony_ci # less empirical basis, because the rational-function case will 836bbbf1280Sopenharmony_ci # not converge reliably if you set it to only the square root. 837bbbf1280Sopenharmony_ci # (Repeatable by using the --test mode.) On the assumption that 838bbbf1280Sopenharmony_ci # input and output error in each iteration can be expected to be 839bbbf1280Sopenharmony_ci # related by a simple power law (because it'll just be down to how 840bbbf1280Sopenharmony_ci # many leading terms of a Taylor series are zero), the cube root 841bbbf1280Sopenharmony_ci # was the next thing to try. 842bbbf1280Sopenharmony_ci threshold = 2^(-epsbits/3) 843bbbf1280Sopenharmony_ci 844bbbf1280Sopenharmony_ci # Main loop. 845bbbf1280Sopenharmony_ci while true 846bbbf1280Sopenharmony_ci # Find all the error extrema we can. 847bbbf1280Sopenharmony_ci function compute_error(x) 848bbbf1280Sopenharmony_ci real_y = f(x) 849bbbf1280Sopenharmony_ci approx_y = ratfn_eval(nc, dc, x) 850bbbf1280Sopenharmony_ci return (approx_y - real_y) * w(x, real_y) 851bbbf1280Sopenharmony_ci end 852bbbf1280Sopenharmony_ci extrema = find_extrema(compute_error, grid) 853bbbf1280Sopenharmony_ci @debug("minimax", "all extrema = ", format_xylist(extrema)) 854bbbf1280Sopenharmony_ci 855bbbf1280Sopenharmony_ci # Winnow the extrema down to the right number, and ensure they 856bbbf1280Sopenharmony_ci # have alternating sign. 857bbbf1280Sopenharmony_ci extrema = winnow_extrema(extrema, n+d+2) 858bbbf1280Sopenharmony_ci @debug("minimax", "winnowed extrema = ", format_xylist(extrema)) 859bbbf1280Sopenharmony_ci 860bbbf1280Sopenharmony_ci # See if we've finished. 861bbbf1280Sopenharmony_ci min_err = minimum([abs(y) for (x,y) = extrema]) 862bbbf1280Sopenharmony_ci max_err = maximum([abs(y) for (x,y) = extrema]) 863bbbf1280Sopenharmony_ci variation = (max_err - min_err) / max_err 864bbbf1280Sopenharmony_ci @debug("minimax", "extremum variation = ", variation) 865bbbf1280Sopenharmony_ci if variation < threshold 866bbbf1280Sopenharmony_ci @debug("minimax", "done!") 867bbbf1280Sopenharmony_ci return nc, dc, max_err, extrema 868bbbf1280Sopenharmony_ci end 869bbbf1280Sopenharmony_ci 870bbbf1280Sopenharmony_ci # If not, refine our function by equalising the error at the 871bbbf1280Sopenharmony_ci # extrema points, and go round again. 872bbbf1280Sopenharmony_ci (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema), 873bbbf1280Sopenharmony_ci n, d, max_err, w) 874bbbf1280Sopenharmony_ci @debug("minimax", "refined approx = ", ratfn_to_string(nc, dc)) 875bbbf1280Sopenharmony_ci end 876bbbf1280Sopenharmony_ciend 877bbbf1280Sopenharmony_ci 878bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 879bbbf1280Sopenharmony_ci# Check if a polynomial is well-conditioned for accurate evaluation in 880bbbf1280Sopenharmony_ci# a given interval by Horner's rule. 881bbbf1280Sopenharmony_ci# 882bbbf1280Sopenharmony_ci# This is true if at every step where Horner's rule computes 883bbbf1280Sopenharmony_ci# (coefficient + x*value_so_far), the constant coefficient you're 884bbbf1280Sopenharmony_ci# adding on is of larger magnitude than the x*value_so_far operand. 885bbbf1280Sopenharmony_ci# And this has to be true for every x in the interval. 886bbbf1280Sopenharmony_ci# 887bbbf1280Sopenharmony_ci# Arguments: 888bbbf1280Sopenharmony_ci# coeffs The coefficients of the polynomial under test. Starts with 889bbbf1280Sopenharmony_ci# the constant term, i.e. coeffs[i] is the coefficient of 890bbbf1280Sopenharmony_ci# x^(i-1) (because Julia arrays are 1-based). 891bbbf1280Sopenharmony_ci# lo, hi The bounds of the interval. 892bbbf1280Sopenharmony_ci# 893bbbf1280Sopenharmony_ci# Return value: the largest ratio (x*value_so_far / coefficient), at 894bbbf1280Sopenharmony_ci# any step of evaluation, for any x in the interval. If this is less 895bbbf1280Sopenharmony_ci# than 1, the polynomial is at least somewhat well-conditioned; 896bbbf1280Sopenharmony_ci# ideally you want it to be more like 1/8 or 1/16 or so, so that the 897bbbf1280Sopenharmony_ci# relative rounding error accumulated at each step are reduced by 898bbbf1280Sopenharmony_ci# several factors of 2 when the next coefficient is added on. 899bbbf1280Sopenharmony_ci 900bbbf1280Sopenharmony_cifunction wellcond(coeffs, lo, hi) 901bbbf1280Sopenharmony_ci x = max(abs(lo), abs(hi)) 902bbbf1280Sopenharmony_ci worst = 0 903bbbf1280Sopenharmony_ci so_far = 0 904bbbf1280Sopenharmony_ci for i = length(coeffs):-1:1 905bbbf1280Sopenharmony_ci coeff = abs(coeffs[i]) 906bbbf1280Sopenharmony_ci so_far *= x 907bbbf1280Sopenharmony_ci if coeff != 0 908bbbf1280Sopenharmony_ci thisval = so_far / coeff 909bbbf1280Sopenharmony_ci worst = max(worst, thisval) 910bbbf1280Sopenharmony_ci so_far += coeff 911bbbf1280Sopenharmony_ci end 912bbbf1280Sopenharmony_ci end 913bbbf1280Sopenharmony_ci return worst 914bbbf1280Sopenharmony_ciend 915bbbf1280Sopenharmony_ci 916bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 917bbbf1280Sopenharmony_ci# Small set of unit tests. 918bbbf1280Sopenharmony_ci 919bbbf1280Sopenharmony_cifunction test() 920bbbf1280Sopenharmony_ci passes = 0 921bbbf1280Sopenharmony_ci fails = 0 922bbbf1280Sopenharmony_ci 923bbbf1280Sopenharmony_ci function approx_eq(x, y, limit=1e-6) 924bbbf1280Sopenharmony_ci return abs(x - y) < limit 925bbbf1280Sopenharmony_ci end 926bbbf1280Sopenharmony_ci 927bbbf1280Sopenharmony_ci function test(condition) 928bbbf1280Sopenharmony_ci if condition 929bbbf1280Sopenharmony_ci passes += 1 930bbbf1280Sopenharmony_ci else 931bbbf1280Sopenharmony_ci println("fail") 932bbbf1280Sopenharmony_ci fails += 1 933bbbf1280Sopenharmony_ci end 934bbbf1280Sopenharmony_ci end 935bbbf1280Sopenharmony_ci 936bbbf1280Sopenharmony_ci # Test Gaussian elimination. 937bbbf1280Sopenharmony_ci println("Gaussian test 1:") 938bbbf1280Sopenharmony_ci m = BigFloat[1 1 2; 3 5 8; 13 34 21] 939bbbf1280Sopenharmony_ci v = BigFloat[1, -1, 2] 940bbbf1280Sopenharmony_ci ret = m \ v 941bbbf1280Sopenharmony_ci println(" ",repr(ret)) 942bbbf1280Sopenharmony_ci test(approx_eq(ret[1], 109/26)) 943bbbf1280Sopenharmony_ci test(approx_eq(ret[2], -105/130)) 944bbbf1280Sopenharmony_ci test(approx_eq(ret[3], -31/26)) 945bbbf1280Sopenharmony_ci 946bbbf1280Sopenharmony_ci # Test leastsquares rational functions. 947bbbf1280Sopenharmony_ci println("Leastsquares test 1:") 948bbbf1280Sopenharmony_ci n = 10000 949bbbf1280Sopenharmony_ci a = array1d(BigFloat, n+1) 950bbbf1280Sopenharmony_ci for i = 0:1:n 951bbbf1280Sopenharmony_ci a[1+i] = i/BigFloat(n) 952bbbf1280Sopenharmony_ci end 953bbbf1280Sopenharmony_ci (nc, dc) = ratfn_leastsquares(x->exp(x), a, 2, 2) 954bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 955bbbf1280Sopenharmony_ci for x = a 956bbbf1280Sopenharmony_ci test(approx_eq(exp(x), ratfn_eval(nc, dc, x), 1e-4)) 957bbbf1280Sopenharmony_ci end 958bbbf1280Sopenharmony_ci 959bbbf1280Sopenharmony_ci # Test golden section search. 960bbbf1280Sopenharmony_ci println("Golden section test 1:") 961bbbf1280Sopenharmony_ci x, y = goldensection(x->sin(x), 962bbbf1280Sopenharmony_ci BigFloat(0), BigFloat(1)/10, BigFloat(4)) 963bbbf1280Sopenharmony_ci println(" ", x, " -> ", y) 964bbbf1280Sopenharmony_ci test(approx_eq(x, asin(BigFloat(1)))) 965bbbf1280Sopenharmony_ci test(approx_eq(y, 1)) 966bbbf1280Sopenharmony_ci 967bbbf1280Sopenharmony_ci # Test extrema-winnowing algorithm. 968bbbf1280Sopenharmony_ci println("Winnow test 1:") 969bbbf1280Sopenharmony_ci extrema = [(x, sin(20*x)*sin(197*x)) 970bbbf1280Sopenharmony_ci for x in BigFloat(0):BigFloat(1)/1000:BigFloat(1)] 971bbbf1280Sopenharmony_ci winnowed = winnow_extrema(extrema, 12) 972bbbf1280Sopenharmony_ci println(" ret = ", format_xylist(winnowed)) 973bbbf1280Sopenharmony_ci prevx, prevy = -1, 0 974bbbf1280Sopenharmony_ci for (x,y) = winnowed 975bbbf1280Sopenharmony_ci test(x > prevx) 976bbbf1280Sopenharmony_ci test(y != 0) 977bbbf1280Sopenharmony_ci test(prevy * y <= 0) # tolerates initial prevx having no sign 978bbbf1280Sopenharmony_ci test(abs(y) > 0.9) 979bbbf1280Sopenharmony_ci prevx, prevy = x, y 980bbbf1280Sopenharmony_ci end 981bbbf1280Sopenharmony_ci 982bbbf1280Sopenharmony_ci # Test actual minimax approximation. 983bbbf1280Sopenharmony_ci println("Minimax test 1 (polynomial):") 984bbbf1280Sopenharmony_ci (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0) 985bbbf1280Sopenharmony_ci println(" ",e) 986bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 987bbbf1280Sopenharmony_ci test(0 < e < 1e-3) 988bbbf1280Sopenharmony_ci for x = 0:BigFloat(1)/1000:1 989bbbf1280Sopenharmony_ci test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001) 990bbbf1280Sopenharmony_ci end 991bbbf1280Sopenharmony_ci 992bbbf1280Sopenharmony_ci println("Minimax test 2 (rational):") 993bbbf1280Sopenharmony_ci (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2) 994bbbf1280Sopenharmony_ci println(" ",e) 995bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 996bbbf1280Sopenharmony_ci test(0 < e < 1e-3) 997bbbf1280Sopenharmony_ci for x = 0:BigFloat(1)/1000:1 998bbbf1280Sopenharmony_ci test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001) 999bbbf1280Sopenharmony_ci end 1000bbbf1280Sopenharmony_ci 1001bbbf1280Sopenharmony_ci println("Minimax test 3 (polynomial, weighted):") 1002bbbf1280Sopenharmony_ci (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0, 1003bbbf1280Sopenharmony_ci (x,y)->1/y) 1004bbbf1280Sopenharmony_ci println(" ",e) 1005bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 1006bbbf1280Sopenharmony_ci test(0 < e < 1e-3) 1007bbbf1280Sopenharmony_ci for x = 0:BigFloat(1)/1000:1 1008bbbf1280Sopenharmony_ci test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001) 1009bbbf1280Sopenharmony_ci end 1010bbbf1280Sopenharmony_ci 1011bbbf1280Sopenharmony_ci println("Minimax test 4 (rational, weighted):") 1012bbbf1280Sopenharmony_ci (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2, 1013bbbf1280Sopenharmony_ci (x,y)->1/y) 1014bbbf1280Sopenharmony_ci println(" ",e) 1015bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 1016bbbf1280Sopenharmony_ci test(0 < e < 1e-3) 1017bbbf1280Sopenharmony_ci for x = 0:BigFloat(1)/1000:1 1018bbbf1280Sopenharmony_ci test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001) 1019bbbf1280Sopenharmony_ci end 1020bbbf1280Sopenharmony_ci 1021bbbf1280Sopenharmony_ci println("Minimax test 5 (rational, weighted, odd degree):") 1022bbbf1280Sopenharmony_ci (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 1, 1023bbbf1280Sopenharmony_ci (x,y)->1/y) 1024bbbf1280Sopenharmony_ci println(" ",e) 1025bbbf1280Sopenharmony_ci println(" ",ratfn_to_string(nc, dc)) 1026bbbf1280Sopenharmony_ci test(0 < e < 1e-3) 1027bbbf1280Sopenharmony_ci for x = 0:BigFloat(1)/1000:1 1028bbbf1280Sopenharmony_ci test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001) 1029bbbf1280Sopenharmony_ci end 1030bbbf1280Sopenharmony_ci 1031bbbf1280Sopenharmony_ci total = passes + fails 1032bbbf1280Sopenharmony_ci println(passes, " passes ", fails, " fails ", total, " total") 1033bbbf1280Sopenharmony_ciend 1034bbbf1280Sopenharmony_ci 1035bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 1036bbbf1280Sopenharmony_ci# Online help. 1037bbbf1280Sopenharmony_cifunction help() 1038bbbf1280Sopenharmony_ci print(""" 1039bbbf1280Sopenharmony_ciUsage: 1040bbbf1280Sopenharmony_ci 1041bbbf1280Sopenharmony_ci remez.jl [options] <lo> <hi> <n> <d> <expr> [<weight>] 1042bbbf1280Sopenharmony_ci 1043bbbf1280Sopenharmony_ciArguments: 1044bbbf1280Sopenharmony_ci 1045bbbf1280Sopenharmony_ci <lo>, <hi> 1046bbbf1280Sopenharmony_ci 1047bbbf1280Sopenharmony_ci Bounds of the interval on which to approximate the target 1048bbbf1280Sopenharmony_ci function. These are parsed and evaluated as Julia expressions, 1049bbbf1280Sopenharmony_ci so you can write things like '1/BigFloat(6)' to get an 1050bbbf1280Sopenharmony_ci accurate representation of 1/6, or '4*atan(BigFloat(1))' to 1051bbbf1280Sopenharmony_ci get pi. (Unfortunately, the obvious 'BigFloat(pi)' doesn't 1052bbbf1280Sopenharmony_ci work in Julia.) 1053bbbf1280Sopenharmony_ci 1054bbbf1280Sopenharmony_ci <n>, <d> 1055bbbf1280Sopenharmony_ci 1056bbbf1280Sopenharmony_ci The desired degree of polynomial(s) you want for your 1057bbbf1280Sopenharmony_ci approximation. These should be non-negative integers. If you 1058bbbf1280Sopenharmony_ci want a rational function as output, set <n> to the degree of 1059bbbf1280Sopenharmony_ci the numerator, and <d> the denominator. If you just want an 1060bbbf1280Sopenharmony_ci ordinary polynomial, set <d> to 0, and <n> to the degree of 1061bbbf1280Sopenharmony_ci the polynomial you want. 1062bbbf1280Sopenharmony_ci 1063bbbf1280Sopenharmony_ci <expr> 1064bbbf1280Sopenharmony_ci 1065bbbf1280Sopenharmony_ci A Julia expression giving the function to be approximated on 1066bbbf1280Sopenharmony_ci the interval. The input value is predefined as 'x' when this 1067bbbf1280Sopenharmony_ci expression is evaluated, so you should write something along 1068bbbf1280Sopenharmony_ci the lines of 'sin(x)' or 'sqrt(1+tan(x)^2)' etc. 1069bbbf1280Sopenharmony_ci 1070bbbf1280Sopenharmony_ci <weight> 1071bbbf1280Sopenharmony_ci 1072bbbf1280Sopenharmony_ci If provided, a Julia expression giving the weighting factor 1073bbbf1280Sopenharmony_ci for the approximation error. The output polynomial will 1074bbbf1280Sopenharmony_ci minimise the largest absolute value of (P-f) * w at any point 1075bbbf1280Sopenharmony_ci in the interval, where P is the value of the polynomial, f is 1076bbbf1280Sopenharmony_ci the value of the target function given by <expr>, and w is the 1077bbbf1280Sopenharmony_ci weight given by this function. 1078bbbf1280Sopenharmony_ci 1079bbbf1280Sopenharmony_ci When this expression is evaluated, the input value to P and f 1080bbbf1280Sopenharmony_ci is predefined as 'x', and also the true output value f(x) is 1081bbbf1280Sopenharmony_ci predefined as 'y'. So you can minimise the relative error by 1082bbbf1280Sopenharmony_ci simply writing '1/y'. 1083bbbf1280Sopenharmony_ci 1084bbbf1280Sopenharmony_ci If the <weight> argument is not provided, the default 1085bbbf1280Sopenharmony_ci weighting function always returns 1, so that the polynomial 1086bbbf1280Sopenharmony_ci will minimise the maximum absolute error |P-f|. 1087bbbf1280Sopenharmony_ci 1088bbbf1280Sopenharmony_ciComputation options: 1089bbbf1280Sopenharmony_ci 1090bbbf1280Sopenharmony_ci --pre=<predef_expr> 1091bbbf1280Sopenharmony_ci 1092bbbf1280Sopenharmony_ci Evaluate the Julia expression <predef_expr> before starting 1093bbbf1280Sopenharmony_ci the computation. This permits you to pre-define variables or 1094bbbf1280Sopenharmony_ci functions which the Julia expressions in your main arguments 1095bbbf1280Sopenharmony_ci can refer to. All of <lo>, <hi>, <expr> and <weight> can make 1096bbbf1280Sopenharmony_ci use of things defined by <predef_expr>. 1097bbbf1280Sopenharmony_ci 1098bbbf1280Sopenharmony_ci One internal remez.jl function that you might sometimes find 1099bbbf1280Sopenharmony_ci useful in this expression is 'goldensection', which finds the 1100bbbf1280Sopenharmony_ci location and value of a maximum of a function. For example, 1101bbbf1280Sopenharmony_ci one implementation strategy for the gamma function involves 1102bbbf1280Sopenharmony_ci translating it to put its unique local minimum at the origin, 1103bbbf1280Sopenharmony_ci in which case you can write something like this 1104bbbf1280Sopenharmony_ci 1105bbbf1280Sopenharmony_ci --pre='(m,my) = goldensection(x -> -gamma(x), 1106bbbf1280Sopenharmony_ci BigFloat(1), BigFloat(1.5), BigFloat(2))' 1107bbbf1280Sopenharmony_ci 1108bbbf1280Sopenharmony_ci to predefine 'm' as the location of gamma's minimum, and 'my' 1109bbbf1280Sopenharmony_ci as the (negated) value that gamma actually takes at that 1110bbbf1280Sopenharmony_ci point, i.e. -gamma(m). 1111bbbf1280Sopenharmony_ci 1112bbbf1280Sopenharmony_ci (Since 'goldensection' always finds a maximum, we had to 1113bbbf1280Sopenharmony_ci negate gamma in the input function to make it find a minimum 1114bbbf1280Sopenharmony_ci instead. Consult the comments in the source for more details 1115bbbf1280Sopenharmony_ci on the use of this function.) 1116bbbf1280Sopenharmony_ci 1117bbbf1280Sopenharmony_ci If you use this option more than once, all the expressions you 1118bbbf1280Sopenharmony_ci provide will be run in sequence. 1119bbbf1280Sopenharmony_ci 1120bbbf1280Sopenharmony_ci --bits=<bits> 1121bbbf1280Sopenharmony_ci 1122bbbf1280Sopenharmony_ci Specify the accuracy to which you want the output polynomial, 1123bbbf1280Sopenharmony_ci in bits. Default 256, which should be more than enough. 1124bbbf1280Sopenharmony_ci 1125bbbf1280Sopenharmony_ci --bigfloatbits=<bits> 1126bbbf1280Sopenharmony_ci 1127bbbf1280Sopenharmony_ci Turn up the precision used by Julia for its BigFloat 1128bbbf1280Sopenharmony_ci evaluation. Default is Julia's default (also 256). You might 1129bbbf1280Sopenharmony_ci want to try setting this higher than the --bits value if the 1130bbbf1280Sopenharmony_ci algorithm is failing to converge for some reason. 1131bbbf1280Sopenharmony_ci 1132bbbf1280Sopenharmony_ciOutput options: 1133bbbf1280Sopenharmony_ci 1134bbbf1280Sopenharmony_ci --full 1135bbbf1280Sopenharmony_ci 1136bbbf1280Sopenharmony_ci Instead of just printing the approximation function itself, 1137bbbf1280Sopenharmony_ci also print auxiliary information: 1138bbbf1280Sopenharmony_ci - the locations of the error extrema, and the actual 1139bbbf1280Sopenharmony_ci (weighted) error at each of those locations 1140bbbf1280Sopenharmony_ci - the overall maximum error of the function 1141bbbf1280Sopenharmony_ci - a 'well-conditioning quotient', giving the worst-case ratio 1142bbbf1280Sopenharmony_ci between any polynomial coefficient and the largest possible 1143bbbf1280Sopenharmony_ci value of the higher-order terms it will be added to. 1144bbbf1280Sopenharmony_ci 1145bbbf1280Sopenharmony_ci The well-conditioning quotient should be less than 1, ideally 1146bbbf1280Sopenharmony_ci by several factors of two, for accurate evaluation in the 1147bbbf1280Sopenharmony_ci target precision. If you request a rational function, a 1148bbbf1280Sopenharmony_ci separate well-conditioning quotient will be printed for the 1149bbbf1280Sopenharmony_ci numerator and denominator. 1150bbbf1280Sopenharmony_ci 1151bbbf1280Sopenharmony_ci Use this option when deciding how wide an interval to 1152bbbf1280Sopenharmony_ci approximate your function on, and what degree of polynomial 1153bbbf1280Sopenharmony_ci you need. 1154bbbf1280Sopenharmony_ci 1155bbbf1280Sopenharmony_ci --variable=<identifier> 1156bbbf1280Sopenharmony_ci 1157bbbf1280Sopenharmony_ci When writing the output polynomial or rational function in its 1158bbbf1280Sopenharmony_ci usual form as an arithmetic expression, use <identifier> as 1159bbbf1280Sopenharmony_ci the name of the input variable. Default is 'x'. 1160bbbf1280Sopenharmony_ci 1161bbbf1280Sopenharmony_ci --suffix=<suffix> 1162bbbf1280Sopenharmony_ci 1163bbbf1280Sopenharmony_ci When writing the output polynomial or rational function in its 1164bbbf1280Sopenharmony_ci usual form as an arithmetic expression, write <suffix> after 1165bbbf1280Sopenharmony_ci every floating-point literal. For example, '--suffix=F' will 1166bbbf1280Sopenharmony_ci generate a C expression in which the coefficients are literals 1167bbbf1280Sopenharmony_ci of type 'float' rather than 'double'. 1168bbbf1280Sopenharmony_ci 1169bbbf1280Sopenharmony_ci --array 1170bbbf1280Sopenharmony_ci 1171bbbf1280Sopenharmony_ci Instead of writing the output polynomial as an arithmetic 1172bbbf1280Sopenharmony_ci expression in Horner's rule form, write out just its 1173bbbf1280Sopenharmony_ci coefficients, one per line, each with a trailing comma. 1174bbbf1280Sopenharmony_ci Suitable for pasting into a C array declaration. 1175bbbf1280Sopenharmony_ci 1176bbbf1280Sopenharmony_ci This option is not currently supported if the output is a 1177bbbf1280Sopenharmony_ci rational function, because you'd need two separate arrays for 1178bbbf1280Sopenharmony_ci the numerator and denominator coefficients and there's no 1179bbbf1280Sopenharmony_ci obviously right way to provide both of those together. 1180bbbf1280Sopenharmony_ci 1181bbbf1280Sopenharmony_ciDebug and test options: 1182bbbf1280Sopenharmony_ci 1183bbbf1280Sopenharmony_ci --debug=<facility> 1184bbbf1280Sopenharmony_ci 1185bbbf1280Sopenharmony_ci Enable debugging output from various parts of the Remez 1186bbbf1280Sopenharmony_ci calculation. <facility> should be the name of one of the 1187bbbf1280Sopenharmony_ci classes of diagnostic output implemented in the program. 1188bbbf1280Sopenharmony_ci Useful values include 'gausselim', 'leastsquares', 1189bbbf1280Sopenharmony_ci 'goldensection', 'equaldev', 'minimax'. This is probably 1190bbbf1280Sopenharmony_ci mostly useful to people debugging problems with the script, so 1191bbbf1280Sopenharmony_ci consult the source code for more information about what the 1192bbbf1280Sopenharmony_ci diagnostic output for each of those facilities will be. 1193bbbf1280Sopenharmony_ci 1194bbbf1280Sopenharmony_ci If you want diagnostics from more than one facility, specify 1195bbbf1280Sopenharmony_ci this option multiple times with different arguments. 1196bbbf1280Sopenharmony_ci 1197bbbf1280Sopenharmony_ci --test 1198bbbf1280Sopenharmony_ci 1199bbbf1280Sopenharmony_ci Run remez.jl's internal test suite. No arguments needed. 1200bbbf1280Sopenharmony_ci 1201bbbf1280Sopenharmony_ciMiscellaneous options: 1202bbbf1280Sopenharmony_ci 1203bbbf1280Sopenharmony_ci --help 1204bbbf1280Sopenharmony_ci 1205bbbf1280Sopenharmony_ci Display this text and exit. No arguments needed. 1206bbbf1280Sopenharmony_ci 1207bbbf1280Sopenharmony_ci""") 1208bbbf1280Sopenharmony_ciend 1209bbbf1280Sopenharmony_ci 1210bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 1211bbbf1280Sopenharmony_ci# Main program. 1212bbbf1280Sopenharmony_ci 1213bbbf1280Sopenharmony_cifunction main() 1214bbbf1280Sopenharmony_ci nargs = length(argwords) 1215bbbf1280Sopenharmony_ci if nargs != 5 && nargs != 6 1216bbbf1280Sopenharmony_ci error("usage: remez.jl <lo> <hi> <n> <d> <expr> [<weight>]\n" * 1217bbbf1280Sopenharmony_ci " run 'remez.jl --help' for more help") 1218bbbf1280Sopenharmony_ci end 1219bbbf1280Sopenharmony_ci 1220bbbf1280Sopenharmony_ci for preliminary_command in preliminary_commands 1221bbbf1280Sopenharmony_ci eval(Meta.parse(preliminary_command)) 1222bbbf1280Sopenharmony_ci end 1223bbbf1280Sopenharmony_ci 1224bbbf1280Sopenharmony_ci lo = BigFloat(eval(Meta.parse(argwords[1]))) 1225bbbf1280Sopenharmony_ci hi = BigFloat(eval(Meta.parse(argwords[2]))) 1226bbbf1280Sopenharmony_ci n = parse(Int,argwords[3]) 1227bbbf1280Sopenharmony_ci d = parse(Int,argwords[4]) 1228bbbf1280Sopenharmony_ci f = eval(Meta.parse("x -> " * argwords[5])) 1229bbbf1280Sopenharmony_ci 1230bbbf1280Sopenharmony_ci # Wrap the user-provided function with a function of our own. This 1231bbbf1280Sopenharmony_ci # arranges to detect silly FP values (inf,nan) early and diagnose 1232bbbf1280Sopenharmony_ci # them sensibly, and also lets us log all evaluations of the 1233bbbf1280Sopenharmony_ci # function in case you suspect it's doing the wrong thing at some 1234bbbf1280Sopenharmony_ci # special-case point. 1235bbbf1280Sopenharmony_ci function func(x) 1236bbbf1280Sopenharmony_ci y = run(f,x) 1237bbbf1280Sopenharmony_ci @debug("f", x, " -> ", y) 1238bbbf1280Sopenharmony_ci if !isfinite(y) 1239bbbf1280Sopenharmony_ci error("f(" * string(x) * ") returned non-finite value " * string(y)) 1240bbbf1280Sopenharmony_ci end 1241bbbf1280Sopenharmony_ci return y 1242bbbf1280Sopenharmony_ci end 1243bbbf1280Sopenharmony_ci 1244bbbf1280Sopenharmony_ci if nargs == 6 1245bbbf1280Sopenharmony_ci # Wrap the user-provided weight function similarly. 1246bbbf1280Sopenharmony_ci w = eval(Meta.parse("(x,y) -> " * argwords[6])) 1247bbbf1280Sopenharmony_ci function wrapped_weight(x,y) 1248bbbf1280Sopenharmony_ci ww = run(w,x,y) 1249bbbf1280Sopenharmony_ci if !isfinite(ww) 1250bbbf1280Sopenharmony_ci error("w(" * string(x) * "," * string(y) * 1251bbbf1280Sopenharmony_ci ") returned non-finite value " * string(ww)) 1252bbbf1280Sopenharmony_ci end 1253bbbf1280Sopenharmony_ci return ww 1254bbbf1280Sopenharmony_ci end 1255bbbf1280Sopenharmony_ci weight = wrapped_weight 1256bbbf1280Sopenharmony_ci else 1257bbbf1280Sopenharmony_ci weight = (x,y)->BigFloat(1) 1258bbbf1280Sopenharmony_ci end 1259bbbf1280Sopenharmony_ci 1260bbbf1280Sopenharmony_ci (nc, dc, e, extrema) = ratfn_minimax(func, (lo, hi), n, d, weight) 1261bbbf1280Sopenharmony_ci if array_format 1262bbbf1280Sopenharmony_ci if d == 0 1263bbbf1280Sopenharmony_ci functext = join([string(x)*",\n" for x=nc],"") 1264bbbf1280Sopenharmony_ci else 1265bbbf1280Sopenharmony_ci # It's unclear how you should best format an array of 1266bbbf1280Sopenharmony_ci # coefficients for a rational function, so I'll leave 1267bbbf1280Sopenharmony_ci # implementing this option until I have a use case. 1268bbbf1280Sopenharmony_ci error("--array unsupported for rational functions") 1269bbbf1280Sopenharmony_ci end 1270bbbf1280Sopenharmony_ci else 1271bbbf1280Sopenharmony_ci functext = ratfn_to_string(nc, dc) * "\n" 1272bbbf1280Sopenharmony_ci end 1273bbbf1280Sopenharmony_ci if full_output 1274bbbf1280Sopenharmony_ci # Print everything you might want to know about the function 1275bbbf1280Sopenharmony_ci println("extrema = ", format_xylist(extrema)) 1276bbbf1280Sopenharmony_ci println("maxerror = ", string(e)) 1277bbbf1280Sopenharmony_ci if length(dc) > 1 1278bbbf1280Sopenharmony_ci println("wellconditioning_numerator = ", 1279bbbf1280Sopenharmony_ci string(wellcond(nc, lo, hi))) 1280bbbf1280Sopenharmony_ci println("wellconditioning_denominator = ", 1281bbbf1280Sopenharmony_ci string(wellcond(dc, lo, hi))) 1282bbbf1280Sopenharmony_ci else 1283bbbf1280Sopenharmony_ci println("wellconditioning = ", string(wellcond(nc, lo, hi))) 1284bbbf1280Sopenharmony_ci end 1285bbbf1280Sopenharmony_ci print("function = ", functext) 1286bbbf1280Sopenharmony_ci else 1287bbbf1280Sopenharmony_ci # Just print the text people will want to paste into their code 1288bbbf1280Sopenharmony_ci print(functext) 1289bbbf1280Sopenharmony_ci end 1290bbbf1280Sopenharmony_ciend 1291bbbf1280Sopenharmony_ci 1292bbbf1280Sopenharmony_ci# ---------------------------------------------------------------------- 1293bbbf1280Sopenharmony_ci# Top-level code: parse the argument list and decide what to do. 1294bbbf1280Sopenharmony_ci 1295bbbf1280Sopenharmony_ciwhat_to_do = main 1296bbbf1280Sopenharmony_ci 1297bbbf1280Sopenharmony_cidoing_opts = true 1298bbbf1280Sopenharmony_ciargwords = array1d(String, 0) 1299bbbf1280Sopenharmony_cifor arg = ARGS 1300bbbf1280Sopenharmony_ci global doing_opts, what_to_do, argwords 1301bbbf1280Sopenharmony_ci global full_output, array_format, xvarname, floatsuffix, epsbits 1302bbbf1280Sopenharmony_ci if doing_opts && startswith(arg, "-") 1303bbbf1280Sopenharmony_ci if arg == "--" 1304bbbf1280Sopenharmony_ci doing_opts = false 1305bbbf1280Sopenharmony_ci elseif arg == "--help" 1306bbbf1280Sopenharmony_ci what_to_do = help 1307bbbf1280Sopenharmony_ci elseif arg == "--test" 1308bbbf1280Sopenharmony_ci what_to_do = test 1309bbbf1280Sopenharmony_ci elseif arg == "--full" 1310bbbf1280Sopenharmony_ci full_output = true 1311bbbf1280Sopenharmony_ci elseif arg == "--array" 1312bbbf1280Sopenharmony_ci array_format = true 1313bbbf1280Sopenharmony_ci elseif startswith(arg, "--debug=") 1314bbbf1280Sopenharmony_ci enable_debug(arg[length("--debug=")+1:end]) 1315bbbf1280Sopenharmony_ci elseif startswith(arg, "--variable=") 1316bbbf1280Sopenharmony_ci xvarname = arg[length("--variable=")+1:end] 1317bbbf1280Sopenharmony_ci elseif startswith(arg, "--suffix=") 1318bbbf1280Sopenharmony_ci floatsuffix = arg[length("--suffix=")+1:end] 1319bbbf1280Sopenharmony_ci elseif startswith(arg, "--bits=") 1320bbbf1280Sopenharmony_ci epsbits = parse(Int,arg[length("--bits=")+1:end]) 1321bbbf1280Sopenharmony_ci elseif startswith(arg, "--bigfloatbits=") 1322bbbf1280Sopenharmony_ci set_bigfloat_precision( 1323bbbf1280Sopenharmony_ci parse(Int,arg[length("--bigfloatbits=")+1:end])) 1324bbbf1280Sopenharmony_ci elseif startswith(arg, "--pre=") 1325bbbf1280Sopenharmony_ci push!(preliminary_commands, arg[length("--pre=")+1:end]) 1326bbbf1280Sopenharmony_ci else 1327bbbf1280Sopenharmony_ci error("unrecognised option: ", arg) 1328bbbf1280Sopenharmony_ci end 1329bbbf1280Sopenharmony_ci else 1330bbbf1280Sopenharmony_ci push!(argwords, arg) 1331bbbf1280Sopenharmony_ci end 1332bbbf1280Sopenharmony_ciend 1333bbbf1280Sopenharmony_ci 1334bbbf1280Sopenharmony_ciwhat_to_do() 1335