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

Depends on #17131
Depends on #17133

CC: @burcin @rwst

Component: symbolics

Author: Jeroen Demeyer

Branch: abab222

Reviewer: Ralf Stephan

Issue created by migration from https://trac.sagemath.org/ticket/17130

[kcrisman]

comment:1

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]

comment:2

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]

comment:3

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]

comment:4

The following patch fixes the problem

diff --git a/src/sage/symbolic/function.pyx b/src/sage/symbolic/function.pyx
index 408b6da..5beb01d 100644
--- a/src/sage/symbolic/function.pyx
+++ b/src/sage/symbolic/function.pyx
@@ -914,38 +914,6 @@ cdef class BuiltinFunction(Function):
         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]

Description changed:

--- 
+++ 
@@ -11,3 +11,10 @@
 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
+```

[jdemeyer]

comment:6

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]

comment:7

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]

comment:8

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]

comment:9

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]

comment:10

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]

comment:11

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.