#1173: implement numerical evaluation of erf at complex arguments
----------------------------------------------------------------------+-----
Reporter: was |
Owner: was
Type: enhancement |
Status: needs_work
Priority: major |
Milestone: sage-4.7.2
Component: calculus |
Keywords:
Work_issues: add erf(sqrt(2)) and erf(45000).n(), formatting, speed |
Upstream: N/A
Reviewer: Burcin Erocal |
Author: D. S. McNeil
Merged: |
Dependencies: #11513
----------------------------------------------------------------------+-----
Comment(by dsm):
Now I remember why I found this so frustrating. I rewrote the patch
following -- plagiarizing is more like it! -- the pattern used in #11143,
but the speed issue hasn't gone away.
{{{
# 4.7.1, OS X 10.6.8
sage: timeit('erf(0.1)',number=1000)
1000 loops, best of 3: 82.9 µs per loop
sage: timeit('erf(10.0)',number=1000)
1000 loops, best of 3: 72.8 µs per loop
sage: timeit('erf(100.0)',number=1000)
1000 loops, best of 3: 73.5 µs per loop
# 4.7.2 before patch
sage: timeit('erf(0.1)',number=1000)
1000 loops, best of 3: 69.4 µs per loop
sage: timeit('erf(10.0)',number=1000)
1000 loops, best of 3: 62.6 µs per loop
sage: timeit('erf(100.0)',number=1000)
1000 loops, best of 3: 62 µs per loop
# 4.7.2 after patch
sage: timeit('erf(0.1)',number=1000)
1000 loops, best of 3: 137 µs per loop
sage: timeit('erf(10.0)',number=1000)
1000 loops, best of 3: 116 µs per loop
sage: timeit('erf(100.0)',number=1000)
1000 loops, best of 3: 116 µs per loop
sage: import mpmath
sage: timeit('mpmath.erf(0.1)')
625 loops, best of 3: 95 µs per loop
sage: timeit('mpmath.erf(10.0)')
625 loops, best of 3: 75.4 µs per loop
sage: timeit('mpmath.erf(100.0)')
625 loops, best of 3: 76.2 µs per loop
}}}
That is, there's about a factor of two penalty in speed for the standard
case, but it's not because the underlying mpmath code is slow:
{{{
sage: timeit('mpmath.erf(0.1r)')
625 loops, best of 3: 27.7 µs per loop
sage: timeit('mpmath.erf(10.0r)')
625 loops, best of 3: 9.85 µs per loop
sage: timeit('mpmath.erf(100.0r)')
625 loops, best of 3: 9.95 µs per loop
}}}
In fact, mpmath isn't that much slower than MPFR:
{{{
sage: z = RR(2); timeit('z.erf()')
625 loops, best of 3: 20 µs per loop
sage: z = 2.0r; timeit('mpmath.erf(z)')
625 loops, best of 3: 57.9 µs per loop
sage: timeit('erf(2.0)') # new code
625 loops, best of 3: 181 µs per loop
}}}
So not much of the total time is spent actually doing any calculations:
it's all overhead. :-( This affects #11143 as well:
{{{
sage: timeit('exp_integral_e1(2.0)')
625 loops, best of 3: 165 µs per loop
sage: z = exp_integral_e1(2); timeit('z.n()')
625 loops, best of 3: 143 µs per loop
sage: timeit('exponential_integral_1(2.0)')
625 loops, best of 3: 44.7 µs per loop
sage: timeit('mpmath.e1(2.0)')
625 loops, best of 3: 123 µs per loop
sage: timeit('mpmath.e1(2.0r)')
625 loops, best of 3: 51 µs per loop
}}}
To wrap up:
(1) Both this patch and #11143 suffer a significant slowdown relative to
PARI, and have major overheads relative to calling mpmath. Some of that's
unavoidable given the symbolic path, but ISTM we should be able to do
better. Ideally there would be a reasonably efficient general special
function implementation pattern, along the lines of what Benjamin used,
that intermittent developers like me could be pointed to as a reference.
(2) In the case of erf and erfc, mpfr offers a very fast evaluation for
real numbers, fast enough that it might be worth using as the default in
those cases (although Python-level branching is slow in my experience,
maybe slow enough to wash away the benefits). Once we settle on an
approach for erf I'll do the same for erfc.
(3) Should I move this out of other.py to special.py, where the
complementary error_fcn function lives now? It would seem a more natural
location for it. We also have some unfortunate naming (erf and error_fcn)
it might be worth addressing.
I don't have numbers for #11948 -- too many dependencies -- but it's
probably considerably faster than this approach. I figure it's probably
worth getting the #11143 -style mpmath wrapper to be faster, though, even
if we went #11948 instead for erf/erfc.
[There are a few doctest failures:
"devel/sage/doc/en/bordeaux_2008/l_series.rst", which I think is
unrelated; a timeout in devel/sage/sage/modules/free_module.py, again
probably unrelated.]
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/1173#comment:26>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
