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]

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.

[jdemeyer]

Author: Jeroen Demeyer

[jdemeyer]

Branch: u/jdemeyer/ticket/17130

[jdemeyer]

New commits:

6d10729

Simplify numerical evaluation of BuiltinFunctions

[jdemeyer]

Commit: 6d10729

[jdemeyer]

Dependencies: #17131

[jdemeyer]

Changed dependencies from #17131 to #17131, #17133

[sagetrac-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

[sagetrac-git]

Changed commit from 6d10729 to b6e1ed4

[kcrisman]

comment:21

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]

comment:22

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]

comment:23

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.

[sagetrac-git]

Changed commit from a486db2 to 9d3cbbd

[sagetrac-git]

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

9d3cbbd

Fix numerical noise

[rwst]

comment:51

Still passes tests in functions/ here.

[vbraun]

comment:52

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]

[sagetrac-git]

Changed commit from 9d3cbbd to abab222

[sagetrac-git]

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

abab222

Fix more numerical noise

[rwst]

comment:55

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.

[vbraun]

Changed branch from u/jdemeyer/ticket/17130 to abab222

[kcrisman]

Changed commit from abab222 to none

[kcrisman]

comment:57

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]

comment:58

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.