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:

Status badges

Description (last modified by zimmerma)

----------------------------------------------------------------------
| 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 and airy_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 as spherical_bessel_J)
  • spherical_harmonic
  • the elliptic integrals
  • spherical_hankel1 and spherical_hankel2
  • hypergeometric_U (PARI is implemented and works with arbitrary precision, but it's not the default)
  • bessel_K and bessel_Y (now fixed with #4102)
  • inverse_jacobi

Change History (28)

comment:1 Changed 12 years ago by ylchapuy

one solution would be to use mpmath:

sage: import mpmath
sage: mpmath.mp.dps = 100
sage: mpmath.airyai(1)
mpf('0.1352924163128814155241474235154663061749441429883307060091020547576335348022657236634871099087486832138')

comment:2 Changed 9 years ago by eviatarbach

  • Cc eviatarbach added

comment:3 Changed 9 years ago by eviatarbach

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 bigfloats.

Other functions do know how to operate on bigfloats:

(%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?

Last edited 9 years ago by eviatarbach (previous) (diff)

comment:4 Changed 9 years ago by kcrisman

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 eviatarbach

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 eviatarbach

  • 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 eviatarbach

  • Description modified (diff)

comment:8 Changed 9 years ago by kcrisman

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 kcrisman

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: Changed 9 years ago by eviatarbach

  • 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 kcrisman

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 jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:13 Changed 9 years ago by zimmerma

  • Description modified (diff)

comment:14 Changed 9 years ago by zimmerma

  • Description modified (diff)

comment:15 Changed 9 years ago by zimmerma

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 zimmerma

  • 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 vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:18 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:19 Changed 8 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:20 Changed 7 years ago by jdemeyer

  • 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 jdemeyer

  • Status changed from needs_review to positive_review

comment:22 Changed 7 years ago by jdemeyer

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 jdemeyer

And the issue with hypergeometric_U is #14896.

comment:24 Changed 7 years ago by zimmerma

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 zimmerma

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 jdemeyer

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 kcrisman

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 vbraun

  • Resolution set to worksforme
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.