#11948 closed defect (fixed)
Fix numeric evaluation of error function
Reported by: | jdemeyer | Owned by: | jdemeyer |
---|---|---|---|
Priority: | major | Milestone: | sage-4.8 |
Component: | numerical | Keywords: | erf erfc |
Cc: | Merged in: | sage-4.8.alpha6 | |
Authors: | Jeroen Demeyer | Reviewers: | Karl-Dieter Crisman |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #11130, #11321, #11854, #11891, #11890, #11836, #11952 | Stopgaps: |
Description
Because the argument is converted to float
first, the error function erf()
cannot numerically be evaluated for complex arguments:
sage: erf(pi - 1/2*I).n() ERROR: An unexpected error occurred while tokenizing input The following traceback may be corrupted or invalid The error message is: ('EOF in multi-line statement', (1108, 0)) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /usr/local/src/sage-4.7.2.alpha2/<ipython console> in <module>() /usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (sage/symbolic/expression.cpp:17793)() /usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (sage/symbolic/expression.cpp:5004)() /usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/functions/other.pyc in _evalf_(self, x, parent) 79 if prec > 53: 80 raise NotImplementedError, "erf not implemented for precision higher than 53" ---> 81 return parent(1 - pari(float(x)).erfc()) 82 83 def _derivative_(self, x, diff_param=None): /usr/local/src/sage-4.7.2.alpha2/local/lib/python2.6/site-packages/sage/rings/complex_number.so in sage.rings.complex_number.ComplexNumber.__float__ (sage/rings/complex_number.c:7233)() TypeError: Unable to convert 3.14159265358979 - 0.500000000000000*I to float; use abs() or real_part() as desired
Attachments (2)
Change History (20)
comment:1 Changed 8 years ago by
comment:2 Changed 8 years ago by
- Status changed from new to needs_review
comment:3 Changed 8 years ago by
comment:4 Changed 8 years ago by
- Status changed from needs_review to needs_work
comment:5 Changed 8 years ago by
- Dependencies changed from #11130 to #11130, #11321, #11854, #11891, #11890, #11836, #11904, #11952
- Status changed from needs_work to needs_review
comment:6 follow-up: ↓ 7 Changed 8 years ago by
Thanks for that info; I didn't realize it was the new Pari that made this possible.
I'd say that if someone reviews this, it should go in. #1173 is probably not going to get finished immediately, and we can always switch to mpmath if necessary or convenient or speedy at that time.
comment:7 in reply to: ↑ 6 Changed 8 years ago by
- Owner changed from jason, jkantor to jdemeyer
comment:8 Changed 8 years ago by
- Dependencies changed from #11130, #11321, #11854, #11891, #11890, #11836, #11904, #11952 to #11130, #11321, #11854, #11891, #11890, #11836, #11952
Changed 8 years ago by
comment:9 follow-up: ↓ 10 Changed 8 years ago by
- Reviewers set to Karl-Dieter Crisman
- Status changed from needs_review to positive_review
Looks good. I've added a tiny bit of doc, because erf?
gives something very short otherwise now, and fixed some grammar in the place where Jeroen did something I can't figure out with whitespace. Docs still look fine there, so all seems well.
Positive review, unless the author finds something he doesn't like about the review patch.
comment:10 in reply to: ↑ 9 Changed 8 years ago by
Replying to kcrisman:
Jeroen did something I can't figure out with whitespace.
Probably changing TABS to spaces.
comment:11 Changed 8 years ago by
- Milestone changed from sage-pending to sage-5.0
comment:12 Changed 8 years ago by
- Milestone changed from sage-5.0 to sage-4.8
comment:13 Changed 8 years ago by
- Merged in set to sage-4.8.alpha6
- Resolution set to fixed
- Status changed from positive_review to closed
comment:14 Changed 8 years ago by
Unfortunately, this causes a nasty problem.
This example in the documentation is fine - comparing mpmath:
sage: mpmath.erf(pi-1/2*i) mpc(real='1.0000111669099367', imag='1.6332655417638451e-6')
But anything along the imaginary axis seems to be off by exactly 1:
sage: for z in [3,33,333,3333,33333]: ....: mpmath.erf(i*z); erf(n(z)*i) ....: mpc(real='0.0', imag='1629.9946226015657') 1.00000000000000 + 1629.86732385786*I mpc(real='0.0', imag='1.5128697751040891e+471') 1.00000000000000 + 1.51286977510409e471*I mpc(real='0.0', imag='5.1260939089106243e+48155') 1.00000000000000 + 5.12609390891062e48155*I mpc(real='0.0', imag='2.6385510598470926e+4824525') 1.00000000000000 + 2.63855105984709e4824525*I
But other values seem ok.
sage: for z in [3,33,333,3333]: mpmath.erf(1+i*z); erf(1.+n(z)*i) ....: mpc(real='-330.81538696857206', imag='443.38888183939281') -330.815386892947 + 443.388881909712*I mpc(real='2.0957487368415288e+468', imag='-5.5629367605580166e+470') 2.09574873684153e468 - 5.56293676055802e470*I mpc(real='-3.8930178706420656e+48153', imag='1.8853741770265906e+48155') -3.89301787064206e48153 + 1.88537417702659e48155*I mpc(real='-4.3084905090066053e+4824524', imag='8.6980843586535772e+4824524') -4.30849050900660e4824524 + 8.69808435865358e4824524*I sage: for z in [3,33,333,3333,33333]: mpmath.erf(-1-i*z); erf(-1.-n(z)*i) ....: mpc(real='330.81538696857206', imag='-443.38888183939281') 330.815386968572 - 443.388881839393*I mpc(real='-2.0957487368415288e+468', imag='5.5629367605580166e+470') -2.09574873684153e468 + 5.56293676055802e470*I mpc(real='3.8930178706420656e+48153', imag='-1.8853741770265906e+48155') 3.89301787064207e48153 - 1.88537417702659e48155*I mpc(real='4.3084905090066053e+4824524', imag='-8.6980843586535772e+4824524') 4.30849050900661e4824524 - 8.69808435865358e4824524*I
comment:15 Changed 8 years ago by
More precisely:
sage: pari(3*i).erfc() -1.76710569338983 E-16 - 1629.86732385786*I sage: mpmath.erfc(3*i) mpc(real='1.0', imag='-1629.9946226015657') sage: 1-pari(3*i).erfc() 1.00000000000000 + 1629.86732385786*I sage: mpmath.erf(3*i) mpc(real='0.0', imag='1629.9946226015657')
comment:16 Changed 8 years ago by
Now, in 5.1.beta6's gp I get
$ Downloads/sage-5.1.beta6/sage -gp GP/PARI CALCULATOR Version 2.5.1 (development git-5c5e253) i386 running darwin (x86-64/GMP-5.0.2 kernel) 64-bit version compiled: Jun 25 2012, gcc-4.6.3 (GCC) (readline v6.2 enabled, extended help enabled) Copyright (C) 2000-2011 The PARI Group PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER. Type ? for help, \q to quit. Type ?12 for how to get moral (and possibly technical) support. parisize = 8000000, primelimit = 500509 ? erfc(3*I) %1 = 0.E-35 - 1629.5516567497094550267455532288372861*I
so I see that this is still in "our" version of Pari. Is this fixed upstream? Is it possible that this is not actually a bug in Pari, but in our use of it? Wolfram Alpha gives the mpmath answer, for what it's worth.
comment:17 Changed 8 years ago by
I should also point out that we even doctest that it is wrong - due to my reviewer patch on this ticket :( so this is certainly a mea culpa among others.
comment:18 Changed 8 years ago by
And apparently even weirder, courtesy of Jeff Denny.
sage: erf(i*1.42) 1.00000000000000 + 4.03986343036907*I sage: import mpmath sage: mpmath.erf(i*1.42) mpc(real='0.0', imag='3.8217653554366318')
I'm opening a ticket for this, even if it gets fixed by one of the logical others (#1173 or #13050) because it is just mathematically wrong. This is #13193; continue discussion there.
#1173 is closely related, and #8983 is somewhat related. I haven't tried whether #1173 solves this issue or not.