libs/math/doc/internals/minimax.qbk

   1 [section:minimax Minimax Approximations and the Remez Algorithm]
   2
   3 The directory `libs/math/minimax` contains an interactive command-line driven
   4 program for the generation of minimax approximations using the Remez
   5 algorithm.  Both polynomial and rational approximations are supported,
   6 although the latter are tricky to converge: it is not uncommon for
   7 convergence of rational forms to fail.  No such limitations are present
   8 for polynomial approximations which should always converge smoothly.
   9
  10 It's worth stressing that developing rational approximations to functions
  11 is often not an easy task, and one to which many books have been devoted.
  12 To use this tool, you will need to have a reasonable grasp of what the Remez
  13 algorithm is, and the general form of the approximation you want to achieve.
  14
  15 Unless you already familar 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
  18 The 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
  25 Therefore to use this tool, you must modify f.cpp to return the function to
  26 approximate.  The tools supports multiple function approximations within
  27 the same compiled program: each as a separate variant:
  28
  29    NTL::RR f(const NTL::RR& x, int variant);
  30
  31 Returns the value of the function /variant/ at point /x/.  So if you
  32 wish you can just add the function to approximate as a new variant
  33 after the existing examples.
  34
  35 In addition to those two files, the program needs to be linked to
  36 a [link math_toolkit.high_precision.use_ntl patched NTL library to compile].
  37
  38 Note that the function /f/ must return the rational part of the
  39 approximation: 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
  44 where /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
  46 absolute error compared to |Y|.
  47
  48 In this case you would define /f/ to return [role serif-italic f(x)/g(x)] and then set the
  49 y-offset of the approximation to /Y/ (see command line options below).
  50
  51 Many other forms are possible, but in all cases the objective is to
  52 split /f(x)/ into a dominant part that you can evaluate easily using
  53 standard math functions, and a smooth and slowly changing rational approximation
  54 part.  Refer to your favourite textbook for more examples.
  55
  56 Command 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 ]