gh-40532: Implement compositional inverse
    
fix #40286


### 📝 Checklist

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

- [ ] The title is concise and informative.
- [ ] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [ ] 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: #40532
Reported by: user202729
Reviewer(s): Martin Rubey, Travis Scrimshaw, user202729

Author: Release Manager

src/sage/modules/free_module_element.pyx

@@ -4271,6 +4271,50 @@ cdef class FreeModuleElement(Vector):   # abstract base class
             return vector(ring, coeffs)
         return vector(coeffs)
 
+    def compositional_inverse(self, allow_multivalued_inverse=True, **kwargs):
+        """
+        Find the compositional inverse of this symbolic function.
+
+        INPUT: see :meth:`sage.symbolic.expression.Expression.compositional_inverse`
+
+        .. SEEALSO::
+
+            :meth:`sage.symbolic.expression.Expression.compositional_inverse`
+
+        EXAMPLES::
+
+            sage: f(x, y, z) = (y, z, x)
+            sage: f.compositional_inverse()
+            (x, y, z) |--> (z, x, y)
+
+        TESTS::
+
+            sage: f.change_ring(SR)
+            (y, z, x)
+            sage: f.change_ring(SR).compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: base ring must be a symbolic expression ring
+        """
+        from sage.rings.abc import CallableSymbolicExpressionRing
+        if not isinstance(self.base_ring(), CallableSymbolicExpressionRing):
+            raise ValueError("base ring must be a symbolic expression ring")
+        from sage.symbolic.ring import SR
+        tmp_vars = [SR.symbol() for _ in range(self.parent().dimension())]
+        input_vars = self.base_ring().args()
+        from sage.symbolic.relation import solve
+        l = solve([a == b for a, b in zip(self.change_ring(SR), tmp_vars)], input_vars, solution_dict=True, **kwargs)
+        if not l:
+            raise ValueError("cannot find an inverse")
+        if len(l) > 1 and not allow_multivalued_inverse:
+            raise ValueError("inverse is multivalued, pass allow_multivalued_inverse=True to bypass")
+        d = l[0]
+        subs_dict = dict(zip(tmp_vars, input_vars))
+        for x in input_vars:
+            if set(d[x].variables()) & set(input_vars):
+                raise ValueError("cannot find an inverse")
+        return self.parent()([d[x].subs(subs_dict) for x in input_vars])
+
 
 # ############################################
 # Generic dense element

src/sage/symbolic/expression.pyx

@@ -13443,6 +13443,83 @@ cdef class Expression(Expression_abc):
             raise TypeError("this expression must be a relation")
         return self / x
 
+    def compositional_inverse(self, allow_multivalued_inverse=True, **kwargs):
+        """
+        Find the compositional inverse of this symbolic function.
+
+        INPUT:
+
+        - ``allow_multivalued_inverse`` -- (default: ``True``); see example below
+        - ``**kwargs`` -- additional keyword arguments passed to :func:`sage.symbolic.relation.solve`.
+
+        .. SEEALSO::
+
+            :meth:`sage.modules.free_module_element.FreeModuleElement.compositional_inverse`.
+
+        EXAMPLES::
+
+            sage: f(x) = x+1
+            sage: f.compositional_inverse()
+            x |--> x - 1
+            sage: var("y")
+            y
+            sage: f(x) = x+y
+            sage: f.compositional_inverse()
+            x |--> x - y
+            sage: f(x) = x^2
+            sage: f.compositional_inverse()
+            x |--> -sqrt(x)
+
+        When ``allow_multivalued_inverse=False``, there is some additional checking::
+
+            sage: f(x) = x^2
+            sage: f.compositional_inverse(allow_multivalued_inverse=False)
+            Traceback (most recent call last):
+            ...
+            ValueError: inverse is multivalued, pass allow_multivalued_inverse=True to bypass
+
+        Nonetheless, the checking is not always foolproof (``x |--> log(x) + 2*pi*I`` is another possibility)::
+
+            sage: f(x) = exp(x)
+            sage: f.compositional_inverse(allow_multivalued_inverse=False)
+            x |--> log(x)
+
+        Sometimes passing ``kwargs`` is useful, for example ``algorithm`` can be used
+        when the default solver fails::
+
+            sage: f(x) = (2/3)^x
+            sage: f.compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            KeyError: x
+            sage: f.compositional_inverse(algorithm="giac")                             # needs sage.libs.giac
+            x |--> -log(x)/(log(3) - log(2))
+
+        TESTS::
+
+            sage: f(x) = x+exp(x)
+            sage: f.compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: cannot find an inverse
+            sage: f(x) = 0
+            sage: f.compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: cannot find an inverse
+            sage: f(x, y) = (x, x)
+            sage: f.compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: cannot find an inverse
+            sage: (x+1).compositional_inverse()
+            Traceback (most recent call last):
+            ...
+            ValueError: base ring must be a symbolic expression ring
+        """
+        from sage.modules.free_module_element import vector
+        return vector([self]).compositional_inverse(allow_multivalued_inverse=allow_multivalued_inverse, **kwargs)[0]
+
     def implicit_derivative(self, Y, X, n=1):
         """
         Return the `n`-th derivative of `Y` with respect to `X` given