Fix coercion bugs in symbolic functions

[jdemeyer]

This uses coercion correctly:

sage: bessel_Y._eval_(RealField(300)(1), 1.0)
-0.781212821300289

However, it seems that __call__() coerces this result back to the first parent, giving false precision:

sage: bessel_Y(RealField(300)(1), 1.0)
-0.781212821300288684511770043172873556613922119140625000000000000000000000000000000000000000

Same issue with functions which are evaluated using Maxima, which does not support arbitrary precision:

sage: R=RealField(300); elliptic_eu(R(1/2), R(1/8))
0.495073732023201484864216581627260893583297729492187500000000000000000000000000000000000000

The gamma_inc() function also mishandles parents:

sage: gamma_inc(float(0), float(1))
AttributeError: type object 'float' has no attribute 'precision'

---

Apart from this, this branch also removes lots of boilerplate from _eval_ like

if not isinstance(x, Expression) and not isinstance(y, Expression) and \
        (is_inexact(x) or is_inexact(y)):
    x, y = coercion_model.canonical_coercion(x, y)
    return self._evalf_(x, y, s_parent(x))

by wrapping _eval_ inside the new method _evalf_or_eval_ which automatically does this boilerplate.

---

Possible follow-ups: #10133, #14766, #16587, #17122, #15200

[kcrisman]

Just to clarify, you mean that the precision of the less precise input (1.0) is used, but then for some reason this is sent back to the precision of the first one?

My guess is that ​this line is to blame.

        # we want to convert the result to the original parent if the input
        # is not exact, so we store the parent here
        org_parent = parent_c(args[0])

I don't have time right now to check this out, but I would not be at all surprised. Or it's the super call just above that, I don't myself 100% understand how the inheritance works, and I keep forgetting when eval versus call are called in the functions. You could try other functions with multiple arguments and see if the same problem happens, or try to reverse the precisions and see what happens to test this hypothesis.

[jdemeyer]

Replying to kcrisman:

Just to clarify, you mean that the precision of the less precise input (1.0) is used, but then for some reason this is sent back to the precision of the first one?

Yes, the less-precise result is converted to the more-precise parent.

[jdemeyer]

Replying to kcrisman:

        # we want to convert the result to the original parent if the input
        # is not exact, so we store the parent here
        org_parent = parent_c(args[0])

Especially the args[0] is very suspicious indeed.

[jdemeyer]

The following patch fixes the problem

src/sage/symbolic/function.pyx

@@ -914 +914 @@
         res = super(BuiltinFunction, self).__call__(
                         *args, coerce=coerce, hold=hold)
 
-        # we want to convert the result to the original parent if the input
-        # is not exact, so we store the parent here
-        org_parent = parent_c(args[0])
-
-        # convert the result back to the original parent previously stored
-        # otherwise we end up with
-        #     sage: arctan(RR(1))
-        #     1/4*pi
-        # which is surprising, to say the least...
-        if org_parent is not SR and \
-                (org_parent is float or org_parent is complex or \
-                (PY_TYPE_CHECK(org_parent, Parent) and \
-                    not org_parent.is_exact())):
-            try:
-                return org_parent(res)
-            except (TypeError, ValueError):
-                pass
-
-            # conversion to the original parent failed
-            # we try if it works with the corresponding complex domain
-            if org_parent is float or org_parent is complex:
-                try:
-                    from sage.rings.complex_double import CDF
-                    return complex(CDF(res))
-                except (TypeError, ValueError):
-                    pass
-            elif hasattr(org_parent, 'complex_field'):
-                try:
-                    return org_parent.complex_field()(res)
-                except (TypeError, ValueError):
-                    pass
-
         return res
 
     cdef _is_registered(self):

but with quite a bit of doctest failures.

[jdemeyer]

I think a correct result with the wrong parent is better than a wrong result with the correct parent, so I'm inclined to really remove that block of code and fix individual functions instead.

[kcrisman]

I'm not surprised - the problem is that this kind of code was added at various times to handle certain cases where the output was not the same type as the input (e.g., complex output for real/float input). For instance, does the example mentioned in the code

arctan(RR(1)) 

work properly with this patch? Again, I'm sorry I don't have time to do more than read code for five minutes at this point :(

Edit: I see your most recent comment - yes, that would certainly be helpful, though I have a feeling a LOT of functions would have to be fixed... because some probably implicitly rely on this block but aren't thoroughly tested for unusual input/output.

Edit by jdemeyer: Sorry, edit instead of reply

[kcrisman]

Perhaps one could take the coercion mutual parents of all args instead for a minimal change? (Would that work?) I think in principle it's better to handle this all at once rather than having things in each individual function - if possible, of course.

[jdemeyer]

Replying to kcrisman:

arctan(RR(1)) 

work properly with this patch?

This works as it should, many doctest failures involve Python types:

sage: arctan(float(1))
1/4*pi

(I guess this goes through Pynac)

[jdemeyer]

Replying to kcrisman:

Perhaps one could take the coercion mutual parents of all args instead for a minimal change?

That still wouldn't fix the case where the function is computed to less precision than the inputs.

I think the following would fix most issues:

• before calling the function, convert all Python types to the corresponding Sage types (float -> RDF and so on)
• after calling the function, convert back if needed

[vbraun]

Not all arguments need to be floating point types, there could be integer or string parameters. So you can't indiscriminately coerce.

IMHO we should just add a clause that if output is float and of lower precision then use that lower precision. In an ideal world we would be able to evaluate everything to arbitrary precision, of course.

[jdemeyer]

New commits:

​6d10729

Simplify numerical evaluation of BuiltinFunctions

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0001941	Improve accuracy of polytopes.regular_polygon()
​24020fe	Partition().to_exp() should return Sage Integers
​b6e1ed4	Merge remote-tracking branches 'origin/u/jdemeyer/ticket/17131' and 'origin/u/jdemeyer/ticket/17133' into ticket/17130

[kcrisman]

Wow, that first commit is a nice patch bomb. I hesitate to make so much change to how symbolic functions are evaluated without taking a pretty close look, my apologies for not trying to do that immediately. Also, I think there was an actual reason for calling it inexact and not numerical, though I don't remember offhand.

[jdemeyer]

Replying to kcrisman:

Wow, that first commit is a nice patch bomb.

Would you be more willing to review the patch if it would be split up? I could try to do that but only if I have confirmation that it will increase the chances of this ticket being reviewed.

[vbraun]

I don't think it is overly long, nor is there a good way of splitting it up. Patch looks good to me.

Its somewhat confusing that is_numerical is completely different from GiNaC::numeric. How about is_approximate, is_limited_precision, or is_floating_point (I'd be happy to call padics "floating point" in this case).

Is it possible to replace is_inexact with is_numerical everywhere? The difference is

sage: is_inexact(pi)
True
sage: sin._is_numerical(pi)
False

I'm not super happy with pi being "inexact", though perhaps there is a technical reason for why we need it.

[kcrisman]

• Deprecation, folks.
• No, it isn't long, but I at least need to think about the mechanism by which things are evaluated and want to make sure we don't miss any odd cases. _eval_ has been fairly standard for quite some time.
• I agree that the naming is problematic no matter what name you choose. I think the rings folks had some good reasons for their choices:

      cpdef bint is_exact(self) except -2:
          """
          Test whether the ring is exact.
  
          .. NOTE::
  
              This defaults to true, so even if it does return ``True``
              you have no guarantee (unless the ring has properly
              overloaded this).
  
          OUTPUT:
  
          Return True if elements of this ring are represented exactly, i.e.,
          there is no precision loss when doing arithmetic.
  
          EXAMPLES::
  
              sage: QQ.is_exact()
              True
              sage: ZZ.is_exact()
              True
              sage: Qp(7).is_exact()
              False
              sage: Zp(7, type='capped-abs').is_exact()
              False
          """
          return True
  

  but the symbolic ring might have inexact stuff in it, so

  Definition:     SR.is_exact(self)
  Source:
      cpdef bint is_exact(self) except -2:
          """
          Return False, because there are approximate elements in the
          symbolic ring.
  
          EXAMPLES::
  
              sage: SR.is_exact()
              False
  
          Here is an inexact element.
  
          ::
  
              sage: SR(1.9393)
              1.93930000000000
          """
          return False
  

  does that make sense?

[jdemeyer]

Replying to kcrisman:

• Deprecation, folks.

What should be deprecated? Is there something which used to work and no longer works now and should be deprecated?

[kcrisman]

• Deprecation, folks.

What should be deprecated? Is there something which used to work and no longer works now and should be deprecated?

I just meant if one changed the name completely, that's all.

[jdemeyer]

Replying to vbraun:

Its somewhat confusing that is_numerical is completely different from GiNaC::numeric. How about is_approximate, is_limited_precision, or is_floating_point (I'd be happy to call padics "floating point" in this case).

I like none of those 3 alternatives, as there are valid cases for exact integers being "numerical": this isn't yet implemented on this ticket, but for number theoretical functions (binomial, factorial) it absolutely makes sense to consider integers "numerical". So the name should not obviously exclude integers.

In fact, this was one of the reasons to make _is_numerical a method of the function and not a global function (possible support of p-adics is indeed another reason).

[jdemeyer]

Replying to kcrisman:

I just meant if one changed the name completely, that's all.

If you show me a concrete instance where a deprecation is needed, I'll happily add it.

[kcrisman]

I just meant if one changed the name completely, that's all.

If you show me a concrete instance where a deprecation is needed, I'll happily add it.

So far it isn't! I was being a little preemptive in terms of the is_* guys. No worries.

[rws]

I have implemented #17151 on top of this patch. Could you please comment on if the code in Func_laguerre._evalf_() needed to determine real_parent can be reduced further, maybe by adding it to this patch, too?

[jdemeyer]

Replying to rws:

I have implemented #17151 on top of this patch. Could you please comment on if the code in Func_laguerre._evalf_() needed to determine real_parent can be reduced further

Most likely it can be reduced further (why not simply take the parent of x?). I also don't like the name real_parent since it can be complex too.

maybe by adding it to this patch, too?

No. People are already complaining that this patch is too big. But a follow-up ticket certainly makes sense.

[kcrisman]

Something else to note (very minor) is that if the top of ​http://trac.sagemath.org/wiki/symbolics/functions needed some slight updating in order to still be useful, it would be good to do that at the time this was resolved.

[jdemeyer]

Replying to kcrisman:

Something else to note (very minor) is that if the top of ​http://trac.sagemath.org/wiki/symbolics/functions needed some slight updating in order to still be useful, it would be good to do that at the time this was resolved.

I don't think anything needs updating. But it would be good though to think more generally about what the API should be for symbolic functions. I consider this ticket as a first step in the right direction.

[rws]

This is very useful for upcoming code of symbolic functions and contains fixes to functions too. I think I can follow what you do in symbolic. So apart from the doctest fails it's positive from me.

make ptestlong:

sage -t --long src/sage/modular/modform_hecketriangle/abstract_space.py  # 12 doctests failed
sage -t --long src/sage/modular/modform_hecketriangle/readme.py  # 7 doctests failed
sage -t --long src/sage/functions/other.py  # 3 doctests failed

---

New commits:

​382f97a	Merge branch 'u/jdemeyer/ticket/17130' of trac.sagemath.org:sage into 6.5beta1
​7265989	17130: reviewer's patch: fix typo

[jdemeyer]

New commits:

​c47dbd4	Fix Trac #17328 again in a better way
​a486db2	Call the factorial() method of an Integer

[rws]

Passes previously failed tests.

[rws]

Patchbot failed because of [tutorial ] IOError: [Errno 28] No space left on device

[vbraun]

I get more numerical noise

sage -t --long src/sage/functions/hyperbolic.py
**********************************************************************
File "src/sage/functions/hyperbolic.py", line 546, in sage.functions.hyperbolic.Function_arccoth.__init__
Failed example:
    float(arccoth(2))
Expected:
    0.5493061443340548
Got:
    0.5493061443340549
**********************************************************************
1 item had failures:
   1 of   9 in sage.functions.hyperbolic.Function_arccoth.__init__
    [156 tests, 1 failure, 1.22 s]
sage -t --long src/sage/functions/hypergeometric.py
    [143 tests, 5.23 s]
sage -t --long src/sage/functions/jacobi.py
    [188 tests, 4.02 s]
sage -t --long src/sage/functions/log.py
**********************************************************************
File "src/sage/functions/log.py", line 661, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(RDF(-1))
Expected:
    -0.3181315052047642 + 1.3372357014306895*I
Got:
    -0.31813150520476413 + 1.3372357014306895*I
**********************************************************************
File "src/sage/functions/log.py", line 663, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(float(-1))
Expected:
    (-0.3181315052047642+1.3372357014306895j)
Got:
    (-0.31813150520476413+1.3372357014306895j)
**********************************************************************
1 item had failures:
   2 of   9 in sage.functions.log.Function_lambert_w._evalf_
sage -t --long src/sage/functions/special.py
**********************************************************************
File "src/sage/functions/special.py", line 355, in sage.functions.special.MaximaFunction._evalf_
Failed example:
    f._evalf_(1, I, parent=CDF)
Expected:
    0.8483795707591759 - 0.0742924734216079*I
Got:
    0.8483795707591759 - 0.07429247342160791*I
**********************************************************************
1 item had failures:
   1 of  13 in sage.functions.special.MaximaFunction._evalf_
    [116 tests, 1 failure, 2.35 s]

[rws]

Replying to vbraun:

I get more numerical noise

Appears to be system dependent, as all tests in functions/ pass here (OpenSUSE 12.3)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​9d3cbbd

Fix numerical noise

[rws]

Still passes tests in functions/ here.

[vbraun]

sage -t --long src/sage/functions/log.py
**********************************************************************
File "src/sage/functions/log.py", line 661, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(RDF(-1))  # abs tol 1e-16
Expected:
    -0.31813150520476413 + 1.3372357014306895*I
Got:
    -0.31813150520476413 + 1.3372357014306893*I
Tolerance exceeded in 1 of 2:
    + 1.3372357014306895 vs + 1.3372357014306893, tolerance 2e-16 > 1e-16
**********************************************************************
File "src/sage/functions/log.py", line 663, in sage.functions.log.Function_lambert_w._evalf_
Failed example:
    lambert_w(float(-1))  # abs tol 1e-16
Expected:
    (-0.31813150520476413+1.3372357014306895j)
Got:
    (-0.31813150520476413+1.3372357014306893j)
Tolerance exceeded in 1 of 2:
    +1.3372357014306895 vs +1.3372357014306893, tolerance 2e-16 > 1e-16
**********************************************************************
1 item had failures:
   2 of   9 in sage.functions.log.Function_lambert_w._evalf_
    [171 tests, 2 failures, 2.51 s]

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​abab222

Fix more numerical noise

[rws]

There may be a followup in a later ticket to deal with an inconsistency, if you agree that the following discrepancy between functions with one and two parameters needs addressing:

sage: i = ComplexField(30).0
sage: gamma(i)
-0.15494983 - 0.49801567*I
sage: gamma(1.0*I)
-0.154949828301811 - 0.498015668118356*I
sage: jacobi_am(2., 2)
0.549820282407325
sage: jacobi_am(2., i)
jacobi_am(2.00000000000000, 1.0000000*I)
sage: jacobi_am(2., 1.0*I)
jacobi_am(2.00000000000000, 1.00000000000000*I)

EDIT: Ah, never mind, jacobi_am takes only reals.

[kcrisman]

Question for Jeroen and Ralf - any updating needed in the documentation in sage/symbolic/function_factory.py function function as a result of this ticket? I'd be happy to open a followup if need be.

[jdemeyer]

Replying to kcrisman:

Question for Jeroen and Ralf - any updating needed in the documentation in sage/symbolic/function_factory.py function function as a result of this ticket?

No, this is only about BuiltinFunctions, there are no functional changes to user-defined functions.