1[section:minimax Minimax Approximations and the Remez Algorithm] 2 3The directory `libs/math/minimax` contains an interactive command-line driven 4program for the generation of minimax approximations using the Remez 5algorithm. Both polynomial and rational approximations are supported, 6although the latter are tricky to converge: it is not uncommon for 7convergence of rational forms to fail. No such limitations are present 8for polynomial approximations which should always converge smoothly. 9 10It's worth stressing that developing rational approximations to functions 11is often not an easy task, and one to which many books have been devoted. 12To use this tool, you will need to have a reasonable grasp of what the Remez 13algorithm is, and the general form of the approximation you want to achieve. 14 15Unless you already familiar with the Remez method, you should first read the 16[link math_toolkit.remez brief background article explaining the principles behind the Remez algorithm]. 17 18The program consists of two parts: 19 20[variablelist 21[[main.cpp][Contains the command line parser, and all the calls to the Remez code.]] 22[[f.cpp][Contains the function to approximate.]] 23] 24 25Therefore to use this tool, you must modify f.cpp to return the function to 26approximate. The tools supports multiple function approximations within 27the same compiled program: each as a separate variant: 28 29 NTL::RR f(const NTL::RR& x, int variant); 30 31Returns the value of the function /variant/ at point /x/. So if you 32wish you can just add the function to approximate as a new variant 33after the existing examples. 34 35In addition to those two files, the program needs to be linked to 36a [link math_toolkit.high_precision.use_ntl patched NTL library to compile]. 37 38Note that the function /f/ must return the rational part of the 39approximation: for example if you are approximating a function 40/f(x)/ then it is quite common to use: 41 42[expression f(x) = g(x)(Y + R(x))] 43 44where /g(x)/ is the dominant part of /f(x)/, /Y/ is some constant, and 45/R(x)/ is the rational approximation part, usually optimised for a low 46absolute error compared to |Y|. 47 48In this case you would define /f/ to return [role serif-italic f(x)/g(x)] and then set the 49y-offset of the approximation to /Y/ (see command line options below). 50 51Many other forms are possible, but in all cases the objective is to 52split /f(x)/ into a dominant part that you can evaluate easily using 53standard math functions, and a smooth and slowly changing rational approximation 54part. Refer to your favourite textbook for more examples. 55 56Command line options for the program are as follows: 57 58[variablelist 59[[variant N][Sets the current function variant to N. This allows multiple functions 60 that are to be approximated to be compiled into the same executable. 61 Defaults to 0.]] 62[[range a b][Sets the domain for the approximation to the range \[a,b\], defaults 63 to \[0,1\].]] 64[[relative][Sets the Remez code to optimise for relative error. This is the default 65 at program startup. Note that relative error can only be used 66 if f(x) has no roots over the range being optimised.]] 67[[absolute][Sets the Remez code to optimise for absolute error.]] 68[[pin \[true|false\]]["Pins" the code so that the rational approximation 69 passes through the origin. Obviously only set this to 70 /true/ if R(0) must be zero. This is typically used when 71 trying to preserve a root at \[0,0\] while also optimising 72 for relative error.]] 73[[order N D][Sets the order of the approximation to /N/ in the numerator and /D/ 74 in the denominator. If /D/ is zero then the result will be a polynomial 75 approximation. There will be N+D+2 coefficients in total, the first 76 coefficient of the numerator is zero if /pin/ was set to true, and the 77 first coefficient of the denominator is always one.]] 78[[working-precision N][Sets the working precision of NTL::RR to /N/ binary digits. Defaults to 250.]] 79[[target-precision N][Sets the precision of printed output to /N/ binary digits: 80 set to the same number of digits as the type that will be used to 81 evaluate the approximation. Defaults to 53 (for double precision).]] 82[[skew val]["Skews" the initial interpolated control points towards one 83 end or the other of the range. Positive values skew the 84 initial control points towards the left hand side of the 85 range, and negative values towards the right hand side. 86 If an approximation won't converge (a common situation) 87 try adjusting the skew parameter until the first step yields 88 the smallest possible error. /val/ should be in the range 89 \[-100,+100\], the default is zero.]] 90[[brake val][Sets a brake on each step so that the change in the 91 control points is braked by /val%/. Defaults to 50, 92 try a higher value if an approximation won't converge, 93 or a lower value to get speedier convergence.]] 94[[x-offset val][Sets the x-offset to /val/: the approximation will 95 be generated for `f(S * (x + X)) + Y` where /X/ is the 96 x-offset, /S/ is the x-scale 97 and /Y/ is the y-offset. Defaults to zero. To avoid 98 rounding errors, take care to specify a value that can 99 be exactly represented as a floating point number.]] 100[[x-scale val][Sets the x-scale to /val/: the approximation will 101 be generated for `f(S * (x + X)) + Y` where /S/ is the 102 x-scale, /X/ is the x-offset 103 and /Y/ is the y-offset. Defaults to one. To avoid 104 rounding errors, take care to specify a value that can 105 be exactly represented as a floating point number.]] 106[[y-offset val][Sets the y-offset to /val/: the approximation will 107 be generated for `f(S * (x + X)) + Y` where /X/ 108 is the x-offset, /S/ is the x-scale 109 and /Y/ is the y-offset. Defaults to zero. To avoid 110 rounding errors, take care to specify a value that can 111 be exactly represented as a floating point number.]] 112[[y-offset auto][Sets the y-offset to the average value of f(x) 113 evaluated at the two endpoints of the range plus the midpoint 114 of the range. The calculated value is deliberately truncated 115 to /float/ precision (and should be stored as a /float/ 116 in your code). The approximation will 117 be generated for `f(x + X) + Y` where /X/ is the x-offset 118 and /Y/ is the y-offset. Defaults to zero.]] 119[[graph N][Prints N evaluations of f(x) at evenly spaced points over the 120 range being optimised. If unspecified then /N/ defaults 121 to 3. Use to check that f(x) is indeed smooth over the range 122 of interest.]] 123[[step N][Performs /N/ steps, or one step if /N/ is unspecified. 124 After each step prints: the peek error at the extrema of 125 the error function of the approximation, 126 the theoretical error term solved for on the last step, 127 and the maximum relative change in the location of the 128 Chebyshev control points. The approximation is converged on the 129 minimax solution when the two error terms are (approximately) 130 equal, and the change in the control points has decreased to 131 a suitably small value.]] 132[[test \[float|double|long\]][Tests the current approximation at float, 133 double, or long double precision. Useful to check for rounding 134 errors in evaluating the approximation at fixed precision. 135 Tests are conducted at the extrema of the error function of the 136 approximation, and at the zeros of the error function.]] 137[[test \[float|double|long\] N] [Tests the current approximation at float, 138 double, or long double precision. Useful to check for rounding 139 errors in evaluating the approximation at fixed precision. 140 Tests are conducted at N evenly spaced points over the range 141 of the approximation. If none of \[float|double|long\] are specified 142 then tests using NTL::RR, this can be used to obtain the error 143 function of the approximation.]] 144[[rescale a b][Takes the current Chebeshev control points, and rescales them 145 over a new interval \[a,b\]. Sometimes this can be used to obtain 146 starting control points for an approximation that can not otherwise be 147 converged.]] 148[[rotate][Moves one term from the numerator to the denominator, but keeps the 149 Chebyshev control points the same. Sometimes this can be used to obtain 150 starting control points for an approximation that can not otherwise be 151 converged.]] 152[[info][Prints out the current approximation: the location of the zeros of the 153 error function, the location of the Chebyshev control points, the 154 x and y offsets, and of course the coefficients of the polynomials.]] 155] 156 157[endsect] [/section:minimax Minimax Approximations and the Remez Algorithm] 158 159[/ 160 Copyright 2006 John Maddock and Paul A. Bristow. 161 Distributed under the Boost Software License, Version 1.0. 162 (See accompanying file LICENSE_1_0.txt or copy at 163 http://www.boost.org/LICENSE_1_0.txt). 164] 165