In Python3.2, calling math.erfc with a value in [-27.2, -30) raises
an OverflowError: math range error. This is inconsistent with the
erfc function from scipy (scipy.special.erfc) as well as with the C99
function by the same name, both of which return 2. I suspect that
this is the result of the cutoff for the use of the continuing fraction
approximation of erfc beginning when abs(x) > 30, but I'm not sure.

Thanks for the report.  What platform are you on?  I'm not seeing this behaviour on OS X:

Python 3.2a0 (py3k:81935M, Jun 12 2010, 10:01:38) 
[GCC 4.2.1 (Apple Inc. build 5659)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from math import erfc
>>> erfc(-28.0)
2.0

It should be easy to fix, since erfc is already 2.0 (to within machine accuracy) in that argument range, but I'd like to understand where the problem is coming from first.

Looking at the source, I think I know where this is coming from:  one of the terms in the expression for erfc(x) is exp(-x*x).  For abs(x) >= 27.284 or so, this underflows to zero.

So the likely cause is that whatever platform you're on is setting errno to ERANGE on underflow to 0;  that errno value later gets interpreted as representing an overflow.

And for x >= 30.0 a constant expression is used instead, so there's no OverflowError any more.

Okay:  solved in principle!

Geremy, can you verify that the attached patch fixes the problem?  If so, I'll add tests and commit.

On Sun, Jun 13, 2010 at 5:14 AM, Mark Dickinson <report@bugs.python.org> wrote:
>
> Mark Dickinson <dickinsm@gmail.com> added the comment:
>
> Thanks for the report.  What platform are you on?  I'm not seeing this behaviour on OS X:

I'm on Ubuntu 10.04 64-bit.

Geremy Condra

> Python 3.2a0 (py3k:81935M, Jun 12 2010, 10:01:38)
> [GCC 4.2.1 (Apple Inc. build 5659)] on darwin
> Type "help", "copyright", "credits" or "license" for more information.
>>>> from math import erfc
>>>> erfc(-28.0)
> 2.0
>
> It should be easy to fix, since erfc is already 2.0 (to within machine accuracy) in that argument range, but I'd like to understand where the problem is coming from first.
>
> ----------
>
> _______________________________________
> Python tracker <report@bugs.python.org>
> <http://bugs.python.org/issue8986>
> _______________________________________
>

On Sun, Jun 13, 2010 at 5:47 AM, Mark Dickinson <report@bugs.python.org> wrote:
>
> Mark Dickinson <dickinsm@gmail.com> added the comment:
>
> Geremy, can you verify that the attached patch fixes the problem?  If so, I'll add tests and commit.

I've only tested it from 2**16 + 1 to -2**16 at unit and half-unit
intervals, but it doesn't blow up
here. Thanks- and I appreciate the quick response.

Geremy Condra

> ----------
> keywords: +patch
> Added file: http://bugs.python.org/file17654/issue8986.patch
>
> _______________________________________
> Python tracker <report@bugs.python.org>
> <http://bugs.python.org/issue8986>
> _______________________________________
>

Fixed in r81967 (trunk) and r81968 (py3k).

I had to weaken the tests for erfc:  its accuracy for largish arguments (25.0 or so) is not ideal---I was seeing errors of 100 ulps and more.  However, I think this level of error is acceptable in practice, since for arguments this size the difference between erfc(x) and erfc(next_float_up(x)) is already several hundred ulps.

If anyone wants to look at improving the accuracy for large arguments, it's the calculation of exp(-x*x) that's the source of the error.  The way to fix it, if anyone cared, would be to do the multiplication x*x in double-double precision, giving x*x = y + z for doubles y and z, and then use exp(-y) * exp(-z) in place of exp(-x*x).

Thanks for catching this, Geremy!