From 4f2b298ead27c7cefc88ca585ef87be55ae9262e Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sat, 12 Sep 2026 20:33:48 +0300
Subject: [PATCH 1/2] Add test for crashing expr.n()

---
 symengine/tests/test_sympy_conv.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/symengine/tests/test_sympy_conv.py b/symengine/tests/test_sympy_conv.py
index 1d9790a7..598a5ff4 100644
--- a/symengine/tests/test_sympy_conv.py
+++ b/symengine/tests/test_sympy_conv.py
@@ -9,7 +9,8 @@
         KroneckerDelta, LeviCivita, erf, erfc, lowergamma, uppergamma,
         loggamma, beta, polygamma, sign, floor, ceiling, conjugate, And,
         Or, Not, Xor, Piecewise, Interval, EmptySet, FiniteSet, Contains,
-        Union, Complement, UniversalSet, Reals, Rationals, Integers)
+        Union, Complement, UniversalSet, Reals, Rationals, Integers,
+        PyFunction)
 import unittest
 
 # Note: We test _sympy_() for SymEngine -> SymPy conversion, as those are
@@ -854,3 +855,11 @@ def test_sympy_roundtrip():
     _check_sympy_roundtrip(x+y)
     _check_sympy_roundtrip(x**y)
     _check_sympy_roundtrip(d)
+
+
+@unittest.skipIf(not have_sympy, "SymPy not installed")
+def test_evalf_symbolic_pyfunction():
+    x = Symbol("x")
+    f = sympify(sympy.besselj(x, 1))
+    assert isinstance(f, PyFunction)
+    assert f.n() == f

From 0b76eaa88061a8ede02edccf1855765a0b5b84e4 Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Sat, 12 Sep 2026 20:34:19 +0300
Subject: [PATCH 2/2] Throw exception when PyFunction cannot be evalf-ed

---
 symengine/lib/pywrapper.cpp            | 23 ++++++++++++++++++++++-
 symengine/lib/symengine.pxd            |  1 +
 symengine/lib/symengine_wrapper.in.pyx |  7 ++++++-
 3 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/symengine/lib/pywrapper.cpp b/symengine/lib/pywrapper.cpp
index 5ab5e586..198fe70e 100644
--- a/symengine/lib/pywrapper.cpp
+++ b/symengine/lib/pywrapper.cpp
@@ -1,5 +1,6 @@
 #include "pywrapper.h"
 #include <symengine/serialize-cereal.h>
+#include <symengine/symengine_exception.h>
 
 #if PY_MAJOR_VERSION >= 3
 #define PyInt_FromLong PyLong_FromLong
@@ -243,7 +244,27 @@ RCP<const Basic> PyFunction::create(const vec_basic &x) const {
 }
 
 RCP<const Number> PyFunction::eval(long bits) const {
-    return pyfunction_class_->get_py_module()->eval_(pyobject_, bits);
+    RCP<const Number> result
+        = pyfunction_class_->get_py_module()->eval_(pyobject_, bits);
+    if (!result.is_null() && !PyErr_Occurred()) {
+        return result;
+    }
+
+    bool not_implemented = result.is_null();
+    if (PyErr_Occurred()) {
+        not_implemented
+            = PyErr_ExceptionMatches(PyExc_NotImplementedError) != 0;
+    }
+    PyErr_Clear();
+
+    if (not_implemented) {
+        // The Python function has no numerical value (e.g. it has free
+        // symbols). Signal the caller that it should stay unevaluated.
+        throw NotImplementedError(
+            "PyFunction cannot be evaluated to a number");
+    }
+
+    throw SymEngineException("Python function evaluation failed");
 }
 
 RCP<const Basic> PyFunction::diff_impl(const RCP<const Symbol> &s) const {
diff --git a/symengine/lib/symengine.pxd b/symengine/lib/symengine.pxd
index 81e88b2a..cc7cd26f 100644
--- a/symengine/lib/symengine.pxd
+++ b/symengine/lib/symengine.pxd
@@ -213,6 +213,7 @@ cdef extern from "<symengine/number.h>" namespace "SymEngine":
         pass
     cdef cppclass NumberWrapper(Basic):
         pass
+    bool is_a_Number(const Basic &b) nogil
     cdef tribool is_zero(const Basic &x) nogil
     cdef tribool is_positive(const Basic &x) nogil
     cdef tribool is_negative(const Basic &x) nogil
diff --git a/symengine/lib/symengine_wrapper.in.pyx b/symengine/lib/symengine_wrapper.in.pyx
index 552c9b74..0612f4fa 100644
--- a/symengine/lib/symengine_wrapper.in.pyx
+++ b/symengine/lib/symengine_wrapper.in.pyx
@@ -2900,8 +2900,13 @@ cdef PyObject* symengine_to_sympy(rcp_const_basic o1):
     return <PyObject*>(t)
 
 cdef RCP[const symengine.Number] sympy_eval(PyObject* o1, long bits):
+    cdef Basic X
     prec = max(1, int(round(bits/3.3219280948873626)-1))
-    cdef Number X = sympify((<object>o1).n(prec))
+    result = (<object>o1).n(prec)
+    X = sympify(result)
+    if not symengine.is_a_Number(deref(X.thisptr)):
+        raise NotImplementedError(
+            "expression cannot be evaluated to a number: %r" % (result,))
     return symengine.rcp_static_cast_Number(X.thisptr)
 
 cdef RCP[const symengine.Number] sage_eval(PyObject* o1, long bits):