gh-40628: support ideals over finite fields where p > 2^29
    
<!-- ^ Please provide a concise and informative title. -->
Currently, doing `Ideal(...).variety()` on an ideal over a multivariate
polynomial over a finite field with p > 2^29, gives us this:

```
RuntimeError: error in Singular function call 'triangLfak':
not implemented for rings with rings as coeffients
error occurred in or before triang.lib::faktorisiere_letzten line 1091:
`        factors = factorize(h[nh],1);`
leaving triang.lib::faktorisiere_letzten (1087)
leaving triang.lib::triangLbas (262)
```

This PR uses the toy implementation of `.variety` in that cases

<!-- ^ 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". -->



### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #40628
Reported by: Arthur
Reviewer(s): Arthur, user202729, Vincent Macri

Author: Release Manager

src/sage/rings/polynomial/multi_polynomial_ideal.py

@@ -2760,6 +2760,21 @@ def _variety_triangular_decomposition(self, ring):
             sage: I = Ideal([ x*y - 1, (x-2)^2 + (y-1)^2 - 1])
             sage: len(I.variety())
             4
+
+        We can compute the variety of an Ideal over a Multivariate Polynomial Ring over a finite field with characteristic `> 2^{29}`, which Singular doesn't support. ::
+
+            sage: p = 626526183415513590854084217045148362725470879166437124330209
+            sage: F = GF(p)
+            sage: R.<x,y,z> = F[]
+            sage: I = Ideal([x^2 - 5*y + z, x*21 + y - z, 20*x + 20*y - 15*z + 20])
+            sage: I.variety()
+            verbose 0 (...: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation.
+            [{z: 475236874226935968499387140880743357810239093941490264140142,
+            y: 303497730986201757454241700121162099180641015844366285478588,
+            x: 366193016391757014347340764061969600539773744194969974791622},
+            {z: 151289309188577622354697076164405004915231785224946860207259,
+            y: 323028452429311833399842516923986263544829863322070838864298,
+            x: 260333167023756576506743452983178762185697134971467149538802}]
         """
 
         def _variety(T, V, v=None):
@@ -2802,10 +2817,14 @@ def _variety(T, V, v=None):
         P = self.ring()
         if ring is not None:
             P = P.change_ring(ring)
-        try:
-            TI = self.triangular_decomposition('singular:triangLfak')
-            T = [list(each.gens()) for each in TI]
-        except TypeError:  # conversion to Singular not supported
+        T = None
+        if P.characteristic() < 2**29:
+            try:
+                TI = self.triangular_decomposition('singular:triangLfak')
+                T = [list(each.gens()) for each in TI]
+            except TypeError:  # conversion to Singular not supported
+                pass
+        if T is None:
             if self.ring().term_order().is_global():
                 verbose("Warning: falling back to very slow toy implementation.", level=0, caller_name='variety')
                 T = toy_variety.triangular_factorization(self.groebner_basis())

src/sage/rings/polynomial/multi_polynomial_libsingular.pyx

@@ -4662,6 +4662,8 @@ cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
                 sig_off()
             finally:
                 if check_error():
+                    if len(self.variables()) == 1:
+                        return Factorization([(self._parent(p), e) for p, e in self.univariate_polynomial().factor()])
                     raise NotImplementedError("Factorization of multivariate polynomials over prime fields with characteristic > 2^29 is not implemented.")
 
             ivv = iv.ivGetVec()