restart: with(LinearAlgebra): with(combinat):
#
# SetA337516 - compute set of resistances such that
#              A337516(n) = card(SetA337516(n))
#              Author Rainer Rosenthal October 2020
#
# Generators Series, Parallel, Bridge:
#
Ser := (x1,x2) -> x1 + x2:
Par := (x1,x2) -> 1 / (1/x1 + 1/x2):
Bri := (x1,x2,x3,x4,x5) -> (x2*x1*x4+x2*x1*x5+x5*x4*x1+x5*x4*x2+x3*(x2+x5)*(x1+x4))/
#                          --------------------------------------------------------
                                     ((x1+x2)*(x4+x5)+x3*(x1+x4+x2+x5)):
Gab := (R1,R2,R3,R4,R5,R6,R7) ->        (R1*R3*R4*R7+R1*R3*R4*R5+R1*R3*R2*R7+R1*R3*R2*R5+R2*R4*R3*R7+R2*R4*R3*R5+
        R2*R4*R1*R7+R2*R4*R1*R5+R5*R7*R1*R3+R5*R7*R1*R4+R5*R7*R2*R3+R5*R7*R2*R4+
        R6*R1*R3*R7+R6*R1*R3*R2+R6*R1*R3*R4+R6*R5*R7*R3+R6*R5*R2*R3+R6*R3*R4*R5+
        R6*R3*R4*R7+R6*R1*R4*R7+R6*R5*R7*R4+R6*R2*R4*R3+R6*R2*R4*R1+R6*R5*R2*R4)/
#   ---------------------------------------------------------------------------------
            (R3*R4*R7+R3*R4*R5+R2*R3*R7+R5*R2*R3+R1*R4*R7+R5*R1*R4+R1*R2*R7+
             R1*R2*R5+R5*R7*R3+R5*R7*R4+R5*R7*R1+R5*R7*R2+R6*R3*R7+R6*R2*R3+
             R6*R3*R4+R6*R1*R7+R6*R1*R2+R6*R1*R4+R6*R5*R7+R6*R5*R2+R6*R4*R5):                                     
#
# Derived set functions in this way:  SetFkt(A,B) = {Fkt(a,b) | a in A and b in B}
#
SetSer := (X1,X2) -> {seq(seq(
           Ser(X1[i1],X2[i2]),
                         i2=1..nops(X2)),
                  i1=1..nops(X1))}:
SetPar := (X1,X2) -> {seq(seq(
           Par(X1[i1],X2[i2]),
                         i2=1..nops(X2)),
                  i1=1..nops(X1))}:
SetBri := (X1,X2,X3,X4,X5) -> {seq(seq(seq(seq(seq(
           Bri(X1[i1],X2[i2],X3[i3],X4[i4],X5[i5]),
                                              i5=1..nops(X5)),
                                       i4=1..nops(X4)),
                                i3=1..nops(X3)),
                         i2=1..nops(X2)),
                  i1=1..nops(X1))}:
SetGab := (X1,X2,X3,X4,X5,X6,X7) -> {seq(seq(seq(seq(seq(seq(seq(
           Gab(X1[i1],X2[i2],X3[i3],X4[i4],X5[i5],X6[i6],X7[i7]),
                                                            i7=1..nops(X7)),
                                                     i6=1..nops(X6)),
                                              i5=1..nops(X5)),
                                       i4=1..nops(X4)),
                                i3=1..nops(X3)),
                         i2=1..nops(X2)),
                  i1=1..nops(X1))}:
SetA337516 := proc(n) option remember;
   local allres,k,w1,w2,xload,partload,i,addresi,
         permad,numad,wheremore,pw,increm,plan,ix,setsplanned;
   if n = 1 then return {1}; fi;
   allres := {}:
   allres := `union`(allres,seq(SetSer(SetA337516(k),SetA337516(n-k)),k=1..n/2));
   allres := `union`(allres,seq(SetPar(SetA337516(k),SetA337516(n-k)),k=1..n/2));
   #
   # Bridge combinations
   #
   xload := n - 5;     # 5 arms in a bridge
   if xload >= 0 then
      partload := partition(xload);
      for i to nops(partload) do
         addresi := partload[i];
         permad := permute(addresi);
         numad := nops(addresi);
         wheremore := choose(5,numad);
         for pw in wheremore do
            for increm in permad do
               plan := [1,1,1,1,1];
               setsplanned := [{1},{1},{1},{1},{1}];
               for ix to nops(pw) do
                  setsplanned[pw[ix]] := SetA337516(plan[pw[ix]] + increm[ix]);
               od;
               allres := allres union SetBri(op(setsplanned));
            od:
         od;
      od:
   fi;
   #
   # Gabel combinations
   #
   xload := n - 7;     # 7 arms in a Gabel
   if xload >= 0 then
      partload := partition(xload);
      for i to nops(partload) do
         addresi := partload[i];
         permad := permute(addresi);
         numad := nops(addresi);
         wheremore := choose(7,numad);
         for pw in wheremore do
            for increm in permad do
               plan := [seq(1,j=1..7)];
               setsplanned := [seq({1},j=1..7)];
               for ix to nops(pw) do
                  setsplanned[pw[ix]] := SetA337516(plan[pw[ix]] + increm[ix]);
               od;
               allres := allres union SetGab(op(setsplanned));
            od:
         od;
      od:
   fi;
   return allres;
end:
A337516 := n -> nops(SetA337516(n)):