Re: complexity of numerical software

This is a repost of a message which appeared at my news server,
but hasn't showed up in Google groups for several hours. My
apologies if it's actually duplicated.

Jean-Claude Arbaut wrote:
> Jaap Spies wrote:
>
>> Jean-Claude Arbaut wrote:
>>
>>> Jaap Spies wrote:
>>>
>>>> David N. Williams wrote in a reaction to carlos@xxxxxxxxxxxx:
>>>>
>>>>> I'm really puzzled by your classification of numerical software.
>>>>> I would have put it in the truly challenging category!?
>>>>>
>>>>> -- David

>>>> Numerical software? What kind of numerical software you are talking
>>>> about?

>>> Maybe it would be intersting to know your definition of "numerical
>>> software"...
>>>
>>
>> Yes. We have the same question!
>
>
> Okay. I thought we were talking about code like SLATEC, NSWC, CMLIB,
> etc. They have (almost) nothing to do with algorithms to find subsets,
> which I would classify with discrete or combinatorics algorithms.

SLATEC, etc., are indeed examples of the kind of thing I had in
mind.

> It's easy to write poor numerical algorithms (with my definition),
> the usual example is "Numerical Recipes", but I find much more
> difficult to write *good* (read /accrurate and fast/) ones. Just
> try to write a correctly rounded sine (or even a square root !).

Jaap Spies also wrote:
> [...]
>
> Testing and debugging of numeric algorithms is easy. Just check the
> results.

Actually, I would include testing numerical analysis software in
the truly challenging category. A classic example is William
Cody's celefunt suite of accuracy tests for complex elementary
functions. It makes clever use of mathematical identities, but
while useful it is still limited, as Cody points out in his ACM
article on Algorithm 714. Surely it is nontrivial in general to
analyze where inaccuracies for a numerical program are likely to
occur, to even begin testing.

Jean-Claude mentioned the square root. Here's a quote (up to
math symbols) from one of those truly wonderful old IBM manuals
for VS FORTRAN:

Effect of Argument Error

epsilon ~ delta/2

Accuracy

The accuracy of the algorithm is perfect; for all (15*2^27)
nonnegative normalized short floating-point numbers x, SQRT
returns the correctly rounded value...

There's clearly nontrivial error analysis underlying IBM's
implementation of SQRT, which you have to worry about if you
don't have access to a trusted implementation like theirs. Then
the known relationship between the input error, delta, and the
output error, epsilon, maybe gives you some hope of
understanding the results of accuracy tests on code that uses
SQRT.

And the subtleties of floating-point arithmetic feed in as well.

-- David

.

Relevant Pages

  • Re: Fixed Point Square Root Improvements?
    ... non-LUT solution which is about 2-2.5 times slower than sqrt() on my Athlon ... I compared Newton-Rhapson to these algorithms also, ... bsr eax, eax ...
    (comp.lang.asm.x86)
  • L1 - Median (a multivariate median)
    ... There are several papers that give algorithms that are described as ... implementations are not robust to finite precision arithmetic. ... where epsilon is a small number. ...
    (sci.image.processing)
  • Re: Cascading different algorithms?
    ... >Is the idea that you would have SQRT states and they would be running ... that quantum computers may be easier to build than predicted. ... there is doubt developing new QC resistant algorithms. ...
    (sci.crypt)