Grid/lib/algorithms/approx

-----------------------------------------------------------------------------------

PAB. Took Mike Clark's AlgRemez from GitHub and (modified a little) include.
This is open source and license and readme and comments are preserved consistent
with the license. Mike, thankyou!
-----------------------------------------------------------------------------------
-----------------------------------------------------------------------------------
AlgRemez

The archive downloadable here contains an implementation of the Remez
algorithm which calculates optimal rational (and polynomial)
approximations to the nth root over a given spectral range.  The Remez
algorithm, although in principle is extremely straightforward to
program, is quite difficult to get completely correct, e.g., the Maple
implementation of the algorithm does not always converge to the
correct answer.

To use this algorithm you need to install GMP, the GNU Multiple
Precision Library, and when configuring the install, you must include
the --enable-mpfr option (see the GMP manual for more details).  You
also have to edit the Makefile for AlgRemez appropriately for your
system, namely to point corrrectly to the location of the GMP library.

The simple main program included with this archive invokes the
AlgRemez class to calculate an approximation given by command line
arguments.  It is invoked by the following

./test y z n d lambda_low lambda_high precision,

where the function to be approximated is f(x) = x^(y/z), with degree
(n,d) over the spectral range [lambda_low, lambda_high], using
precision digits of precision in the arithmetic.  So an example would
be

./test 1 2 5 5 0.0004 64 40

which corresponds to constructing a rational approximation to the
square root function, with degree (5,5) over the range [0.0004,64]
with 40 digits of precision used for the arithmetic.  The parameters y
and z must be positive, the approximation to f(x) = x^(-y/z) is simply
the inverse of the approximation to f(x) = x^(y/z).  After the
approximation has been constructed, the roots and poles of the
rational function are found, and then the partial fraction expansion
of both the rational function and it's inverse are found, the results
of which are output to a file called "approx.dat".  In addition, the
error function of the approximation is output to "error.dat", where it
can be checked that the resultant approximation satisfies Chebychev's
criterion, namely all error maxima are equal in magnitude, and
adjacent maxima are oppostie in sign.  There are some caveats here
however, the optimal polynomial approximation has complex roots, and
the root finding implemented here cannot (yet) handle complex roots.
In addition, the partial fraction expansion of rational approximations
is only found for the case n = d, i.e., the degree of numerator
polynomial equals that of the denominator polynomial.  The convention
for the partial fraction expansion is that polar shifts are always
written added to x, not subtracted.

To do list

1.  Include an exponential dampening factor in the function to be
approximated.  This may sound trivial to implement, but for some
parameters, the algorithm seems to breakdown.  Also, the roots in the
rational approximation sometimes become complex, which currently
breaks the stupidly simple root finding code.

2. Make the algorithm faster - it's too slow when running on qcdoc.

3. Add complex root finding.

4. Add more options for error minimisation - currently the code
minimises the relative error, should add options for absolute error,
and other norms.

There will be a forthcoming publication concerning the results
generated by this software, but in the meantime, if you use this
software, please cite it as
"M.A. Clark and A.D. Kennedy, https://github.com/mikeaclark/AlgRemez, 2005".

If you have any problems using the software, then please email scientist.mike@gmail.com.