Example:

restart;
myTable:=table(["a"=1,"b"=-1,"c"=1]);
select(i->myTable[i[]]=1,{indices(myTable)});

                                               {["a"], ["c"]}

Edit.

If we define  f  as a procedure, then we get a more convenient and compact syntax for specific calculations: 

restart;
with(orthopoly):
f:=unapply(convert(diff( P(n, cos(theta)), theta, theta),diff), n):
int( f(2), theta=0..Pi/4);
int( f(3), theta=0..Pi/4);

As an alternative you can use  diff  instead of  D :

restart;
g:=x->x^(1/3);
f:=x->convert(g(x),surd);
diff(f(x),x);

It's very simple:

restart;
A:=[1,2,6,7,9,10,15,17]:
A[3..5];

                                             [6, 7, 9]

restart;
f:=(a*x1+b*x2)/x3+c*x1^2/(e*x2+d);
normal(f);
numer(%);

Or

restart;
(a*x1+b*x2)/x3+c*x1^2/(e*x2+d)=f;
normal(%);
%*denom(lhs(%));

Code in Maple 2018.2 (I have no newer versions):

restart;
B:=sqrt( (-4*u^(1/3)+1)*u^(4/3));
A:=1/(-12*u+3*u^(2/3)-3*B);
res:=Int(A,u);
res1:=IntegrationTools:-Change(res,t=u^(1/3));
value(res1);
simplify(%) assuming t>0;
eval(%, t=u^(1/3)); 

You can get several more solutions from one solution if you use transformations that transform a cubic integer lattice into itself. They will be, for example, reflections about the coordinate planes, rotations around the coordinate axes at angles that are multiples of 90 degrees, and some others. But in terms of geometry, you end up with the same pair of vectors. New solutions can be easily obtained using the brute force method. First we find the list L and then from the list  L, the list  L1  is obtained, where none of the coordinates is equal to 0:

restart;
k:=0:
for a from 0 to 9 do
for b from a to 9 do
for c from b to 9 do
for x from -9 to 9 do
for y from -9 to 9 do
for z from -9 to 9 do
if 4*(a*x+b*y+c*z)^2=3*(a^2+b^2+c^2)*(x^2+y^2+z^2) and a^2+b^2+c^2<>0 and x^2+y^2+z^2<>0 and a*x+b*y+c*z>0 then k:=k+1;
L[k]:=[<a,b,c>,<x,y,z>] fi;
od: od: od: od: od: od:
L:=convert(L,list):  
nops(L);
L1:=select(t->`and`(seq(t[1][i]<>0,i=1..3),seq(t[2][i]<>0,i=1..3)),L):
nops(L1);

The very first pair of vectors in the lists  L  and  L1  is your original example, in which the first and third coordinates are reversed (reflection of space relative to the plane  x=z ). To reduce the search, I chose the vector  <a,b,c>  with non-decreasing coordinates.

This can be done in many ways. Here's one:

A := [1,0,3,4,5,5]: B := [0,1,2,6,3,6]:
nops(select(`>`, A - B, 0));

                                           3

The simple recursive procedure does the work:

restart;
a:=proc(n)
local m;
option remember;
if n=1 then return 1 elif n=2 then return 0 else
for m from 1 to n-2 do
if a(n-1)=a(n-1-m) then return m fi; od;
min(seq(abs(a(n-1)-a(k)),k=1..n-2)) fi;
end proc:

Examples of use:

seq(a(n), n=1..100);
plots:-listplot([%], style=point, symbol=solidcircle, color=red);

Edit. The procedure above works about 2 times faster than tomleslie's method.

In real domain the two-argument function  arctan(y,x)  returns the polar angle  alpha  of any point with the coordinates  (x,y)  on the plane in the range  -Pi < alpha <= Pi . Below are all the variants as it expressed through the standard function arctan(x) :

arctan(y,x) assuming x>=0,y>=0;
arctan(y,x) assuming x>=0,y<0;
arctan(y,x) assuming x<0,y>=0;
arctan(y,x) assuming x<0,y<0;

If I understand your problem correctly, then you want to get derivatives of different orders from the original expression  u[0]  using the formulas for differentiating the sum and the product at each step of the for-loop. See the following code for your specific example:

restart;
u[0] := exp(-x^2);
n:=10:
d := r->diff(coeff(r,u[0]),x)*u[0]+coeff(r,u[0])*diff(u[0],x):

for k from 1 to n do 
if type(u[k-1],`+`) then u[k] := map(d,u[k-1]) else 
u[k] := d(u[k-1]) end if; 
print(u[k]); 
end do:

Sx:=1/sqrt(2)*Matrix([[0,1,0],[1,0,1],[0,1,0]]);
lam,v:=LinearAlgebra:-Eigenvectors(Sx);
map(r->Vector(r/~LinearAlgebra:-VectorNorm(Vector(r),2)),convert(v^%T,listlist)); 

The problem is easy to solve if you just cut out regions with asymptotes from the plot:

restart;
func:=unapply(-tan((-4*t+x+2*y)/(2*sqrt(15)))/(2*sqrt(2)), x,y,t):
y1:=solve((-4*t+x+2*y)/(2*sqrt(15))=-Pi/2,y);
y2:=solve((-4*t+x+2*y)/(2*sqrt(15))=Pi/2,y);
y3:=solve((-4*t+x+2*y)/(2*sqrt(15))=3*Pi/2,y);
y4:=solve((-4*t+x+2*y)/(2*sqrt(15))=5*Pi/2,y);
A:=plot3d(func(x,y,4), x=-10..10, y=eval(y1+0.1..y2-0.1,t=4), style=surface):
B:=plot3d(func(x,y,4), x=-10..10, y=eval(y2+0.1..y3-0.1,t=4), style=surface):
C:=plot3d(func(x,y,4), x=-10..10, y=eval(y3+0.1..y4-0.1,t=4), style=surface):
plots:-display(A,B,C, scaling=constrained);

restart;
eq := cos(2*x)-sin(x)*(sin(x)^2+1)^(1/2)+cos(x)^2*sin(x)/(sin(x)^2+1)^(1/2):
solve(subs(x=arcsin(y), eq));

By default, Maple works in the complex domain and, of 3 values for a cube root, returns the value with the smallest argument (the principal value of a root). When working in the real domain, use the  surd  command:

a:=(-8)^(1/3);
evalc(a);
surd(-8,3);
plot(surd(x,3), x=-8..8, size=[1000,250], scaling=constrained, labels=[x,surd(x,3)]);

An alternative for  surd  would be to call the package  RealDomain :

with(RealDomain):
plot(x^(1/3), x=-8..8, size=[1000,250], scaling=constrained, labels=[x, x^(1/3)]);