gh-38721: Type of Z/nZ NTL polynomial evaluation should be scalar
    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->

When a polynomial over Z/n implemented in NTL is evaluated, the result
should be in Z/n.  Currently it's not; it's still an element of the
polynomial ring.  Example:

```
sage: R.<x> = PolynomialRing(Zmod(4), 'x', implementation='NTL')
sage: x.subs(1).parent()
Univariate Polynomial Ring in x over Ring of integers modulo 4 (using
NTL)
```

This patch fixes this behavior.
    
URL: #38721
Reported by: Kyle Hofmann
Reviewer(s): Kwankyu Lee

Author: Release Manager

src/sage/rings/polynomial/polynomial_modn_dense_ntl.pyx

@@ -1253,6 +1253,22 @@ cdef class Polynomial_dense_modn_ntl_zz(Polynomial_dense_mod_n):
             sage: S.<y> = PolynomialRing(Integers(5), implementation='NTL')
             sage: f(y)
             y^3 + 2
+
+        TESTS::
+
+            sage: R.<x> = PolynomialRing(Integers(100), implementation='NTL')
+            sage: f = x^3 + 7
+            sage: f(1).parent() == R.base_ring()
+            True
+            sage: f(int(1)).parent() == R.base_ring()
+            True
+            sage: f(x + 1).parent() == f.parent()
+            True
+
+            sage: R.<x> = PolynomialRing(Zmod(12), 'x', implementation='NTL')
+            sage: u = Zmod(4)(3)
+            sage: x(u).parent() == u.parent()
+            True
         """
         if len(args) != 1 or len(kwds) != 0:
             return Polynomial.__call__(self, *args, **kwds)
@@ -1273,7 +1289,7 @@ cdef class Polynomial_dense_modn_ntl_zz(Polynomial_dense_mod_n):
             return Polynomial.__call__(self, *args, **kwds)
         else:
             zz_pX_eval(fx.x, self.x, x.x)
-            return self._parent(int(fx))
+            return self._parent._base(int(fx))
 
 
 cdef class Polynomial_dense_modn_ntl_ZZ(Polynomial_dense_mod_n):
@@ -1808,6 +1824,25 @@ cdef class Polynomial_dense_modn_ntl_ZZ(Polynomial_dense_mod_n):
             sage: S.<y> = PolynomialRing(Integers(5), implementation='NTL')
             sage: f(y)
             y^3 + 2
+
+        TESTS::
+
+            sage: R.<x> = PolynomialRing(Integers(10^30), implementation='NTL')
+            sage: f = x^3 + 7
+            sage: f(1).parent() == R.base_ring()
+            True
+            sage: f(int(1)).parent() == R.base_ring()
+            True
+            sage: f(x + 1).parent() == f.parent()
+            True
+
+            sage: R.<x> = PolynomialRing(Zmod(10^30), 'x', implementation='NTL')
+            sage: u = Zmod(10^29)(3)
+            sage: x(u).parent() == u.parent()
+            True
+            sage: v = Zmod(10)(3)
+            sage: x(v).parent() == v.parent()
+            True
         """
         if len(args) != 1 or len(kwds) != 0:
             return Polynomial.__call__(self, *args, **kwds)
@@ -1826,7 +1861,7 @@ cdef class Polynomial_dense_modn_ntl_ZZ(Polynomial_dense_mod_n):
             return Polynomial.__call__(self, *args, **kwds)
         else:
             ZZ_pX_eval(fx.x, self.x, x.x)
-            return self._parent(fx._integer_())
+            return self._parent._base(fx._integer_())
 
 
 cdef class Polynomial_dense_mod_p(Polynomial_dense_mod_n):