Opened 12 years ago
Closed 7 years ago
#9263 closed defect (worksforme)
Some special functions don't work with arbitrary precision
Reported by: | zimmerma | Owned by: | jason, jkantor |
---|---|---|---|
Priority: | major | Milestone: | sage-duplicate/invalid/wontfix |
Component: | numerical | Keywords: | |
Cc: | eviatarbach | Merged in: | |
Authors: | Reviewers: | Jeroen Demeyer | |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description (last modified by )
---------------------------------------------------------------------- | Sage Version 4.4.2, Release Date: 2010-05-19 | | Type notebook() for the GUI, and license() for information. | ---------------------------------------------------------------------- sage: n(airy_ai(1),digits=100) 0.1352924163128813861423083153567858971655368804931640625000000000000000000000000000000000000000000000
Clearly the last digits are wrong. This is due to Maxima; currently we are not giving it bigfloats, and in any case some of the functions don't work with arbitrary precision.
The following functions should be changed to a different default backend (mpmath probably):
airy_ai
andairy_bi
(will be fixed with #12455)spherical_bessel_J
(it also has a SciPy? option, but it also truncates the precision)spherical_bessel_Y
(same asspherical_bessel_J
)spherical_harmonic
- the elliptic integrals
spherical_hankel1
andspherical_hankel2
hypergeometric_U
(PARI is implemented and works with arbitrary precision, but it's not the default)bessel_K
andbessel_Y
(now fixed with #4102)inverse_jacobi
Change History (28)
comment:1 Changed 12 years ago by
comment:2 Changed 9 years ago by
- Cc eviatarbach added
comment:3 Changed 9 years ago by
There appear to be two different problems related to numerical evaluation with Maxima. First, that some functions are locked to float precision. In Maxima:
(%i15) airy_ai(bfloat(%pi)),fpprec:20; (%o15) airy_ai(3.1415926535897932385b0)
I think it's returning unevaluated because airy_ai
doesn't know how to operate on bigfloat
s.
Other functions do know how to operate on bigfloat
s:
(%i18) bfloat(spherical_bessel_j(4, bfloat(%pi))),fpprec:200; (%o18) 6.471630031847746208103870635408583211756194941699504852294921875b-2
But, the interface is losing precision:
sage: spherical_bessel_J(4, pi.n(digits=1000)).n(digits=100) 0.06471630031847745712081376723290304653346538543701171875000000000000000000000000000000000000000000000
This is because Maxima truncates to float precision:
(%i20) 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825,fpprec:200; (%o20) 3.141592653589793
One way of avoiding this is converting to an exact rational and passing it to bfloat
, which may not be worth the overhead. Maybe something like the patch in #11643?
comment:4 Changed 9 years ago by
I believe the "correct" solution (for Sage) for now is to use #12455 and mpmath for (at least default) evaluation. I've cc:ed you on that ticket.
comment:5 Changed 9 years ago by
Oh, thank you. I guess there's not much use working on the Maxima problems when mpmath works well. I'll change this ticket to apply to all Maxima-evaluated functions.
comment:6 Changed 9 years ago by
- Component changed from basic arithmetic to numerical
- Description modified (diff)
- Owner changed from AlexGhitza to jason, jkantor
- Summary changed from airy_ai yields wrong results in arbitrary precision to Maxima-evaluated functions don't work with arbitrary precision
comment:7 Changed 9 years ago by
- Description modified (diff)
comment:8 Changed 9 years ago by
That is a good change. However, we don't have that many of those left, and so probably it will have to wait until we reliably have a keyword for evaluation algorithm.
comment:9 Changed 9 years ago by
Maybe there is a way to sense whether we are passing in something other than regular precision and then ask Maxima to use a bfloat.
comment:10 follow-up: ↓ 11 Changed 9 years ago by
- Description modified (diff)
- Summary changed from Maxima-evaluated functions don't work with arbitrary precision to Some special functions don't work with arbitrary precision
Do you mean #12289? What would be the problem with changing the backend before that's implemented though?
Actually, this ticket applies to many more functions than I thought initially. I added them to the description.
I also changed it again to apply to all special functions that don't work with arbitrary precision.
comment:11 in reply to: ↑ 10 Changed 9 years ago by
Do you mean #12289? What would be the problem with changing the backend before that's implemented though?
No problem at all, I just meant that it would make more sense to switch them to mpmath first, and then worry about getting Maxima to have the right precision after that ticket. Though I guess even spherical Bessel hasn't been implemented in mpmath yet...
comment:12 Changed 9 years ago by
- Milestone changed from sage-5.11 to sage-5.12
comment:13 Changed 9 years ago by
- Description modified (diff)
comment:14 Changed 9 years ago by
- Description modified (diff)
comment:15 Changed 9 years ago by
about hypergeometric_U
, the documentation (in Sage 5.11) says that Pari is the default, contrary to what is said in the description of this ticket. Which one is correct?
Paul
comment:16 Changed 9 years ago by
- Description modified (diff)
bessel_K
and bessel_Y
are fixed in Sage 5.11 thanks to #4102, thus I update the description:
sage: n(bessel_K(1,2), prec=100) 0.13986588181652242728459880704 sage: n(bessel_Y(1,2), prec=100) -0.10703243154093754688837077228
Paul
comment:17 Changed 9 years ago by
- Milestone changed from sage-6.1 to sage-6.2
comment:18 Changed 8 years ago by
- Milestone changed from sage-6.2 to sage-6.3
comment:19 Changed 8 years ago by
- Milestone changed from sage-6.3 to sage-6.4
comment:20 Changed 7 years ago by
- Milestone changed from sage-6.4 to sage-duplicate/invalid/wontfix
- Reviewers set to Jeroen Demeyer
- Status changed from new to needs_review
Works for me.
comment:21 Changed 7 years ago by
- Status changed from needs_review to positive_review
comment:22 Changed 7 years ago by
hypergeometric_U
is strange, since the answer is always 53 bits:
sage: hypergeometric_U(1,2,3) 0.333333333333333 sage: R = RealField(1000) sage: hypergeometric_U(R(1), R(2), R(3)) 0.333333333333333 sage: hypergeometric_U(R(1), R(2), R(3)).n(100) ... TypeError: cannot approximate to a precision of 100 bits, use at most 53 bits
The others work properly.
comment:23 Changed 7 years ago by
And the issue with hypergeometric_U
is #14896.
comment:24 Changed 7 years ago by
with Sage 6.0 I get:
sage: R = RealField(1000) sage: hypergeometric_U(R(1), R(2), R(3)).n(100) 0.33333333333333331482961625625 sage: hypergeometric_U(1,2,3).n(100) 0.33333333333333331482961625625
Maybe a regression?
comment:25 Changed 7 years ago by
is the issue with airy_ai
fixed? I still get with Sage 6.0:
sage: n(airy_ai(1),digits=75) 0.135292416312881413897883930985699407756328582763671875000000000000000000000
comment:26 Changed 7 years ago by
Well, Sage 6.0 is ancient.
With 6.7.beta4:
sage: n(airy_ai(1),digits=75) 0.135292416312881415524147423515466306174944142988330706009102054757633534802
and
sage: hypergeometric_U(1,2,3).n(100) --------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-4-7c7c457a86be> in <module>() ----> 1 hypergeometric_U(Integer(1),Integer(2),Integer(3)).n(Integer(100)) /usr/local/src/sage-config/src/sage/structure/element.pyx in sage.structure.element.Element.numerical_approx (build/cythonized/sage/structure/element.c:7437)() 744 """ 745 from sage.misc.functional import numerical_approx --> 746 return numerical_approx(self, prec=prec, digits=digits, 747 algorithm=algorithm) 748 n = numerical_approx /usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/misc/functional.pyc in numerical_approx(x, prec, digits, algorithm) 1329 1330 if prec > inprec: -> 1331 raise TypeError("cannot approximate to a precision of %s bits, use at most %s bits" % (prec, inprec)) 1332 1333 # The issue is not precision, try conversion instead TypeError: cannot approximate to a precision of 100 bits, use at most 53 bits
comment:27 Changed 7 years ago by
Note that Maxima functions we didn't usually have a good precision interface with, if that is what is going on with the hgu.
comment:28 Changed 7 years ago by
- Resolution set to worksforme
- Status changed from positive_review to closed
one solution would be to use mpmath: