Note that gmpy2’s factorial function requires integer arguments and returns a float result. Instead gmpy2 has a fac function that is analogous to math.factorial (or imath.factorial) and that computes and returns exact integer values. The interface of gmpy2 is like having imath, fractions, math and cmath all merged together in a single module so it is probably not a great comparison for what you would want the imath module’s interface to be.
Also I think aiming for compatibility with gmpy2 specifically is misguided in the sense that really what you want is ecosystem wide compatibility and ideally even compatibility beyond Python with comparable functions in other languages as well.
Let’s take math.comb as an example. There is gmpy2.comb which has the same name and is almost the same function but has a slightly different domain. I assume that the matching name here is an afterthought though because there is also an identical gmpy2.bincoef function.
The domain for math.comb(n, k) is nonnegative integer n and k but gmpy2.comb allows n to be negative. In mpmath and SymPy this function is called binomial and is defined for arbitrary complex n and k although specifically for negative integers n and k SymPy and mpmath disagree about the value of the function because there are two different conventions (I think mpmath has the right one).
SciPy calls the function comb and matches SymPy’s definition for negative integer values. The SciPy docs say that n and k have to be integers but it accepts and always returns floats unless exact=True is passed in which case it requires and returns ints. It doesn’t look like SciPy provides a complex extension of this function.
Looking beyond Python comb is called nchoosek in Matlab and requires nonnegative integers. In Mathematica and Maple the function is called binomial and is defined for complex arguments I think in a way that is consistent with mpmath. In Julia the function is called binomial. I might be using an outdated version of Julia here because the online docs say it is for Integer but I can pass float/complex for n but k needs to be an integer. For negative integer arguments Julia’s binomial matches the SymPy/SciPy convention rather than mpmath/Mathematica.
Ideally the authors of all Python libraries including the stdlib would consider this kind of compatibility when choosing a name, domain and definition of a mathematical function. You want the names to be consistent so that everyone agrees whether it is called e.g. atan or arctan. You also want the domains and definitions to be consistent as well though.
In the case of comb I think that means it should be called binomial although there is still the question of the negative integer convention. The math and gmpy2 decision to require nonnegative k might seem like it dodges that question but really it does not because it just creates a different kind of incompatibility that the different implementations of seemingly the same functions don’t have the same domains.
The case of atan vs arctan has to some extent been resolved as part of the array API since I think NumPy (and possibly imitators) was the only part of the Python ecosystem that used arctan and they have now all agreed to provide atan. It would be good to have some more standardisation/coordination around these things for functions that go beyond the current array API but I don’t think that choosing any one particular library like gmpy2 (or the stdlib or NumPy or …) is a good basis for defining any conventions.