Lines Matching +full:two +full:- +full:dimensional
2 # -*- julia -*-
4 # remez.jl - implementation of the Remez algorithm for polynomial approximation
6 # Copyright (c) 2015-2019, Arm Limited.
7 # SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
11 # ----------------------------------------------------------------------
33 # ----------------------------------------------------------------------
34 # Global variables configured by command-line options.
35 floatsuffix = "" # adjusted by --floatsuffix
36 xvarname = "x" # adjusted by --variable
37 epsbits = 256 # adjusted by --bits
38 debug_facilities = Set() # adjusted by --debug
39 full_output = false # adjusted by --full
40 array_format = false # adjusted by --array
41 preliminary_commands = array1d(String, 0) # adjusted by --pre
43 # ----------------------------------------------------------------------
87 # coefficient of x^(i-1) (because Julia arrays are 1-based).
106 # Starts with the constant term, and 1-based, as above.
108 # Starts with the constant term, and 1-based, as above.
131 # Starts with the constant term, and 1-based, as above.
150 # Starts with the constant term, and 1-based, as above.
152 # Starts with the constant term, and 1-based, as above.
173 join([" "*string(x)*" -> "*string(y) for (x,y) in xys], "\n") *
177 # ----------------------------------------------------------------------
178 # Matrix-equation solver for matrices of BigFloat.
180 # I had hoped that Julia's type-genericity would allow me to solve the
193 # matrix_in 2-dimensional array of BigFloats, representing a matrix M
194 # in row-first order, i.e. matrix_in[r,c] represents the
196 # vector_in 1-dimensional array of BigFloats, representing a vector V.
198 # Return value: a 1-dimensional array X of BigFloats, satisfying M X = V.
218 # non-zero entry in column i (and in a row we haven't sorted
261 M[j,k] = M[j,k] - M[i,k] * factor
263 V[j] = V[j] - V[i] * factor
279 # ----------------------------------------------------------------------
280 # Least-squares fitting of a rational function to a set of (x,y)
290 # Least-squares fitting of a _polynomial_ is actually a sensible thing
293 # minimise the error function (P/Q-f)^2 that you actually wanted;
294 # instead it minimises (P-fQ)^2. But that should be good enough to
309 # E = \sum (p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d)^2
314 # 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) x^j
318 # 0 = \sum 2(p_0 + p_1 x + ... + p_n x^n - y - y q_1 x - ... - y q_d x^d) y x^j
329 # f The function to be approximated. Maps BigFloat -> BigFloat.
330 # xvals Array of BigFloats, giving the list of x-coordinates at which
336 # w Error-weighting function. Takes two BigFloat arguments x,y
345 # coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
348 w = (x,y)->BigFloat(1))
350 # Again because Julia arrays are 1-based, we'll have sums[i,j]
351 # being the sum of x^(i-1) y^(j-1).
359 sums[i,j] += x^(i-1) * y^(j-1) * weight
379 matrix[row, 1+n+j] = -sums[1+i+j, 2]
391 matrix[row, 1+n+j] = -sums[1+i+j, 3]
404 # And marshal the results into two separate polynomial vectors to
411 # ----------------------------------------------------------------------
412 # Golden-section search to find a maximum of a function.
415 # f Function to be maximised/minimised. Maps BigFloat -> BigFloat.
427 threshold = abs(c-a) * 2^(-epsbits/2)
436 if abs(b-a) < abs(c-a)
447 while abs(c-a) > threshold
448 @debug("goldensection", "a: ", a, " -> ", fa)
449 @debug("goldensection", "b: ", b, " -> ", fb)
450 @debug("goldensection", "c: ", c, " -> ", fc)
458 d = a + (b-a) * one_over_phi
466 # two cases: either f(d) > f(b), or vice versa. In either
471 @debug("goldensection", "d: ", d, " -> ", fd)
483 @debug("goldensection", "done: ", b, " -> ", fb)
487 # ----------------------------------------------------------------------
491 # f The function to be approximated. Maps BigFloat -> BigFloat.
503 # the interval, then the only way we can do that is to set two
505 prev = max(1, i-1)
517 x, y = goldensection(x->-f(x), xp, xi, xn)
518 push!(extrema, (x, -y))
524 # ----------------------------------------------------------------------
530 # We do this using a dynamic-programming approach. We work along the
556 # A zero-valued extremum cannot possibly contribute to any
565 # Special case: we don't store a best zero-length
566 # sequence :-)
569 othersign = 3-sign # map 1->2 and 2->1
570 oldscore, oldlist = best[i-1, othersign]
596 # ----------------------------------------------------------------------
597 # Construct a rational-function approximation with equal and
598 # alternating weighted deviation at a specific set of x-coordinates.
601 # f The function to be approximated. Maps BigFloat -> BigFloat.
602 # coords An array of BigFloats giving the x-coordinates. There should
609 # w Error-weighting function. Takes two BigFloat arguments x,y
612 # maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
613 # parameter, defaulting to the always-return-1 function.
618 # coefficient of x^(i-1) (because Julia arrays are 1-based). D[1] will
622 w = (x,y)->BigFloat(1))
650 currsign = -currsign
669 # --------------------------------- ± e/w(x) = f(x)
683 # equation will be anti-monotonic in the input e, so we can
684 # just binary-search until we get the two as close to equal as
699 y_adj = y - currsign * e / w(x,y)
704 matrix[i,1+n+j] = -x^j * y_adj
707 currsign = -currsign
723 last_e = (ratfn_eval(ncoeffs, dcoeffs, x) - y) * w(x,y) * -currsign
730 threshold = 2^(-epsbits/2) # convergence threshold
733 # value (e0) will be one end of our binary-search interval,
740 if abs(e1-e0) <= threshold
745 s = sign(e1-e0)
754 @assert(sign(e2-e1) == -s)
756 # Now binary-search until our two endpoints narrow enough.
758 while abs(e1-e0) > threshold
761 if sign(enew-emid) == s
773 # ----------------------------------------------------------------------
774 # Top-level function to find a minimax rational-function approximation.
777 # f The function to be approximated. Maps BigFloat -> BigFloat.
782 # w Error-weighting function. Takes two BigFloat arguments x,y
785 # maximum value of abs((f(x)-R(x)) * w(x,f(x))). Optional
786 # parameter, defaulting to the always-return-1 function.
793 # N[i] is the coefficient of x^(i-1) (because Julia
794 # arrays are 1-based). D[1] will be 1.
803 w = (x,y)->BigFloat(1))
804 # We start off by finding a least-squares approximation. This
808 # Least-squares approximations tend to look nicer in a minimax
819 halfwid = (hi-lo)/2
822 grid = [ mid - halfwid * cos(pi*i/nnodes) for i=0:1:nnodes ]
825 # Find the initial least-squares approximation.
835 # less empirical basis, because the rational-function case will
837 # (Repeatable by using the --test mode.) On the assumption that
842 threshold = 2^(-epsbits/3)
850 return (approx_y - real_y) * w(x, real_y)
863 variation = (max_err - min_err) / max_err
872 (nc, dc) = ratfn_equal_deviation(f, map(x->x[1], extrema),
878 # ----------------------------------------------------------------------
879 # Check if a polynomial is well-conditioned for accurate evaluation in
890 # x^(i-1) (because Julia arrays are 1-based).
895 # than 1, the polynomial is at least somewhat well-conditioned;
904 for i = length(coeffs):-1:1
916 # ----------------------------------------------------------------------
923 function approx_eq(x, y, limit=1e-6)
924 return abs(x - y) < limit
939 v = BigFloat[1, -1, 2]
943 test(approx_eq(ret[2], -105/130))
944 test(approx_eq(ret[3], -31/26))
953 (nc, dc) = ratfn_leastsquares(x->exp(x), a, 2, 2)
956 test(approx_eq(exp(x), ratfn_eval(nc, dc, x), 1e-4))
961 x, y = goldensection(x->sin(x),
963 println(" ", x, " -> ", y)
967 # Test extrema-winnowing algorithm.
973 prevx, prevy = -1, 0
984 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0)
987 test(0 < e < 1e-3)
989 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)
993 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2)
996 test(0 < e < 1e-3)
998 test(abs(ratfn_eval(nc, dc, x) - exp(x)) <= e * 1.0000001)
1002 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 4, 0,
1003 (x,y)->1/y)
1006 test(0 < e < 1e-3)
1008 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1012 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 2,
1013 (x,y)->1/y)
1016 test(0 < e < 1e-3)
1018 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1022 (nc, dc, e, x) = ratfn_minimax(x->exp(x), (BigFloat(0), BigFloat(1)), 2, 1,
1023 (x,y)->1/y)
1026 test(0 < e < 1e-3)
1028 test(abs(ratfn_eval(nc, dc, x) - exp(x))/exp(x) <= e * 1.0000001)
1035 # ----------------------------------------------------------------------
1057 approximation. These should be non-negative integers. If you
1074 minimise the largest absolute value of (P-f) * w at any point
1086 will minimise the maximum absolute error |P-f|.
1090 --pre=<predef_expr>
1093 the computation. This permits you to pre-define variables or
1105 --pre='(m,my) = goldensection(x -> -gamma(x),
1110 point, i.e. -gamma(m).
1120 --bits=<bits>
1125 --bigfloatbits=<bits>
1129 want to try setting this higher than the --bits value if the
1134 --full
1138 - the locations of the error extrema, and the actual
1140 - the overall maximum error of the function
1141 - a 'well-conditioning quotient', giving the worst-case ratio
1143 value of the higher-order terms it will be added to.
1145 The well-conditioning quotient should be less than 1, ideally
1146 by several factors of two, for accurate evaluation in the
1148 separate well-conditioning quotient will be printed for the
1155 --variable=<identifier>
1161 --suffix=<suffix>
1165 every floating-point literal. For example, '--suffix=F' will
1169 --array
1177 rational function, because you'd need two separate arrays for
1183 --debug=<facility>
1197 --test
1203 --help
1210 # ----------------------------------------------------------------------
1217 " run 'remez.jl --help' for more help")
1228 f = eval(Meta.parse("x -> " * argwords[5]))
1230 # Wrap the user-provided function with a function of our own. This
1234 # special-case point.
1237 @debug("f", x, " -> ", y)
1239 error("f(" * string(x) * ") returned non-finite value " * string(y))
1245 # Wrap the user-provided weight function similarly.
1246 w = eval(Meta.parse("(x,y) -> " * argwords[6]))
1251 ") returned non-finite value " * string(ww))
1257 weight = (x,y)->BigFloat(1)
1268 error("--array unsupported for rational functions")
1292 # ----------------------------------------------------------------------
1293 # Top-level code: parse the argument list and decide what to do.
1302 if doing_opts && startswith(arg, "-")
1303 if arg == "--"
1305 elseif arg == "--help"
1307 elseif arg == "--test"
1309 elseif arg == "--full"
1311 elseif arg == "--array"
1313 elseif startswith(arg, "--debug=")
1314 enable_debug(arg[length("--debug=")+1:end])
1315 elseif startswith(arg, "--variable=")
1316 xvarname = arg[length("--variable=")+1:end]
1317 elseif startswith(arg, "--suffix=")
1318 floatsuffix = arg[length("--suffix=")+1:end]
1319 elseif startswith(arg, "--bits=")
1320 epsbits = parse(Int,arg[length("--bits=")+1:end])
1321 elseif startswith(arg, "--bigfloatbits=")
1323 parse(Int,arg[length("--bigfloatbits=")+1:end]))
1324 elseif startswith(arg, "--pre=")
1325 push!(preliminary_commands, arg[length("--pre=")+1:end])