For example, here are two equations containing trigonometric functions (Note that they do not form one system!): 

restart; # There are more examples, yet for the sake of briefness, they are omitted here. 
eqn__0 := cos(x)*cos(y)*cos(x + y) = 2*(sin(x)*sin(y) - 1)*2*(sin(x)*sin(x + y) - 1)*2*(sin(y)*sin(x + y) - 1):
eqn__1 := (cos(x + y) - (cos(x) + cos(y)) + 1)**2 + 2*cos(x)*cos(y)*cos(x + y) = 0:

Unfortunately, none of 

(* Tag０ *) RealDomain:-solve(eqn__0, {y, x}):
(* Tag１ *) solve(eqn__0, {y, x}) assuming y + x >= 0, (y, x) <=~ Pi:
(* Tag２ *) RealDomain:-solve(eqn__1, {y, x}):
(* Tag３ *) solve(eqn__1, {y, x}) assuming y + x >= 0, (y, x) <=~ Pi:

outputs concise solutions.
Using `plot3d`, it is easy to check that when "And(y + x >= 0, (y, x) <=~ Pi)", “{y = Pi/2, x = 0}, {y = Pi/3, x = Pi/3}, {y = 0, x = Pi/2}, {y = Pi/2, x = Pi/2}” is both the only solution to "eqn__0" and the only solution to "eqn__1". But how to get Maple to do so without manual intervention?

Edit. The main purpose is to automatically find the generic solutions to each of the two equations (Tag０ and Tag２) (separately). Now that the cosine and sine functions are both periodic with period 2π and both (lhs - rhs)(eqn__0) and (lhs - rhs)(eqn__1) are even symmetric, it is enough to focus only on the region y + x ≥ 0 ∧ (y, x) ≤~ Pi. So, in theory, a second-best workaround should be Tag１ and Tag３. However, why is Maple still unable to find the four exact solutions above?

The help page mentions: 

But what about the so-called q-hypergeometric function? Though there exist QDifferenceEquations:-IsQHypergeometricTerm and QDifferenceEquations:-QHypergeometricSolution in Maple, they do not seem to represent the function itself. 
For example, how to type the q-Gauss sum (cf. DLMF's §17.6(i)) or verify the last “simple series expression” given in Basic hypergeometric series - Wikipedia? In Mma, one may achieve these with something like 

, 

while 

convert("QHypergeometricPFQ[{a, b}, {c}, q, c/(a b)]", 'FromMma', 'evaluate');
 = 
                              /                 c \
            QHypergeometricPFQ|[a, b], [c], q, ---|
                              \                a b/

So has the q-hypergeometric function been implemented in Maple?

I would like to simulate the evolution of the so-called B, C, K, W system and SKI combinator calculus in Maple.
The rewrite rules of them are simple: 

However, if I try to evaluate the following example given in the  article, 

Maple will only return an unchanged result: 

restart;
applyrule([cS(x::anything)(y::anything)(z::anything) = x(z)(y(z)), 
   cK(x::anything)(y::anything) = x, cI(x::anything) = x], 
  cS(cK(cS(cI)))(cS(cK(cK))(cI))(x)(y)); # Unable to reduce??? 
 = 
              cS(cK(cS(cI)))(cS(cK(cK))(cI))(x)(y)

I believe that this is not an outlier.
Here are two additional instances: 

> rls := [cS(x::anything)(y::anything)(z::anything) = x(z)(y(z)), cK(x::anything)(y::anything) = x]:
> map2(applyrule, rls, [cS(cS(cS)(cS))(cS)(cS(cS))(cK), cS(cS(cS))(cS)(cS)(cS)(cS(cS)(cK(cK)))]);
                       [cS(cS(cS)(cS))(cS)(cS(cS))(cK), cS(cS(cS))(cS)(cS)(cS)(cS(cS)(cK(cK)))]

>

So why can't `applyrule` apply rules as desired? Meanwhile, how to automatically and thoroughly (like :-eval['recurse'] or MmaTranslator:-Mma:-ReplaceRepeated) apply those transformation rules to 

1. , 
2. , and
3. ? 

I have read something like How to apply a recursive rule in an expression? - MaplePrimes, but they are not the same issue.

According to the documentation of MmaTranslator:-Mma:-PolynomialReduce, this command yields . However, 

restart;
MmaTranslator:-Mma:-PolynomialReduce(x**2+y**2,{x-y,y+a});
 = 
                       [         2    2]
                       [[0, 0], x  + y ]

In[1]:= PolynomialReduce[x^2+y^2,{x-y,y+a}](*Mathematica*)

Out[1]= {{x + y, -2 a + 2 y}, 2 a^2}

In SymPy and in MuPAD: 

The output of both is the same as that of Mma; only the result given by Maple is inconsistent with Mathematica's. 

The example above is so simple that the desired result can be found simply by hand. Here is a larger example: 
Given two polynomials .txt and .txt, as well as a list of polynomials .txt, I would like to evaluate 

# Suppose that one has downloaded these three files. 
poly1, poly2 := fscanf("poly1.txt", "%a")[], fscanf("poly2.txt", "%a")[]:
pList := MmaTranslator:-Mma:-ReadList("pList.txt"):
MmaTranslator:-Mma:-PolynomialReduce((a - poly1)*(a - poly2), pList);

 But its result is just “[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], …]”, while when a＝0 it should be “[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 1, 1, 1, 1, 0, 0, 0, 1, 2, 1, 0, 2, 2, 3, 1, 1, 1, 2, 1, 0, 0, 0, 1, 1, 2, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0], 0]”.
So why does  return a distinct value?

In an old question, @mbras asked for a "partial" `convert/elsymfun`. However, SymPy's sympy.polys.rings.PolyElement.symmetrize seems to provide more examples that cannot be handled by the program that appeared in that question.
For instance, 

>>> from sympy import var
>>> var('x:z,p:r')
(x, y, z, p, q, r)
>>> from sympy.polys.polyfuncs import symmetrize
>>> symmetrize(x**2-(y**2+2**z),[y,x],formal=True,symbols=[p+p,q*q])[0]
-2**z - 4*p**2 + 2*q**2
>>> symmetrize(x*x*y+y*y*z+z*z*x,[y,x,z],formal=True,symbols=[p,q,r])
(0, x**2*y + x*z**2 + y**2*z, [(p, x + y + z), (q, x*y + x*z + y*z), (r, x*y*z)])

Though I can , can't the built-in  be generalized to such expressions (in other words, write the polynomial part of input as a symmetric part and a remainder with (named, if need be) elementary symmetric polynomials)?

Besides, since any symmetric polynomial can also be expressed in terms of the complete symmetric polynomials, is there a similar  command in Maple?