Re: [OctDev] complex error function

[Steven G. Johnson]

On 11/20/12 5:33 PM, Daniel J Sebald wrote
> How much of this is the Faddeeva function and how much of it is
> supporting code? One concern I would have, just looking at the code for
> a few minutes, is all the repetitive routines for square root, sinc,
> etc. I would think the preference should be to utilize such algorithms
> from Octave's core or utilize the library that Octave is using. That way
> things remain consistent.

I don't provide my own implementation of sqrt, so I'm not sure what you 
are talking about.

Regarding sinc, you are welcome to replace it with Octave's 
implementation, but that will be slower since I only use it in a case 
where sin(x) has already been computed and I can take advantage of that 
if x is not small.   Furthermore, we are talking about literally two 
lines of code, so I fail to see why including an additional two lines of 
code for the sake of performance is a dealbreaker for you:

// return sinc(x) = sin(x)/x, given both x and sin(x)
// [since we only use this in cases where sin(x) has already been computed]
static inline double sinc(double x, double sinx) {
  return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x;
}

The only other "standard" math function that I replace is a version of 
sinh(x) that is specialized for small |x| (since I only use it in that 
case).   Again, you could replace this with the libc sinh function at 
the cost of performance.  Again, you would save only three lines of code 
this way, so I don't see the benefit.

Everything else is specialized routines, needed for the various complex 
error functions and their special cases, that are not redundant anything 
in Octave as far as I can tell.

> Certainly error function and complimentary error function find general
> use, but Faddeeva function seems like something particular to wave
> theory or Fourier analysis. Whether that belongs in the special
> functions category or a wave-theory (or similar) package, I don't know.

Faddeeva is just a scaled error function of complex arguments.  Error 
functions of complex arguments have lots of applications as mentioned in 
my first email in this thread.


> This has an MIT license, correct? Is that compatible with GPL? From my
> brief browsing of the code, it does seem to be original and reference
> up-to-date papers on the algorithm and its limitations.

Yes, the MIT license is GPL compatible.  It is what FSF calls the "X11 
license":
        http://www.gnu.org/licenses/license-list.html#GPLCompatibleLicenses

And yes, it is original code by me.  (I originally used a public-domain 
function from SLATEC, but have replaced that with a faster algorithm 
written by me.)

> It seems to me the first question is whether to replace erf/erfc in the
> core source with something that handles complex inputs. I'm not sure
> Faddeeva function on its own warrants core code. But there might be a
> better way to put Faddeeva function in a package, script code being one
> possibility if that turns out to not be too slow.

The only difference between the Faddeeva function and erfc is the 
scaling.  I would strongly argue that you should provide either the 
Faddeeva function or the equivalent erfcx function (= Faddeeva function 
with an i factor).  The reason for the scaling is that the erfc function 
quickly underflows:  for x > 28, erfc(x) gives 0, which loses *all* 
significant digits.   In contrast, erfcx(x) = exp(x^2) erfc(x) never 
underflows for large x, so it allows you to look accurately at the tail 
of the error function for large arguments.

You already have an erfcx function in OF, although mine is faster.  Mine 
extends this to complex arguments.

--SGJ