@Carl Love Many thanks.  I changed your example to one of my problems.  I think I got it right but I would like to hear from you.

parms := [alpha[1, 0], alpha[1, 1], alpha[1, 2], alpha[1, 3], alpha[1, 4], alpha[1, 5], alpha[1, 6], alpha[1, 7], alpha[1, 8], alpha[1, 9], alpha[2, 0], alpha[2, 1], alpha[2, 2], alpha[2, 3], alpha[2, 4], alpha[2, 5], alpha[2, 6], alpha[2, 7], alpha[2, 8], alpha[2, 9], alpha[3, 0], alpha[3, 1], alpha[3, 2], alpha[3, 3], alpha[3, 4], alpha[3, 5], alpha[3, 6], alpha[3, 7], alpha[3, 8], alpha[3, 9]]

S:= parms:
B:=s-> evalb(numelems(convert(map2(op, 1, s),set))=3):
n:=4:
# The Array method:
A:= Array(1..0): 
for C in Iterator:-Combination(nops(S),n) do
    if B((s:= S[[seq(C+~1)]])) then A,= s fi
od:
L1:= [seq(A)]:
# The embedded loop method: 
L2:= [
    for C in Iterator:-Combination(nops(S), n) do
        if B((s:= S[[seq(C+~1)]])) then s fi
    od
]:

Could you tell me what ",=" means? Many thanks. 

@Carl Love Many thanks.  Please give some time to completely understand your solution.   Meanwhile, can you give away the major syntax additions? I am also running Maple 2020.  

The lines 

B1:= (s::set)-> `if`(B(s), s, NULL);
L:= [seq(B1(S[[seq(C+~1)]]), C= Iterator:-Combination(nops(S), 6))]

don't seem to work. They return

L:=[ifelse(B({22424170465, 193139816415, 386408307450, 395718860534, 412286285840, 427552056869}), {22424170465, 193139816415, 386408307450, 395718860534, 412286285840, 427552056869}, NULL) .....

@Carl Love Many thanks.  Two questions if I may:

1) I don't know I get it right, but can you do the same using seq (option 2)?

2) Suppose I want some save some of them using a certain criterion.   Would it faster to save them to a table?  (I know you answered a question a few years ago showing the difference between list, table and (I completely forget) but I can't find it). Many thanks.

@Kitonum That phrase was referring to the use of subsets.   The main point is to get the subsets with a fixed number of elements one by one without generating them all up front.   When I sent the question I was not aware that nextcomb was available.  Iterator seems a good option too.  Many thanks. 

@acer Many thanks.  That is exactly what I need.   It seems that Iterator:-Combination is another good option.   

@Kitonum Many thanks but let us consider combinat:-numbcomb(30,10) = 30045015.  Depending on my previous calculations I will have to test them all and I am not sure if they will fit the memory.  

@acer Since you opened the door, here is the problem.  I will take the 3D quadratic form as an illustration.

f := [x^2*alpha[1, 4]+x*y*alpha[1, 5]+x*z*alpha[1, 6]+y^2*alpha[1, 7]+y*z*alpha[1, 8]+z^2*alpha[1, 9]+x*alpha[1, 1]+y*alpha[1, 2]+z*alpha[1, 3]+alpha[1, 0], x^2*alpha[2, 4]+x*y*alpha[2, 5]+x*z*alpha[2, 6]+y^2*alpha[2, 7]+y*z*alpha[2, 8]+z^2*alpha[2, 9]+x*alpha[2, 1]+y*alpha[2, 2]+z*alpha[2, 3]+alpha[2, 0], x^2*alpha[3, 4]+x*y*alpha[3, 5]+x*z*alpha[3, 6]+y^2*alpha[3, 7]+y*z*alpha[3, 8]+z^2*alpha[3, 9]+x*alpha[3, 1]+y*alpha[3, 2]+z*alpha[3, 3]+alpha[3, 0]]

To test the symmetry (another problem) I need to swap y with z and vice-versa. That was the problem in the head of this thread.  Suppose now the test for the symmetry y-z in the 3D cubic form 

f1 := [x^3*alpha[1, 10]+x^2*y*alpha[1, 11]+x^2*z*alpha[1, 12]+x*y^2*alpha[1, 13]+x*y*z*alpha[1, 14]+x*z^2*alpha[1, 15]+y^3*alpha[1, 16]+y^2*z*alpha[1, 17]+y*z^2*alpha[1, 18]+z^3*alpha[1, 19]+x^2*alpha[1, 4]+x*y*alpha[1, 5]+x*z*alpha[1, 6]+y^2*alpha[1, 7]+y*z*alpha[1, 8]+z^2*alpha[1, 9]+x*alpha[1, 1]+y*alpha[1, 2]+z*alpha[1, 3]+alpha[1, 0], x^3*alpha[2, 10]+x^2*y*alpha[2, 11]+x^2*z*alpha[2, 12]+x*y^2*alpha[2, 13]+x*y*z*alpha[2, 14]+x*z^2*alpha[2, 15]+y^3*alpha[2, 16]+y^2*z*alpha[2, 17]+y*z^2*alpha[2, 18]+z^3*alpha[2, 19]+x^2*alpha[2, 4]+x*y*alpha[2, 5]+x*z*alpha[2, 6]+y^2*alpha[2, 7]+y*z*alpha[2, 8]+z^2*alpha[2, 9]+x*alpha[2, 1]+y*alpha[2, 2]+z*alpha[2, 3]+alpha[2, 0], x^3*alpha[3, 10]+x^2*y*alpha[3, 11]+x^2*z*alpha[3, 12]+x*y^2*alpha[3, 13]+x*y*z*alpha[3, 14]+x*z^2*alpha[3, 15]+y^3*alpha[3, 16]+y^2*z*alpha[3, 17]+y*z^2*alpha[3, 18]+z^3*alpha[3, 19]+x^2*alpha[3, 4]+x*y*alpha[3, 5]+x*z*alpha[3, 6]+y^2*alpha[3, 7]+y*z*alpha[3, 8]+z^2*alpha[3, 9]+x*alpha[3, 1]+y*alpha[3, 2]+z*alpha[3, 3]+alpha[3, 0]]

T1 and T2 need to be updated in the routines.   I haven't given much thought on that now but it will be in my mind soon. 

@acer Many many thanks.  I will certainly take into account your advice.  Since my maple proficiency is as bad as you could see by the questions I sent to the list,  I am afraid that I have a long way to find the cleanest and fastest solution to the problem.

@Carl Love Many thanks.  I need help to understand the last line.  specindex is barely touched in the help file (I am running 2017 and 2020).  And the op part as well. 

@acer Many many thanks.  I have thousands of lists like the one in the example.  That one is just to prepare the set of coefficients to be tested for the symmetry y->z in a full 3D quadratic model.  It is my plan to extend that to a higher degree and a higher dimension.  The problem will be now to build T1 and T2 programmatically for all cases.  

@Carl Love Many many thanks. I am carefully checking the line after varA. I know what is going on but I think I will never come up with a solution like that.  

@Carl I found the mistake.   It works flawlessly.  Would you have a shorter solution?   That is, no need to apply your solution four times. 

@Carl Love Many thanks.   It worked but please my answer to @Joe Riel.   I could not figure out how to change your solution to accommodate the remaining substitutions.  I have to confess that I did not completely understand it.  For the extra 3 steps I did:

1st step - Back to A.

T:= table([0=0,1=1,2=2,3=3,4=4,5=5,6=6,7=7,8=8,9=9]):
newA1:=subsindets(
    newA, 
    `?[]`~(B, `[]`~(anything, {indices(T,'nolist')})), 
    a-> A[op(1,a), T[op(2,a)]]
);

2nd step

T:= table([2= 3, 3= 2]):
newA2:=subsindets(
    newA1, 
    `?[]`~(A, `[]`~(anything, {indices(T,'nolist')})), 
    a-> B[T[op(1,a)], op(2,a)]
);

 almost gave me the right answer.  The first B terms were messed up.  

@Joe Riel Many thanks.  It seems that you guessed the next step but I need to apply two more substitutions so the final output is

{Subscript[A, 1, 0], Subscript[A, 1, 1], Subscript[A, 1, 3], \
Subscript[A, 1, 2], Subscript[A, 1, 4], Subscript[A, 1, 6], \
Subscript[A, 1, 5], Subscript[A, 1, 9], Subscript[A, 1, 8], \
Subscript[A, 1, 7], Subscript[A, 3, 0], Subscript[A, 3, 1], \
Subscript[A, 3, 3], Subscript[A, 3, 2], Subscript[A, 3, 4], \
Subscript[A, 3, 6], Subscript[A, 3, 5], Subscript[A, 3, 9], \
Subscript[A, 3, 8], Subscript[A, 3, 7], Subscript[A, 2, 0], \
Subscript[A, 2, 1], Subscript[A, 2, 3], Subscript[A, 2, 2], \
Subscript[A, 2, 4], Subscript[A, 2, 6], Subscript[A, 2, 5], \
Subscript[A, 2, 9], Subscript[A, 2, 8], Subscript[A, 2, 7]}

the extra 3 steps in Mathematica to get the result shown above are:

Subscript[B, m_, n_] -> Subscript[A, m, n] 
{Subscript[A, 2, m_] -> Subscript[B, 3, m],  Subscript[A, 3, m_] -> Subscript[B, 2, m]}
Subscript[B, m_, n_] -> Subscript[A, m, n]

With your solution, I do not need the first step.

I tried to modify the code to do the second step but it didn't work (that means that I did not understand what you did). 

@Carl Love  Many thanks.  For the example shown above the output should be

{Subscript[A, 1, 0], Subscript[A, 1, 1], Subscript[B, 1, 3], \
Subscript[B, 1, 2], Subscript[A, 1, 4], Subscript[B, 1, 6], \
Subscript[B, 1, 5], Subscript[B, 1, 9], Subscript[A, 1, 8], \
Subscript[B, 1, 7], Subscript[A, 2, 0], Subscript[A, 2, 1], \
Subscript[B, 2, 3], Subscript[B, 2, 2], Subscript[A, 2, 4], \
Subscript[B, 2, 6], Subscript[B, 2, 5], Subscript[B, 2, 9], \
Subscript[A, 2, 8], Subscript[B, 2, 7], Subscript[A, 3, 0], \
Subscript[A, 3, 1], Subscript[B, 3, 3], Subscript[B, 3, 2], \
Subscript[A, 3, 4], Subscript[B, 3, 6], Subscript[B, 3, 5], \
Subscript[B, 3, 9], Subscript[A, 3, 8], Subscript[B, 3, 7]}

if I change varA to {Subscript[A, 1, 9], Subscript[A, 2, 9], Subscript[A, 3, 9]}, the output is {Subscript[B, 1, 7], Subscript[B, 2, 7], Subscript[B, 3, 7]}