text:
local reps,tori,labels,I,II,ords;

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],              # E_8
#       [[1,1],[2,1],[10,1],[11,1],[5,1],[6,1],[7,1],[8,1]],            # E_8(a_1)
       [[1,1],[2,1],[10,1],[11,1],[5,1],[6,1],[7,1],[8,1]],            # E_8(a_1)
#       [[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1],[8,1]],          # E_8(a_2)
       [[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1],[8,1]],          # E_8(a_2)
#       [[8,1],[18,1],[7,1],[6,1],[19,1],[1,1],[17,1],[4,1]],           # E_8(a_3)
       [[7,1],[8,1],[9,1],[10,1],[11,1],[12,1],[13,1],[19,1]],         # E_8(a_3)
#       [[18,1],[6,1],[19,1],[1,1],[17,1],[21,1],[4,1],[8,1]],          # E_8(a_4)
       [[9,1],[10,1],[11,1],[12,1],[13,1],[14,1],[15,1],[19,1]],       # E_8(a_4)
       [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1]],                    # E_7
#       [[8,1],[7,1],[14,1],[26,1],[19,1],[1,1],[17,1],[4,1]],          # E_8(b_4)
       [[7,1],[8,1],[9,1],[17,1],[18,1],[19,1],[20,1],[21,1]],         # E_8(b_4)
#       [[18,1],[22,1],[19,1],[1,1],[17,1],[20,1],[21,1],[4,1],[6,1]],  # E_8(a_5)
       [[9,1],[14,1],[15,1],[17,1],[18,1],[19,1],[20,1],[26,1]],       # E_8(a_5)
#       [[1,1],[2,1],[11,1],[10,1],[5,1],[6,1],[7,1]],                  # E_7(a_1)
       [[1,1],[3,1],[10,1],[11,1],[5,1],[6,1],[7,1]],                  # E_7(a_1)
#       [[18,1],[8,1],[31,1],[7,1],[33,1],[4,1],[19,1],[23,1]],         # E_8(b_5)
       [[8,1],[16,1],[17,1],[18,1],[20,1],[21,1],[24,1],[33,1]],       # E_8(b_5)
       [[8,1],[7,1],[6,1],[5,1],[4,1],[2,1],[3,1]],                    # D_7
#       [[18,1],[31,1],[7,1],[33,1],[4,1],[29,1],[23,1],[19,1]],        # E_8(a_6)
       [[4,1],[10,1],[12,1],[14,1],[15,1],[26,1],[27,1],[30,1]],       # E_8(a_6)
#       [[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],                # E_7(a_2)
       [[1,1],[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],                # E_7(a_2)
       [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[8,1]],                    # E_6A_1

#      [[2,1],[3,1],[10,1],[12,1],[6,1],[7,1],[8,1],[83,1]],            # (D_7(a_1))^{(2)}    char 2
#       [[8,1],[7,1],[6,1],[5,1],[10,1],[11,-1],[2,1],[3,1]],            # D_7(a_1)
       [[2,1],[3,1],[10,1],[12,1],[6,1],[7,1],[8,1]],                   # D_7(a_1)
#       [[31,1],[33,1],[18,1],[4,1],[39,-1],[41,1],[8,1],[19,1],[23,1]], # E_8(b_6)
       [[15,1],[22,1],[23,1],[24,1],[25,1],[26,1],[27,1],[28,1]],       # E_8(b_6)
#       [[17,1],[1,1],[19,1],[6,1],[7,1],[18,1],[4,1]],                  # E_7(a_3)
       [[1,1],[4,1],[7,1],[12,1],[13,1],[17,1],[19,1]],                 # E_7(a_3)
#       [[3,1],[5,1],[1,1],[6,1],[11,1],[12,1],[2,1],[8,1]],             # E_6(a_1)A_1
       [[1,1],[2,1],[11,1],[12,1],[5,1],[6,1],[8,1]],                   # E_6(a_1)+A_1
#       [[1,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1],[67,1]],             # (A_7)^{(3)}     char 3
       [[1,1],[3,1],[4,1],[5,1],[6,1],[7,1],[8,1]],                     # A_7
#       [[8,1],[7,1],[13,1],[10,1],[11,-1],[12,1],[2,1],[3,1]],          # D_7(a_2)
       [[2,1],[3,1],[10,1],[12,1],[13,1],[14,1],[8,1]],                 # D_7(a_2)
       [[1,1],[2,1],[3,1],[4,1],[5,1],[6,1]],                           # E_6
       [[7,1],[6,1],[5,1],[4,1],[2,1],[3,1]],                           # D_6
#       [[1,1],[2,1],[3,1],[4,1],[5,1],[7,1],[8,1],[51,1]],             # (D_5+A_2)^{(2)}     char 2
       [[1,1],[3,1],[4,1],[2,1],[5,1],[7,1],[8,1]],                     # D_5A_2
#       [[3,1],[5,1],[1,1],[6,1],[11,1],[12,1],[2,1]],                   # E_6(a_1)
       [[1,1],[2,1],[11,1],[12,1],[5,1],[6,1]],                         # E_6(a_1)
#       [[17,1],[1,1],[19,1],[14,1],[26,1],[7,1],[18,1],[4,1]],          # E_7(a_4)
       [[1,1],[4,1],[7,1],[11,1],[25,1],[26,1],[27,1]],                 # E_7(a_4)
       [[2,1],[4,1],[5,1],[6,1],[7,1],[8,1],[1,1]],                     # A_6A_1
#       [[7,1],[6,1],[5,1],[10,1],[11,-1],[2,1],[3,1]],                  # D_6(a_1)
       [[2,1],[3,1],[10,1],[12,1],[6,1],[7,1]],                         # D_6(a_1)
       [[2,1],[4,1],[5,1],[6,1],[7,1],[8,1]],                           # A_6
#       [[47,1],[49,1],[5,1],[45,1],[44,1],[42,1],[27,1],[46,1]],        # E_8(a_7)
       [[30,1],[33,1],[35,1],[36,1],[40,1],[45,1],[53,1],[58,1]],       # E_8(a_7)

       [[1,1],[3,1],[4,1],[2,1],[5,1],[7,1]],                # D_5A_1
#       [[24,1],[26,1],[7,1],[27,1],[23,1],[4,1],[25,1]],     # E_7(a_5)
       [[4,1],[11,1],[14,1],[16,1],[24,1],[25,1],[26,1]],    # E_7(a_5)
#       [[17,1],[1,1],[19,1],[6,1],[18,1],[4,1],[8,1]],       # E_6(a_3)A_1
       [[4,1],[6,1],[9,1],[11,1],[18,1],[19,1],[8,1]],       # E_6(a_3)+A_1
#       [[7,1],[13,1],[10,1],[11,-1],[12,1],[2,1],[3,1]],     # D_6(a_2)
       [[2,1],[3,1],[10,1],[12,1],[13,1],[14,1]],            # D_6(a_2)
#       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1],[7,1],[8,1]],  # D_5(a_1)A_2
       [[1,1],[10,1],[11,1],[2,1],[5,1],[7,1],[8,1]],        # D_5(a_1)+A_2
       [[2,1],[4,1],[5,1],[6,1],[7,1],[1,1]],                # A_5A_1
       [[1,1],[3,1],[4,1],[2,1],[6,1],[7,1],[8,1]],          # A_4A_3
       [[1,1],[3,1],[4,1],[2,1],[5,1]],                      # D_5
#       [[17,1],[1,1],[19,1],[6,1],[18,1],[4,1]],             # E_6(a_3)
       [[4,1],[6,1],[9,1],[11,1],[18,1],[19,1]],             # E_6(a_3)
#       [[2,1],[3,1],[4,1],[5,1],[7,1],[8,1],[63,1]],        # (D_4+A_2)^{(2)}     char 2
       [[3,1],[4,1],[2,1],[5,1],[7,1],[8,1]],                # D_4A_2
       [[5,1],[6,1],[7,1],[8,1],[1,1],[3,1],[2,1]],          # A_4A_2A_1
#       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1],[7,1]],        # D_5(a_1)A_1
       [[1,1],[10,1],[11,1],[2,1],[5,1],[7,1]],              # D_5(a_1)+A_1
       [[1,1],[3,1],[4,1],[5,1],[6,1]],                      # A_5
       [[1,1],[3,1],[4,1],[2,1],[6,1],[7,1]],                # A_4A_2
       [[1,1],[3,1],[4,1],[2,1],[6,1],[8,1]],                # A_4A_1^2
#       [[1,1],[3,1],[10,1],[12,1],[2,1],[5,1]],              # D_5(a_1)
       [[1,1],[10,1],[11,1],[2,1],[5,1]],                    # D_5(a_1)

       [[1,1],[3,1],[4,1],[6,1],[7,1],[8,1]],          # A_3^2
       [[1,1],[3,1],[4,1],[2,1],[6,1]],                # A_4A_1
#       [[3,1],[10,1],[12,1],[2,1],[5,1],[7,1],[8,1]],  # D_4(a_1)A_2
       [[2,1],[3,1],[10,1],[12,1],[7,1],[8,1]],        # D_4(a_1)+A_2
       [[3,1],[4,1],[2,1],[5,1],[7,1]],                # D_4A_1
       [[5,1],[6,1],[7,1],[1,1],[3,1],[2,1]],          # A_3A_2A_1
       [[1,1],[3,1],[4,1],[2,1]],                      # A_4
#       [[1,1],[3,1],[4,1],[6,1],[7,1],[44,1]],        # (A_3+A_2)^{(2)}     char 2
       [[3,1],[4,1],[2,1],[6,1],[7,1]],                # A_3A_2
#       [[3,1],[10,1],[12,1],[2,1],[5,1],[7,1]],        # D_4(a_1)A_1
       [[2,1],[3,1],[10,1],[12,1],[7,1]],              # D_4(a_1)+A_1
       [[5,1],[6,1],[7,1],[2,1],[3,1]],                # A_3A_1^2
       [[1,1],[3,1],[5,1],[6,1],[2,1],[8,1]],          # A_2^2A_1^2
       [[3,1],[4,1],[2,1],[5,1]],                      # D_4
#       [[3,1],[10,1],[12,1],[2,1],[5,1]],              # D_4(a_1)
       [[2,1],[3,1],[10,1],[12,1]],                    # D_4(a_1)
       [[1,1],[3,1],[4,1],[6,1]],                      # A_3A_1
       [[1,1],[3,1],[5,1],[6,1],[2,1]],                # A_2^2A_1
       [[1,1],[3,1],[5,1],[6,1]],                      # A_2^2
       [[1,1],[3,1],[2,1],[5,1],[7,1]],                # A_2A_1^3
       [[1,1],[3,1],[4,1]],                            # A_3
       [[2,1],[4,1],[1,1],[6,1]],                      # A_2A_1^2
       [[1,1],[3,1],[2,1]],                            # A_2A_1
       [[2,1],[3,1],[5,1],[7,1]],                      # A_1^4
       [[1,1],[3,1]],                                  # A_2
       [[1,1],[4,1],[6,1]],                            # A_1^3
       [[1,1],[6,1]],                                  # A_1^2
       [[1,1]]                                         # A_1
       ];

I:=[[],
    [4],           #??? E_8(a_1)
    [4,6],         #??? E_8(a_2)
    [2,3,5],       # E_8(a_3)
    [2,3,5,7],     # E_8(a_4)
    [],            # E_7
    [2,3,5,6],     # E_8(b_4)
    [2,3,5,6,8],   # E_8(a_5)
    [4],           # E_7(a_1)
    [1,2,3,5,6],   # E_8(b_5)
    [],            # D_7
    [1,2,3,5,6,8], # E_8(a_6)
    [4,6],         #??? E_7(a_2)
    [],            # E_6A_1

#    [4],             # (D_7(a_1))^{(2)}     char 2
    [4],             # D_7(a_1)
    [1,2,3,5,6,7],   # E_8(b_6)
    [2,3,5],         # E_7(a_3)
    [4],             # E_6(a_1)A_1
#    [],              # (A_7)^{(3)}     char 3
    [],              # A_7
    [4,5],           #??? D_7(a_2)
    [],              # E_6
    [],              # D_6
#    [],              # (D_5+A_2)^{(2)}     char 2
    [],              # D_5A_2
    [4],             # E_6(a_1)
    [2,3,5,6],       # E_7(a_4)
    [],              # A_6A_1
    [4],             # D_6(a_1)
    [],              # A_6
    [1,2,3,4,6,7,8], # E_8(a_7)

    [],          # D_5A_1
    [1,2,3,5,6], # E_7(a_5)
    [2,3,5],     # E_6(a_3)A_1
    [4,5],       #??? D_6(a_2)
    [4],         # D_5(a_1)A_2
    [],          # A_5A_1
    [],          # A_4A_3
    [],          # D_5
    [2,3,5],     # E_6(a_3)
#    [],          # (D_4+A_2)^{(2)}     char 2
    [],          # D_4A_2
    [],          # A_4A_2A_1
    [4],         # D_5(a_1)A_1
    [],          # A_5
    [],          # A_4A_2
    [],          # A_4A_1
    [4],         # D_5(a_1) check this also

    [],  # A_3^2
    [],  # A_4A_1
    [4], # D_4(a_1)A_2
    [],  # D_4A_1
    [],  # A_3A_2A_1
    [],  # A_4
#    [],  # (A_3+A_2)^{(2)}     char 2
    [],  # A_3A_2
    [4], # D_4(a_1)A_1
    [],  # A_3A_1^2
    [],  # A_2^2A_1^2
    [],  # D_4
    [4], # D_4(a_1)
    [],  # A_3A_1
    [],  # A_2^2A_1
    [],  # A_2^2
    [],  # A_2A_1^3
    [],  # A_3
    [],  # A_2A_1^2
    [],  # A_2A_1
    [],  # A_1^4
    [],  # A_2
    [],  # A_1^3
    [],  # A_1^2
    []   # A_1
    ];

#
# The parabolic where the representative is distinguished
#
II:=[[1,2,3,4,5,6,7,8], # E_8
     [1,2,3,4,5,6,7,8], # E_8(a_1)
     [1,2,3,4,5,6,7,8], # E_8(a_2)
     [1,2,3,4,5,6,7,8], # E_8(a_3)
     [1,2,3,4,5,6,7,8], # E_8(a_4)
     [1,2,3,4,5,6,7],   # E_7
     [1,2,3,4,5,6,7,8], # E_8(b_4)
     [1,2,3,4,5,6,7,8], # E_8(a_5)
     [1,2,3,4,5,6,7],   # E_7(a_1)
     [1,2,3,4,5,6,7,8], # E_8(b_5)
     [2,3,4,5,6,7,8],   # D_7
     [1,2,3,4,5,6,7,8], # E_8(a_6)
     [1,2,3,4,5,6,7],   # E_7(a_2)
     [1,2,3,4,5,6,8],   # E_6A_1

#     [2,3,4,5,6,7,8],   # D_7(a_1)^{(2)}     char 2
     [2,3,4,5,6,7,8],   # D_7(a_1)
     [1,2,3,4,5,6,7,8], # E_8(b_6)
     [1,2,3,4,5,6,7],   # E_7(a_3)
     [1,2,3,4,5,6,8],   # E_6(a_1)A_1
#     [1,3,4,5,6,7,8],   # A_7^{(3)}     char 3
     [1,3,4,5,6,7,8],   # A_7
     [2,3,4,5,6,7,8],   # D_7(a_2)
     [1,2,3,4,5,6],     # E_6
     [2,3,4,5,6,7],     # D_6
#     [1,2,3,4,5,7,8],   # D_5A_2^{(2)}     char 2
     [1,2,3,4,5,7,8],   # D_5A_2
     [1,2,3,4,5,6],     # E_6(a_1)
     [1,2,3,4,5,6,7],   # E_7(a_4)
     [1,2,4,5,6,7,8],   # A_6A_1
     [2,3,4,5,6,7],     # D_6(a_1)
     [2,4,5,6,7,8],     # A_6
     [1,2,3,4,5,6,7,8], # E_8(a_7)

     [1,2,3,4,5,7],   # D_5A_1
     [1,2,3,4,5,6,7], # E_7(a_5)
     [1,2,3,4,5,6,8], # E_6(a_3)A_1
     [2,3,4,5,6,7],   # D_6(a_2)
     [1,2,3,4,5,7,8], # D_5(a_1)A_2
     [1,2,4,5,6,7],   # A_5A_1
     [1,2,3,4,6,7,8], # A_4A_3
     [1,2,3,4,5],     # D_5
     [1,2,3,4,5,6],   # E_6(a_3)
#     [2,3,4,5,7,8],   # D_4A_2^{(2)}     char 2
     [2,3,4,5,7,8],   # D_4A_2
     [1,2,3,5,6,7,8], # A_4A_2A_1
     [1,2,3,4,5,7],   # D_5(a_1)A_1
     [1,3,4,5,6],     # A_5
     [1,2,3,4,6,7],   # A_4A_2
     [1,2,3,4,6,8],   # A_4A_1^2
     [1,2,3,4,5],     # D_5(a_1)

     [1,3,4,6,7,8], # A_3^2
     [1,2,3,4,6],   # A_4A_1
     [2,3,4,5,7,8], # D_4(a_1)A_2
     [2,3,4,5,7],   # D_4A_1
     [1,2,3,5,6,7], # A_3A_2A_1
     [1,2,3,4],     # A_4
#     [2,3,4,6,7],   # A_3A_2^{(2)}     char 2
     [2,3,4,6,7],   # A_3A_2
     [2,3,4,5,7],   # D_4(a_1)A_1
     [2,3,5,6,7],   # A_3A_1^2
     [1,2,3,5,6,8], # A_2^2A_1^2
     [2,3,4,5],     # D_4
     [2,3,4,5],     # D_4(a_1)
     [1,3,4,6],     # A_3A_1
     [1,2,3,5,6],   # A_2^2A_1
     [1,3,5,6],     # A_2^2
     [1,2,3,5,7],   # A_2A_1^3
     [1,3,4],       # A_3
     [1,2,4,6],     # A_2A_1^2
     [1,2,3],       # A_2A_1
     [2,3,5,7],     # A_1^2
     [1,3],         # A_2
     [1,4,6],       # A_1^3
     [1,6],         # A_1^2
     [1]            # A_1
     ];

labels:=["E_8",
         "E_8(a_1)",
         "E_8(a_2)",
         "E_8(a_3)",
         "E_8(a_4)",
         "E_7",
         "E_8(b_4)",
         "E_8(a_5)",
         "E_7(a_1)",
         "E_8(b_5)",
         "D_7",
         "E_8(a_6)",
         "E_7(a_2)",
         "E_6A_1",

#         "D_7(a_1)^{(2)}",
         "D_7(a_1)",
         "E_8(b_6)",
         "E_7(a_3)",
         "E_6(a_1)A_1",
#         "A_7^{(3)}",
         "A_7",
         "D_7(a_2)",
         "E_6",
         "D_6",
#         "D_5A_2^{(2)}",
         "D_5A_2",
         "E_6(a_1)",
         "E_7(a_4)",
         "A_6A_1",
         "D_6(a_1)",
         "A_6",
         "E_8(a_7)",

         "D_5A_1",
         "E_7(a_5)",
         "E_6(a_3)A_1",
         "D_6(a_2)",
         "D_5(a_1)A_2",
         "A_5A_2",
         "A_4A_3",
         "D_5",
         "E_6(a_3)",
#         "D_4A_2^{(2)}",
         "D_4A_2",
         "A_4A_2A_1",
         "D_5(a_1)A_1",
         "A_5",
         "A_4A_2",
         "A_4A_1^2",
         "D_5(a_1)",

         "A_3^2",
         "A_4A_1",
         "D_4(a_1)A_2",
         "D_4A_1",
         "A_3A_2A_1",
         "A_4",
#         "A_3A_2^{(2)}",
         "A_3A_2",
         "D_4(a_1)A_1",
         "A_3A_1^2",
         "A_2^2A_1^2",
         "D_4",
         "D_4(a_1)",
         "A_3A_1",
         "A_2^2A_1",
         "A_2^2",
         "A_2A_1^3",
         "A_3",
         "A_2A_1^2",
         "A_2A_1",
         "A_1^4",
         "A_2",
         "A_1^3",
         "A_1^2",
         "A_1"
         ];

tori:=[
    [2,2,2,2,2,2,2,2],   # E_8
    [2,2,2,0,2,2,2,2],   # E_8(a_1)
    [2,2,2,0,2,0,2,2],   # E_8(a_2)
    [2,0,0,2,0,2,2,2],   # E_8(a_3)
    [2,0,0,2,0,2,0,2],   # E_8(a_4)
    [2,2,2,2,2,2,2,-27], # E_7
    [2,0,0,2,0,0,2,2],   # E_8(b_4)
    [2,0,0,2,0,0,2,0],   # E_8(a_5)
    [2,2,2,0,2,2,2,-21], # E_7(a_1)
    [0,0,0,2,0,0,2,2],   # E_8(b_5)
    [-21,2,2,2,2,2,2,2], # D_7
    [0,0,0,2,0,0,2,0],   # E_8(a_6)
    [2,2,2,0,2,0,2,-17], # E_7(a_2)
    [2,2,2,2,2,2,-17,2], # E_6A_1

#    [-16,2,0,2,2,2,2,2],   # D_7(a_1)^{(2)}     char 2
    [-16,2,0,2,2,2,2,2],   # D_7(a_1)
    [0,0,0,2,0,0,0,2],     # E_8(b_6)
    [2,0,0,2,0,2,2,-15],   # E_7(a_3)
    [2,2,2,0,2,2,-13,2],   # E_6(a_1)A_1
#    [2,-15,2,2,2,2,2,2],   # A_7^{(3)}     char 4
    [2,-15,2,2,2,2,2,2],   # A_7
    [-13,2,2,0,2,0,2,2],   # D_7(a_2)
    [2,2,2,2,2,2,-16,0],   # E_6
    [-15,2,2,2,2,2,2,-10], # D_6
#    [2,2,2,2,2,-12,2,2],   # D_5A_2^{(2)}     char 2
    [2,2,2,2,2,-12,2,2],   # D_5A_2
    [2,2,2,0,2,2,-12,0],   # E_6(a_1)
    [2,0,0,2,0,0,2,-11],   # E_7(a_4)
    [2,2,-11,2,2,2,2,2],   # A_6(a_1)
    [-11,2,2,0,2,2,2,-8],  # D_6(a_1)
    [0,2,-10,2,2,2,2,2],   # A_6
    [0,0,0,0,2,0,0,0],     # E_8(a_7)

    [2,2,2,2,2,-11,2,-1], # D_5A_1
    [0,0,0,2,0,0,2,-9],   # E_7(a_5)
    [2,0,0,2,0,2,-9,2],   # E_6(a_3)A_1
    [-9,2,2,0,2,0,2,-6],  # D_6(a_2)
    [2,2,2,0,2,-9,2,2],   # D_5(a_1)A_2
    [2,2,-9,2,2,2,2,-5],  # A_5A_1
    [2,2,2,2,-9,2,2,2],   # A_4A_3
    [2,2,2,2,2,-10,0,0],  # D_5
    [2,0,0,2,0,2,-8,0],   # E_6(a_3)
#    [-6,2,2,2,2,-8,2,2],  # D_4A_2^{(2)}     char 2
    [-6,2,2,2,2,-8,2,2],  # D_4A_2
    [2,2,2,-7,2,2,2,2],   # A_4A_2A_1
    [2,2,2,0,2,-8,2,-1],  # D_5(a_1)A_1
    [2,-9,2,2,2,2,-5,0],  # A_5
    [2,2,2,2,-8,2,2,-2],  # A_4A_2
    [2,2,2,2,-7,2,-2,2],  # A_4A_1^2
    [2,2,2,0,2,7,0,0],    # D_5(a_1)

    [2,-3,2,2,-6,2,2,2],   # A_3^2
    [2,2,2,2,-7,2,-1,0],   # A_4A_1
    [-4,2,2,0,2,-6,2,2],   # D_4(a_1)A_2
    [-6,2,2,2,2,-7,2,-1],  # D_4A_1
    [2,2,2,-6,2,2,2,-3],   # A_3A_2A_1
    [2,2,2,2,-6,0,0,0],    # A_4
#    [-3,2,2,2,-6,2,2,-2],  # A_3A_2^{(2)}     char 2
    [-3,2,2,2,-6,2,2,-2],  # A_3A_2
    [-4,2,2,0,2,-5,2,-1],  # D_4(a_1)A_1
    [-1,2,2,-5,2,2,2,-3],  # A_3A_1^2
    [2,2,2,-5,2,2,-3,2],   # A_2^2A_1^2
    [-6,2,2,2,2,-6,0,0],   # D_4
    [-4,2,2,0,2,-4,0,0],   # D_4(a_1)
    [2,-3,2,2,-4,2,-1,0],  # A_3A_1
    [2,2,2,-5,2,2,-2,0],   # A_2^2A_1
    [2,0,2,-4,2,2,-2,0],   # A_2^2
    [2,2,2,-4,2,-2,2,-1],  # A_2A_1^3
    [2,-3,2,2,-3,0,0,0],   # A_3
    [2,2,-3,2,-3,2,-1,0],  # A_2A_1^2
    [2,2,2,-3,0,0,0,0],    # A_2A_1
    [-1,2,2,-3,2,-2,2,-1], # A_1^4
    [2,2,0,-2,0,0,0,0],    # A_2
    [2,-1,-2,2,-2,2,-1,0], # A_1^3
    [2,0,-1,0,-1,2,-1,0],  # A_1^2
    [2,0,-1,0,0,0,0,0]     # A_1
    ];

ords:=[];

return [reps,tori,labels,I,II,ords];


text:
CaesarCipher := function(s, n)
	local r, c, i, lower, upper;
	lower := "abcdefghijklmnopqrstuvwxyz";
	upper := "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
	r := "";
	for c in s do
		i := Position(lower, c);
		if i <> fail then
			Add(r, lower[RemInt(i + n - 1, 26) + 1]);
		else
			i := Position(upper, c);
			if i <> fail then
				Add(r, upper[RemInt(i + n - 1, 26) + 1]);
			else
				Add(r, c);
			fi;
		fi;
	od;
	return r;
end;

CaesarCipher("IBM", 25);
# "HAL"

CaesarCipher("Vgg cphvi wzdibn vmz wjmi amzz viy zlpvg di ydbidot viy mdbcon.", 5);
# "All human beings are born free and equal in dignity and rights."

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
%Define
    $kw_lexer_class /.NoKWLexer./
    $_IDENTIFIER /.0./
%End
%Headers
    --
    -- Additional methods for the action class not provided in the template
    --
    /.
         public class NoKWLexer
        {
            public int[] getKeywordKinds() { return null; }

            public int lexer(int curtok, int lasttok) { return 0; }

            public void setInputChars(char[] inputChars) { }

             int getKind(int c) { return 0; }

            public NoKWLexer(char[] inputChars, int identifierKind) { }
        }
    ./
%End

%Import
    LexerBasicMapF.gi
%End

ext: gi

max_stars_repo_name: kuafuwang/LPG-VScode

max_issues_repo_name: kuafuwang/LPG-VScode

max_forks_repo_name: kuafuwang/LPG-VScode


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# idea : Restrict by dims(@(2)) <= codelet size
#
RewriteRules(RulesSums, rec(
    CommutePrmScatInplace := ARule(Compose, [@(1, [Prm,Scat]), [Inplace, @(2)]],
	e -> [ Inplace(@(1).val * @(2).val * @(1).val.transpose()), Copy(@(1).val) ]),

    CommuteInplaceGath := ARule(Compose, [[Inplace, @(2)], @(1, Gath)],
	e -> [ Copy(@(1).val), Inplace(@(1).val.transpose() * @(2).val * @(1).val) ])
));

Class(RulesInplace, RuleSet);
RewriteRules(RulesInplace, rec(
    CommutePrmScatInplace := RulesSums.rules.CommutePrmScatInplace,
    CommuteInplaceGath := RulesSums.rules.CommuteInplaceGath,

    InplaceSUM := Rule([Inplace, @(1,SUM)], e -> SUM(List(@(1).val.children(), Inplace)))
));


text:
#############################################################################
##
#W  ncgraphs.gi                NCGraph Package                    John Gouwar
##
##  Installation file for functions of the NCGraph package.
##
##

#############################################################################
##
#F  RemoveCenter( <G> ) . . . . . . . .  removes the center of a group
##
InstallGlobalFunction( RemoveCenter, function( G )
  local center, removed_center, g;
  center := Centre(G);
  removed_center := [];
  for g in G do
    if g in center then
      continue;
    else
      Add(removed_center, g);
    fi;
  od;
  return removed_center;
end );

#############################################################################
##
#F  NonCommutingGraph(<G>) . . . . . . .  generate non-commuting graph of <G>
##
InstallGlobalFunction( NonCommutingGraph, function( G )
  local non_com_relation, id;
    id := One(G);
  non_com_relation := function(x, y) return Comm(x, y) <> id; end;
  if IsAbelian(G) then
    return Graph(Group(id), [], OnPoints, non_com_relation);
  fi;
  return Graph(Group(id), RemoveCenter(G), OnPoints, non_com_relation);
end );

#############################################################################
##
#F  AdjacencyMatrix(<G>) . . . . . . .  generate an adjaceny matrix of <G>
##
InstallGlobalFunction( AdjacencyMatrix, function( G )
  local adj_mat, order, adjacencies, i, j;
  order := G.order;
  if order = 0 then
    return [];
  fi;
  adjacencies := G.adjacencies;
  adj_mat := NullMat(order, order);
  for i in [1..order] do
    for j in adjacencies[i] do
      adj_mat[i][j] := 1;
    od;
  od;
  return adj_mat;
end );

#E  ncgraphs.gi  . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: GouwarPower/NCGraph

max_issues_repo_name: GouwarPower/NCGraph

max_forks_repo_name: GouwarPower/NCGraph


text:
local reps,tori,labels,I,II,ords,mizuno,blocks;

mizuno:=[[ 1, 2, 3, 4, 5, 6, fail,
           7, 8, 9, 10,
           11, 6, 12, 13,
           15, 14, 16, 11,
           17, 18, 19, 21,
           20, 16, 22, 23,
           24, 25, 21, 20,
           26, 27, 23, 28,
           25, 29, 30, 27,
           31, 28, 32, 33,
           30, 31, 34, 32,
           33, fail, 35, 34,
           fail, 36, 35, fail,
           36, fail, fail, fail,
           fail, fail, fail, fail ],
          []];

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1]],          # E_6
       [[2,1],[4,1],[5,1],[6,1],[8,1],[16,1]],         # E_6(a_1)
       [[1,1],[2,1],[10,1],[11,1],[12,1]],             # D_5
       [[8,1],[9,1],[12,1],[10,1],[11,1],[16,1]],      # A_5+A_1
       [[8,1],[9,1],[16,1],[12,1],[18,1]],             # D_5(a_1)
       [[1,1],[15,1],[16,1],[17,1],[6,1]],             # A_5
       [[14,1],[15,1],[16,1],[17,1],[12,1]],           # A_4+A_1
       [[2,1],[14,1],[16,1],[18,1]],                   # D_4
       [[14,1],[15,1],[17,1],[12,1]],                  # A_4
       [[14,1],[16,1],[22,1],[24,1]],                  # D_4(a_1)
       [[14,1],[22,1],[23,1],[24,1]],                  # A_3+A_1
       [[20,1],[21,1],[28,1],[23,1],[24,1]],           # 2A_2+A_1
       [[14,1],[22,1],[24,1]],                         # A_3
       [[26,1],[27,1],[28,1],[29,1]],                  # A_2+2A_1
       [[20,1],[21,1],[23,1],[24,1]],                  # 2A_2
       [[26,1],[27,1],[35,1]],                         # A_2+A_1
       [[26,1],[35,1]],                                # A_2
       [[37,1],[38,1],[40,1]],                         # 3A_1
       [[42,1],[43,1]],                                # 2A_1
       [[53,1]]                                        # A_1
       ];

labels:=["E_6",
         "E_6(a_1)",
         "D_5",
         "A_5+A_1",
         "D_5(a_1)",
         "A_5",
         "A_4+A_1",
         "D_4",
         "A_4",
         "D_4(a_1)",
         "A_3+A_1",
         "2A_2+A_1",
         "A_3",
         "A_2+2A_1",
         "2A_2",
         "A_2+A_1",
         "A_2",
         "3A_1",
         "2A_1",
         "A_1",
         ];

I:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];

II:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];

tori:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];

#blocks:=[];
blocks:=[[[16, 4], [8, 1], [6, 1]], [[16, 2], [11, 2], [8, 2], [4, 2]], [[8, 9], [6, 1]], [[8, 8], [4, 2], [3, 2]], [[8, 5], [6, 1], [4, 8]], [[8, 8], [2, 6], [1, 2]], [[8, 4], [6, 2], [5, 2], [4, 4], [2, 4]], [[8, 1], [6, 9], [2, 8]], [[8, 2], [7, 4], [5, 2], [4, 2], [3, 4], [1, 4]], [[4, 18], [3, 2]], [[4, 16], [2, 6], [1, 2]], [[4, 16], [2, 6], [1, 2]], [[4, 16], [2, 4], [1, 6]], [[4, 10], [3, 2], [2, 16]], [[4, 16], [1, 14]], [[4, 8], [3, 6], [2, 10], [1, 8]], [[4, 2], [3, 18], [1, 16]], [[2, 38], [1, 2]], [[2, 32], [1, 14]], [[2, 22], [1, 34]]];

#ords:=[];
ords:=[ 16, 16, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2 ];

return [reps,tori,labels,I,II,ords,blocks,mizuno];


text:
Import(paradigms.smp);

Class(fCubeTranspose, PermClass, rec(
    domain := self >> self.params[1],
    range  := self >> self.params[1],
    def    := (N, stage, cube) -> rec(),
    lambda := self >> let(i := Ind(self.params[1]), errExp(TInt)),
    transpose := self >> self,
    isIdentity := True
));

Class(MPIPrm, Prm);

Class(MPIRCPrm, Prm,
    rec(
        dims := self >> 2*Replicate(2, self.func.params[1])
));

Declare(MPITensor);

Class(MPISum, SMPSum);

Class(MPIGath, Gath);

Class(MPIScat, Scat);

Class(MPITensor, BaseMat, SumsBase, rec(
    dims := self >> self.L.dims() * self.P,
    isReal := self >> true,
    #-----------------------------------------------------------------------
    rChildren := self >> [self.L, self.P],
    rSetChild := rSetChildFields("L", "P"),
    #-----------------------------------------------------------------------
    new := (self, L, P) >> SPL(WithBases(self,
        rec(L   := L,
            P   := P,
            dimensions := L.dims()*P)
    )),

    #-----------------------------------------------------------------------
    transpose := self >> MPITensor(self.L.transpose(), self.P),
    #-----------------------------------------------------------------------

    print := (self,i,is) >> Print(self.name, "(",
        self.L.print(i+is,is), ", ", self.P, ")"),
    #-----------------------------------------------------------------------
    toAMat := self >> Tensor(I(self.P), self.L).toAMat(),
    #-----------------------------------------------------------------------
    isPermutation := False
));

Class(MPISumsgenMixin, rec(
    MPITensor := (self, o, opts) >> let(p:= o.P, i:= Ind(p),
        MPISum(p, i, p,
        MPIScat(fTensor(fBase(i), fId(Rows(o.L)))) * o.L * MPIGath(fTensor(fBase(i), fId(Cols(o.L))))
    ))
));

max_stars_repo_name: spiral-software/spiral-package-mpi

max_issues_repo_name: spiral-software/spiral-package-mpi

max_forks_repo_name: spiral-software/spiral-package-mpi


text:
1 < 2;
# true

2 < 1;
# false

# GAP has also the value fail, which cannot be used as a boolean but may be used i$

1 = fail;
# false

fail = fail;
# true


text:
# In reality we only need a list indexed by triples, so we
# will quite probably end up with a tree where the leaves have
# weights as labels.
InstallGlobalFunction(CyclicSubList,
function(l, pos, len)
    local r, i, j, llen;

    llen := Length(l);
    r := [];
    i := pos; j := 1;
    while j <= len do
        r[j] := l[i];
        i := i + 1; j := j + 1;
        if i > llen then
            i := 1;
        fi;
    od;

    return r;
end);

InstallGlobalFunction( EnterAllSubwords,
function(idx, word, value)
    local i, v, pos, ww;
    pos := [1..Length(word)];
    ww := List(word, x -> __ID(x));
    for i in pos do
        v := idx[ ww{ CyclicSubList(pos, i, 3) } ];
        if value < v then
            idx[ ww{ CyclicSubList(pos, i, 3) } ] := value;
        fi;
    od;
end);

#
InstallGlobalFunction( DigraphDijkstraST,
function(graph, s, t)
    local vertices, dist, prev, queue, u, v, alt;

    dist := [];
    prev := [];
    queue := BinaryHeap({x,y} -> x[1] < y[1]);

    for v in DigraphVertices(graph) do
        dist[v] := infinity;
        prev[v] := -1;
    od;

    dist[s] := 0;
    Push(queue, [0, s]);

    while not IsEmpty(queue) do
        u := Pop(queue);
        u := u[2];
        if u = t then
            return [ dist, prev ];
        fi;
        for v in OutNeighbours(graph)[u] do
            alt := dist[u] + DigraphEdgeLabel(graph, u, v);
            if alt < dist[v] then
                dist[v] := alt;
                prev[v] := u;
                Push(queue, [dist[v], v]);
            fi;
        od;
    od;
    return infinity;
end);

InstallGlobalFunction( DigraphDijkstraS,
function(graph, s)
    local vertices, dist, prev, queue, u, v, alt;

    dist := [];
    prev := [];
    queue := BinaryHeap({x,y} -> x[1] < y[1]);

    for v in DigraphVertices(graph) do
        dist[v] := infinity;
        prev[v] := -1;
    od;

    dist[s] := 0;
    Push(queue, [0, s]);

    while not IsEmpty(queue) do
        u := Pop(queue);
        u := u[2];
        for v in OutNeighbours(graph)[u] do
            alt := dist[u] + DigraphEdgeLabel(graph, u, v);
            if alt < dist[v] then
                dist[v] := alt;
                prev[v] := u;
                Push(queue, [dist[v], v]);
            fi;
        od;
    od;

    return [ dist, prev ];
end);

max_stars_repo_name: RussWoodroofe/walrus

max_issues_repo_name: RussWoodroofe/walrus

max_forks_repo_name: RussWoodroofe/walrus


text:
#! @Title Skeletal vector spaces
#! @Author Sebastian Gutsche and Sebastian Posur
#! @Date 01/08/2018
#! @Chapter HandsOn

#! @Section Objects and morphisms

LoadPackage( "CAP" );

LoadPackage( "MatricesForHomalg" );

#################################
##
## Technical declarations
##
#################################

## Objects
DeclareCategory( "IsQVectorSpace",
                 IsCapCategoryObject );
## Morphisms
DeclareCategory( "IsQVectorSpaceMorphism",
                 IsCapCategoryMorphism );

#################################
##
## Attributes
##
#################################

#! @Description
#! The argument is an object $V$ in the category of rational vector spaces.
#! The output is its dimension.
#! @Returns a non-negative integer
#! @Arguments V
DeclareAttribute( "Dimension",
                  IsQVectorSpace );

#! @Description
#! The argument is a morphism $\alpha$ in the category of rational vector spaces.
#! The output is its underlying matrix.
#! @Returns a homalg matrix
#! @Arguments alpha
DeclareAttribute( "UnderlyingMatrix",
                  IsQVectorSpaceMorphism );

#################################
##
## Operations
##
#################################

#! @Description
#! The argument is a non-negative integer $d$.
#! The output is a rational vector space of dimension $d$.
#! @Returns an object
#! @Arguments d
DeclareOperation( "QVectorSpace",
                  [ IsInt ] );

#! @Description
#! The first argument is a rational vector space $V$.
#! The second argument $A$ is either a homalg matrix defined over the rationals
#! or an input that can be used as the first argument in the HomalgMatrix constructor.
#! The third argument is a rational vector space $W$.
#! The output is vector space morphism from $V$ to $W$ defined by $A$.
#! @Returns an morphism in $\mathrm{Hom}(V,W)$
#! @Arguments V,A,W
DeclareOperation( "QVectorSpaceMorphism",
                  [ IsQVectorSpace, IsObject, IsQVectorSpace ] );

#################################
##
## Creation of category
##
#################################

BindGlobal( "vecspaces", CreateCapCategory( "QVectorSpaces" ) );

AddObjectRepresentation( vecspaces, IsQVectorSpace );

AddMorphismRepresentation( vecspaces, IsQVectorSpaceMorphism );

SetIsAbelianCategory( vecspaces, true );

#################################
##
## Creation of Q
##
#################################

BindGlobal( "VECTORSPACES_FIELD", HomalgFieldOfRationals( ) );

#################################
##
## Constructors for objects and morphisms
##
#################################

##
InstallMethod( QVectorSpace,
               [ IsInt ],

  function( dim )
    local space;

    if dim < 0 then

        Error( "the argument has to be a non-negative integer" );

    fi;

    space := rec( );

    ObjectifyObjectForCAPWithAttributes( space, vecspaces,
                                         Dimension, dim
    );

    return space;

end );

##
InstallMethod( QVectorSpaceMorphism,
                  [ IsQVectorSpace, IsObject, IsQVectorSpace ],

  function( source, matrix, range )
    local morphism;

    if not IsHomalgMatrix( matrix ) then

      matrix := HomalgMatrix( matrix, Dimension( source ), Dimension( range ), VECTORSPACES_FIELD );

    fi;

    morphism := rec( );

    ObjectifyMorphismForCAPWithAttributes( morphism, vecspaces,
                                         Source, source,
                                         Range, range,
                                         UnderlyingMatrix, matrix
    );

    return morphism;

end );

#################################
##
## View
##
#################################

##
InstallMethod( ViewObj,
               [ IsQVectorSpace ],

  function( obj )

    Print( "<A rational vector space of dimension ", String( Dimension( obj ) ), ">" );

end );

##
InstallMethod( ViewObj,
               [ IsQVectorSpaceMorphism ],

  function( mor )

    Print( "A rational vector space homomorphism with matrix: \n" );

    Display( UnderlyingMatrix( mor ) );

end );

#################################
##
## Functions to be added to category
##
#################################

##
identity_morphism := function( obj )

    return QVectorSpaceMorphism( obj, HomalgIdentityMatrix( Dimension( obj ), VECTORSPACES_FIELD ), obj );

end;

AddIdentityMorphism( vecspaces, identity_morphism );

##
pre_compose := function( mor_left, mor_right )
    local composition;

    composition := UnderlyingMatrix( mor_left ) * UnderlyingMatrix( mor_right );

    return QVectorSpaceMorphism( Source( mor_left ), composition, Range( mor_right ) );

end;

AddPreCompose( vecspaces, pre_compose );

##
is_equal_for_objects := function( vecspace_1, vecspace_2 )

    return Dimension( vecspace_1 ) = Dimension( vecspace_2 );

end;

AddIsEqualForObjects( vecspaces, is_equal_for_objects );

##
is_equal_for_morphisms := function( a, b )

    return UnderlyingMatrix( a ) = UnderlyingMatrix( b );

end;

AddIsEqualForMorphisms( vecspaces, is_equal_for_morphisms );

##
kernel_emb := function( morphism )
    local syzygies, kernel_obj;

    syzygies := SyzygiesOfRows( UnderlyingMatrix( morphism ) );

    kernel_obj := QVectorSpace( NrRows( syzygies ) );

    return QVectorSpaceMorphism( kernel_obj, syzygies, Source( morphism ) );

end;

AddKernelEmbedding( vecspaces, kernel_emb );

##
lift := function( alpha, beta )
    local solution;

    solution := RightDivide( UnderlyingMatrix( alpha ), UnderlyingMatrix( beta ) );

    if solution = fail then

        return fail;

    fi;

    return QVectorSpaceMorphism( Source( alpha ),
           solution,
           Source( beta ) );

end;

AddLift( vecspaces, lift );

##
cokernel_proj := function( morphism )
    local syzygies, cokernel_obj;

    syzygies := SyzygiesOfColumns( UnderlyingMatrix( morphism ) );

    cokernel_obj := QVectorSpace( NrColumns( syzygies ) );

    return QVectorSpaceMorphism( Range( morphism ),
              syzygies, cokernel_obj );

end;

AddCokernelProjection( vecspaces, cokernel_proj );

##
colift := function( alpha, beta )
    local solution;

    solution := LeftDivide( UnderlyingMatrix( alpha ), UnderlyingMatrix( beta ) );

    if solution = fail then

        return fail;

    fi;

    return QVectorSpaceMorphism( Range( alpha ),
              solution,
              Range( beta ) );

end;

AddColift( vecspaces, colift );

##
zero_object := function( )

    return QVectorSpace( 0 );

end;

AddZeroObject( vecspaces, zero_object );

##
universal_morphism_into_zero_object := function( source )

    return QVectorSpaceMorphism( source,
              HomalgZeroMatrix( Dimension( source ), 0, VECTORSPACES_FIELD ),
              QVectorSpace( 0 ) );

end;

AddUniversalMorphismIntoZeroObject( vecspaces, universal_morphism_into_zero_object );

##
universal_morphism_into_zero_object_with_given_zero_object := function( source, terminal_object )

    return QVectorSpaceMorphism( source,
              HomalgZeroMatrix( Dimension( source ), 0, VECTORSPACES_FIELD ),
              terminal_object );

end;

AddUniversalMorphismIntoZeroObjectWithGivenZeroObject( vecspaces, universal_morphism_into_zero_object_with_given_zero_object );

##
universal_morphism_from_zero_object := function( sink )

    return QVectorSpaceMorphism( QVectorSpace( 0 ),
              HomalgZeroMatrix( 0, Dimension( sink ), VECTORSPACES_FIELD ),
              sink );

end;

AddUniversalMorphismFromZeroObject( vecspaces, universal_morphism_from_zero_object );

##
universal_morphism_from_zero_object_with_given_zero_object := function( sink, initial_object )

    return QVectorSpaceMorphism( initial_object,
                   HomalgZeroMatrix( 0, Dimension( sink ), VECTORSPACES_FIELD ),
                   sink );

end;

AddUniversalMorphismFromZeroObjectWithGivenZeroObject( vecspaces, universal_morphism_from_zero_object_with_given_zero_object );

##
addition_for_morphisms := function( a, b )

    return QVectorSpaceMorphism( Source( a ),
              UnderlyingMatrix( a ) + UnderlyingMatrix( b ),
              Range( a ) );

end;

AddAdditionForMorphisms( vecspaces, addition_for_morphisms );

##
additive_inverse_for_morphisms := function( a )

    return QVectorSpaceMorphism( Source( a ),
              - UnderlyingMatrix( a ),
              Range( a ) );

end;

AddAdditiveInverseForMorphisms( vecspaces, additive_inverse_for_morphisms );

##
direct_sum := function( object_product_list )
    local dim;

    dim := Sum( List( object_product_list, c -> Dimension( c ) ) );

    return QVectorSpace( dim );

end;

AddDirectSum( vecspaces, direct_sum );

##
injection_of_cofactor_of_direct_sum := function( object_product_list, injection_number )
    local components, dim, dim_pre, dim_post, dim_cofactor, coproduct, number_of_objects, injection_of_cofactor;

    components := object_product_list;

    number_of_objects := Length( components );

    dim := Sum( components, c -> Dimension( c ) );

    dim_pre := Sum( components{ [ 1 .. injection_number - 1 ] }, c -> Dimension( c ) );

    dim_post := Sum( components{ [ injection_number + 1 .. number_of_objects ] }, c -> Dimension( c ) );

    dim_cofactor := Dimension( object_product_list[ injection_number ] );

    coproduct := QVectorSpace( dim );

    injection_of_cofactor := HomalgZeroMatrix( dim_cofactor, dim_pre ,VECTORSPACES_FIELD );

    injection_of_cofactor := UnionOfColumns( injection_of_cofactor,
                                         HomalgIdentityMatrix( dim_cofactor, VECTORSPACES_FIELD ) );

    injection_of_cofactor := UnionOfColumns( injection_of_cofactor,
                                         HomalgZeroMatrix( dim_cofactor, dim_post, VECTORSPACES_FIELD ) );

    return QVectorSpaceMorphism( object_product_list[ injection_number ], injection_of_cofactor, coproduct );

end;

AddInjectionOfCofactorOfDirectSum( vecspaces, injection_of_cofactor_of_direct_sum );

##
universal_morphism_from_direct_sum := function( diagram, sink )
    local dim, coproduct, components, universal_morphism, morphism;

    components := sink;

    dim := Sum( components, c -> Dimension( Source( c ) ) );

    coproduct := QVectorSpace( dim );

    universal_morphism := UnderlyingMatrix( sink[1] );

    for morphism in components{ [ 2 .. Length( components ) ] } do

      universal_morphism := UnionOfRows( universal_morphism, UnderlyingMatrix( morphism ) );

    od;

    return QVectorSpaceMorphism( coproduct, universal_morphism, Range( sink[1] ) );

end;

AddUniversalMorphismFromDirectSum( vecspaces, universal_morphism_from_direct_sum );

##
projection_in_factor_of_direct_sum := function( object_product_list, projection_number )
    local components, dim, dim_pre, dim_post, dim_factor, direct_product, number_of_objects, projection_in_factor;

    components := object_product_list;

    number_of_objects := Length( components );

    dim := Sum( components, c -> Dimension( c ) );

    dim_pre := Sum( components{ [ 1 .. projection_number - 1 ] }, c -> Dimension( c ) );

    dim_post := Sum( components{ [ projection_number + 1 .. number_of_objects ] }, c -> Dimension( c ) );

    dim_factor := Dimension( object_product_list[ projection_number ] );

    direct_product := QVectorSpace( dim );

    projection_in_factor := HomalgZeroMatrix( dim_pre, dim_factor, VECTORSPACES_FIELD );

    projection_in_factor := UnionOfRows( projection_in_factor,
                                         HomalgIdentityMatrix( dim_factor, VECTORSPACES_FIELD ) );

    projection_in_factor := UnionOfRows( projection_in_factor,
                                         HomalgZeroMatrix( dim_post, dim_factor, VECTORSPACES_FIELD ) );

    return QVectorSpaceMorphism( direct_product, projection_in_factor, object_product_list[ projection_number ] );

end;

AddProjectionInFactorOfDirectSum( vecspaces, projection_in_factor_of_direct_sum );

##
universal_morphism_into_direct_sum := function( diagram, sink )
    local dim, direct_product, components, universal_morphism, morphism;

    components := sink;

    dim := Sum( components, c -> Dimension( Range( c ) ) );

    direct_product := QVectorSpace( dim );

    universal_morphism := UnderlyingMatrix( sink[1] );

    for morphism in components{ [ 2 .. Length( components ) ] } do

      universal_morphism := UnionOfColumns( universal_morphism, UnderlyingMatrix( morphism ) );

    od;

    return QVectorSpaceMorphism( Source( sink[1] ), universal_morphism, direct_product );

end;

AddUniversalMorphismIntoDirectSum( vecspaces, universal_morphism_into_direct_sum );

#################################
##
## Finalize category
##
#################################

Finalize( vecspaces );

#################################
##
## Test the basic operations
##
#################################

# Creating objects and morphisms

V := QVectorSpace( 2 );

CapCategory( V );

Dimension( V );

W := QVectorSpace( 3 );

alpha := QVectorSpaceMorphism( V, [ [ 1, 1, 1 ], [ -1, -1, -1 ] ], W );

CapCategory( alpha );

UnderlyingMatrix( alpha );

# Testing the KernelEmbedding

KernelEmbedding( alpha );

KernelObject( alpha );

# Computing an intersection

M1 := QVectorSpace( 2 );

M2 := QVectorSpace( 2 );

N := QVectorSpace( 3 );

iota1 := QVectorSpaceMorphism( M1, [ [ 1, 0, 0 ], [ 0, 1, 1 ] ], N );

IsMonomorphism( iota1 );

iota2 := QVectorSpaceMorphism( M2, [ [ 1, 1, 0 ], [ 0, 0, 1 ] ], N );

IsMonomorphism( iota2 );

pi1 := ProjectionInFactorOfDirectSum( [ M1, M2 ], 1 );

pi2 := ProjectionInFactorOfDirectSum( [ M1, M2 ], 2 );

lambda := PostCompose( iota1, pi1 );

phi := lambda - PostCompose( iota2, pi2 );

kappa := KernelEmbedding( phi );

PostCompose( lambda, kappa );

PreCompose( ProjectionInFactorOfFiberProduct( [ iota1, iota2 ], 1 ), iota1 );

#################################
##
## A function for computing homology
##
#################################

HomologyObject := function( alpha, beta )
  local iota, lambda;

  if not IsZero( PreCompose( alpha, beta ) ) then

      Error( "the composition of the given morphisms has to be zero" );

  fi;

  iota := ImageEmbedding( alpha );

  lambda := KernelLift( beta, iota );

  return CokernelObject( lambda );

end;

HomologyObject( alpha, CokernelProjection( alpha ) );

HomologyObject( KernelEmbedding( alpha ), alpha );

pi1 := ProjectionInFactorOfDirectSum( [ V, V, V ], 1 );

iota1 := InjectionOfCofactorOfDirectSum( [ V, V, V ], 1 );

pi2 := ProjectionInFactorOfDirectSum( [ V, V, V ], 2 );

iota2 := InjectionOfCofactorOfDirectSum( [ V, V, V ], 2 );

HomologyObject( PreCompose( pi1, iota1 ), PreCompose( pi2, iota2 ) );


text:
Class(RulesMPIPromoteNT, RuleSet);

RewriteRules(RulesMPIPromoteNT, rec(
    DAG_remove_mpi := Rule([@(1,TDAG), [@(2, TDAGNode), @(3),
            @(4).cond(e->IsList(e) and Length(e) = 1 and ForAll(Collect(e[1], var), v-> v.id in ["Y", "Yptr"])),
            @(5).cond(e->IsList(e) and Length(e) = 1 and ForAll(Collect(e[1], var), v-> v.id in ["X", "Xptr"])),
            ...],...],
        e-> @(2).val.params[1].withTags([AParMPI(rec(Y := @(4).val, X := @(5).val))])),

    TFunc_mergeTags := Rule([@(1, TFCall, e->e.hasTags()), @(2).cond(e->e.hasTags() and let(tg := e.getTags()[1], IsBound(tg.isParMPI) and tg.isParMPI)), ...],
        e-> ApplyFunc(TFCall, [@(2).val.setTags([])]::Drop(@(1).val.params,1)).withTags(
            [ApplyFunc(ObjId(@(1).val.getTags()[1]), @(1).val.getTags()[1].params::@(2).val.getTags()[1].params)]::Drop(@(1).val.getTags(),1))
    )
));

RewriteRules(RulesRC, rec(
    RC_MPISum := Rule([RC, @(1, MPISum)],
        e -> let(s:=@(1).val, MPISum(s.nthreads, s.tid, s.var, s.domain, RC(s.child(1))))),

    RC_MPIPrm := Rule([RC, @(1,MPIPrm)],
        e -> MPIRCPrm(@(1).val.func)),

    RC_MPIGath := Rule([RC, @(1, MPIGath)], e -> MPIGath(fTensor(@(1).val.func, fId(2)))),

    RC_MPIScat := Rule([RC, @(1, MPIScat)], e -> MPIScat(fTensor(@(1).val.func, fId(2))))

));

Class(LocalizeMPIRules, RuleSet);
RewriteRules(LocalizeMPIRules, rec(
    merge_GPS := ARule(Compose, [@(1, MPIGath), @(2, MPIRCPrm), @(3, MPIScat)],
        e->let(p := @(1).val.free()[1].range, f := @(2).val.func.params,
               [ MPIRCPrm(ApplyFunc(ObjId(@(2).val.func), [f[1]/p]::Drop(f, 1)))])),

    merge_GRCDS := ARule(Compose, [@(1, MPIGath), @(2, RCDiag), @(3, MPIScat)],
        e->let(p := @(1).val.free()[1].range, rcdf := @(2).val.element,
               [ RCDiag(FDataOfs(rcdf.var, rcdf.len/p, rcdf.ofs)) ])),

    drop_NoPull := Rule(@(1, NoPull), e->e.child(1)),

    drop_MPISum := Rule(@(1, MPISum), e->e.child(1)),

    drop_GathScat := ARule(Compose, [@(1, MPIGath), @(2, MPIScat)],
        e-> [  ])
));

Class(LocalizeMPIRulesCleanup, RuleSet);
RewriteRules(LocalizeMPIRulesCleanup, rec(
    drop_MPIGathScath := ARule(Compose, [@(1, [MPIGath, MPIScat])], e-> [])
));

FullLocalizeMPI := [MergedRuleSet(RulesSums, LocalizeMPIRules), LocalizeMPIRulesCleanup];


text:
f := function(n)
  return f(n+1);
end;

# Now loop until an error occurs
f(0);

# Error message :
#   Entering break read-eval-print loop ...
#   you can 'quit;' to quit to outer loop, or
#   you may 'return;' to continue

n;
# 4998

# quit "brk mode" and return to GAP
quit;


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details

FFTE_copyright :=
"/*\n"::
"\n"::
"     FFTE: A FAST FOURIER TRANSFORM PACKAGE\n"::
"\n"::
"     (C) COPYRIGHT SOFTWARE, 2000-2004, 2008-2014, ALL RIGHTS RESERVED\n"::
"                BY\n"::
"         DAISUKE TAKAHASHI\n"::
"         FACULTY OF ENGINEERING, INFORMATION AND SYSTEMS\n"::
"         UNIVERSITY OF TSUKUBA\n"::
"         1-1-1 TENNODAI, TSUKUBA, IBARAKI 305-8573, JAPAN\n"::
"         E-MAIL: daisuke@cs.tsukuba.ac.jp\n"::
"\n"::
"\n"::
"     WRITTEN BY DAISUKE TAKAHASHI\n\n"::
"     THIS KERNEL WAS GENERATED BY SPIRAL "::Version()::"\n"::
"*/\n\n";

max_stars_repo_name: spiral-software/spiral-package-ffte

max_issues_repo_name: spiral-software/spiral-package-ffte

max_forks_repo_name: spiral-software/spiral-package-ffte


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Declare(VTensor, RCVTensor, VTensorInd);

#   A x I_v
Class(VTensor, Tensor, rec(
    new := (self, L) >> SPL(WithBases(self, rec(
        _children := [L[1]],
        dimensions := When(IsBound(L[1].dims), L[1].dims(), L[1].dimensions) * L[2],
        vlen := L[2]))),
    # NOTE: vlen not exposed to rChildren
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.vlen),

    normalizedArithCost := (self) >> self.vlen * self.child(1).normalizedArithCost(),
    print := (self,i,is) >> Print(self.name, "(",
    self.child(1).print(i+is,is), ", ", self.vlen, ")"),
    toAMat := self >> Tensor(self.child(1), I(self.vlen)).toAMat(),
    sums := self >> Inherit(self, rec(_children := [self.child(1).sums()])),
    isPermutation := False,
    dims := self >> self.child(1).dims() * self.vlen,
    needInterleavedLeft := False,
    needInterleavedRight := False,
    transpose := self >> VTensor(self.child(1).transpose(), self.vlen),
    # for BlockSums
    isBlockTransitive := true,
    area := self >> self.child(1).area() * self.vlen,
    #NOTE: check if these two should be set to true or false
    cannotChangeDataFormat := False,
    totallyCannotChangeDataFormat := False
));

#   RC(A x I_v)
Class(RCVTensor, Tensor, rec(
    new := (self, L) >> SPL(WithBases(self, rec(
        _children := [L[1]],
        dimensions := 2 * L[1].dimensions * L[2],
        vlen := L[2]))),
    # NOTE: vlen not exposed to rChildren
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.vlen),
    print := (self,i,is) >> Print(self.name, "(",
    self.child(1).print(i+is,is), ", ", self.vlen, ")"),
    toAMat := self >> RC(Tensor(self.child(1), I(self.vlen)).toAMat()),
    sums := self >> Inherit(self, rec(_children := [self.child(1).sums()])),
    isPermutation := False,
    dims := self >> self.child(1).dims() * self.vlen * 2,
    needInterleavedLeft := True,
    needInterleavedRight := True,
    transpose := self >> RCVTensor(self.child(1).transpose(), self.vlen),
    # for BlockSums
    isBlockTransitive := true
));

Class(VTensorInd, Tensor, rec(
    new := (self, L) >> SPL(WithBases(self, rec(
        _children := [L[1]],
        dimensions := When(IsBound(L[1].dims), L[1].dims(), L[1].dimensions) * L[2].range,
        vlen := L[2]))),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.vlen),

    rChildren := self >> [self._children[1], self.vlen],

    rSetChild := meth ( self, n, what )
            if n=1 then
                self._children[1] := what;
            elif n=2 then
                self.vlen := what;
            else
                Error("VTensorInd only has 2 rChildren");
            fi;
        end,
#    from_rChildren := (self, rch) >> ObjId(self)(rch[1], rch[2]),

    print := (self,i,is) >> Print(self.name, "(",
    self.child(1).print(i+is,is), ", ", self.vlen, ")"),
    toAMat := self >> let(d := self.child(1).dims(), v := self.vlen.range, L(d[1]*v, d[1]).toAMat() * IDirSum(self.vlen, self.child(1)).toAMat() * L(d[2]*v, v).toAMat()),
    sums := self >> Inherit(self, rec(_children := [self.child(1).sums()])),
    isPermutation := False,
    dims := self >> self.child(1).dims() * self.vlen.range,
    needInterleavedLeft := False,
    needInterleavedRight := False,
    transpose := self >> VTensorInd(self.child(1).transpose(), self.vlen),
    # for BlockSums
    isBlockTransitive := true,
    area := self >> self.child(1).area() * self.vlen.range,
    free := meth(self)
        local fvar;
        fvar := self.child(1).free();
        SubtractSet(fvar, Set([self.vlen]));
        return fvar;
    end
));


text:
Load(graph);
Import(graph);
conf := SparseDefaultConf;
opts := conf.getOpts(V(1));
opts.vec := true;
i := Ind(3);
front := var.fresh_t("frontier",TArray(TSemiring_MinTimes(TInt), opts.symbol[1]));
fdos := FDataOfs(front, opts.symbol[1], i);
rv := RowVec(fdos);
spmv := IterVStack(i, opts.symbol[1], rv);
rt := RandomRuleTree(last, opts);
srt := SumsRuleTree(rt, opts);
cs := CodeSums(BFSKernel(PullBFSKernel(srt)), opts);
PrintCode("spmv", cs, opts);

# c -> MarkDefUse(c), c -> FixValueTypes(c), c -> Compile.declareVars(c),
#    (c, opts) -> HashConsts(c, opts)


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

_setN := function(o) o.params[2].nested := true; return o; end;

#F Filt(<n>, <taps>, <v>, <ds>) - FIR filter (linear convolution)
#F
#F Filt is a non-terminal for a linear convolution matrix optionallly
#F downsampled by a factor of <ds>.
#F
#F For ds=1, it is a (n x n+l-1) matrix,  where l is the number of filter
#F coefficients in <taps>.
#F
#F For ds<>1, it is a (n x ds*(n-1)+l) matrix.
#F
#F <n>    - number of outputs desired
#F <taps> -  list of filter taps or a polynomial transfer function
#F <v>    - valuation, optional, but used when <taps> is not a polynomial
#F <ds>   - downsampling factor
#F
#F Examples: Filt(5, [1, 2, 3, 4, 5]) corresponds to the matrix
#F             [ [1, 2, 3, 4, 5, 0, 0, 0, 0]
#F               [0, 1, 2, 3, 4, 5, 0, 0, 0]
#F               [0, 0, 1, 2, 3, 4, 5, 0, 0]
#F               [0, 0, 0, 1, 2, 3, 4, 5, 0]
#F               [0, 0, 0, 0, 1, 2, 3, 4, 5] ]
#F
#F           Filt(3, [1, 2, 3], -2, 2) corresponds to the matrix
#F             [ [1, 2, 3, 0, 0, 0, 0]
#F               [0, 0, 1, 2, 3, 0, 0]
#F               [0, 0, 0, 0, 1, 2, 3] ]
#F
#F  Alternatively for downsampled filters valuation can be omitted:
#F   Filt(n, taps, -Length(taps)+1, ds) is equivalent to
#F   DSFilt(n, taps, ds)
#F
Class(Filt, TaggedNonTerminal, DataNonTerminalMixin, rec(
    abbrevs := [
       (n,L)   -> Checked(IsPosIntSym(n), Concat([n], toFiltFunc(L), [1])),
       (n,L,v) -> Checked(IsPosIntSym(n), Concat([n], toFiltFuncV(L,v), [1])),
       (n,L,v,ds) -> Checked(IsPosIntSym(n), IsPosIntSym(ds), Concat([n], toFiltFuncV(L,v), [ds])),
    ],
    isReal := self >> true,

    filtlen := self >> self.params[2].domain(),

    setData := meth(self, newData) self.params[2] := newData; return self; end,

    dims := self >> let(n:=self.params[1], k:=self.params[2].domain(), ds:=self.params[4],
                        [n, ds*(n-1) + k]),

    terminate := self >> let(
	n    := self.params[1],
	coef := self.params[2].lambda().tolist(),
	ds   := self.params[4],
	k    := Length(coef),
	When(n=1, Mat([coef]),
	          RowTensor(n,k-ds,Mat([coef])))),

    HashId := self >> let(h := [ self.params[1], self.params[2].domain(), self.params[4] ],
    When(IsBound(self.tags), Concatenation(h, self.tags), h)),

    print := (self,i,is) >> Print(self.name, "(", PrintCS(self.params), ")",
    When(self.hasTags(), Print(".withTags(", self.getTags(), ")"), "")),
    hashAs := meth(self)
        local hfilt;
        hfilt := Copy(self);
        hfilt.params[2] := FUnk(self.params[2].domain());
        return hfilt;
    end,
    normalizedArithCost := self >> self.params[1] * (2*self.params[2].domain() - 1)
));

#F DSFilt(n, taps, ds) - creates a downsampled filter nonterminal
#F
#F This command is equivalent to Filt(n, taps, -Length(taps)+1, ds)
#F
DSFilt := (n, taps, ds) -> let(
    func := toFiltFunc(taps),
    Filt(n, func[1], func[2], ds)
);

Conv := function(arg)
    local filt, tr, v, deg;
    filt := ApplyFunc(Filt, arg).transpose();
    deg := filt.params[2].domain();
    v   := filt.params[3];
    filt.params[2] := fCompose(filt.params[2], J(deg));
    filt.params[3] := -(v+deg)+1;
    return filt;
end;

_vecList := (lst, t, vlen) -> let(l := Length(lst), Checked(l mod vlen = 0,
    List([0..l/vlen-1], i -> Value(TVect(t, vlen), lst{[1+vlen*i .. vlen*(i+1)]}))
));

_pruneVecZeros := lst -> List(lst, x->When(ForAll(x.v, v->v=0), 0, x));

#  a  b  c  d  0  0  0
#  0  a  b  c  d  0  0
#  0  0  a  b  c  d  0
#  0  0  0  a  b  c  d

#  a  b  c  d  0  0  0  0  0  0
#  0  0  a  b  c  d  0  0  0  0
#  0  0  0  0  a  b  c  d  0  0
#  0  0  0  0  0  0  a  b  c  d

_vecTapsHoriz := function(filt, vlen)
    local zeros, taps, k, vtaps, t, ds;
    t := filt.params[2].range(); # data type
    ds := filt.params[4];
# spiral> _roundup;
# (n, v) -> v * QuoInt((n + v - 1), v)
# _roundup(k + (v-1) ds,  v)
# k + v ds -ds + v - 1
# (k - ds - 1) + v (ds+1)
# (k - ds - 1) / v + v * (ds+1)
    k := filt.params[2].domain();
    zeros := Replicate(vlen*(ds + 1 + QuoInt(k-ds-1, vlen)) - k, 0);
    taps := Concatenation(List(filt.params[2].tolist(), EvalScalar), zeros);
    vtaps := List([0 .. vlen-1], blk ->
    List([0 .. k-1 + Length(zeros)], j -> taps[1 + ((j - blk*ds) mod (k + Length(zeros)))]));
    vtaps := List(vtaps, x->_vecList(x,t,vlen));
    return vtaps;
end;

_vecTaps := function(filt, vlen)
    local zeros, taps, k, vtaps, t, ds;
    t := filt.params[2].range(); # data type
    ds := filt.params[4];
    zeros := Replicate((vlen-1)*ds, 0);
    taps := Concatenation(zeros, List(filt.params[2].tolist(), EvalScalar), zeros);
    k := filt.params[2].domain();
    vtaps := List((vlen-1)*ds + [1..k+(vlen-1)*ds], i -> Value(TVect(t, vlen),
             List([0..vlen-1], j -> taps[i-ds*j])));
    return vtaps;
end;

##
## Rules
##
NewRulesFor(Filt, rec(
    #F Filt_Base: (base case)
    #F
    #F Computes convolution by definition
    #F
    Filt_Base := rec(
    requiredFirstTag := ANoTag,
    info             := "Filt -> Mat",
    forTransposition := true,
    limit            := -1, # means no limit
    applicable     := (self, t) >> When(self.limit < 0, true,
        t.params[1] <= self.limit or t.params[2].domain() <= self.limit),
    apply := (t, C, Nonterms) -> t.terminate()
    ),

    #F Filt_Nest: splits a wider filter into narrower filters, i.e. diagonal stripes
    #F
    Filt_Nest := rec(
	requiredFirstTag := ANoTag,
	forTransposition := false,

        minLength := 2,
        maxLength := 16,
        minTaps := 4,

        libApplicable := (self, t) >> let(n:=t.params[1], ntaps:=t.params[2].domain(),
            leq(self.minLength, n, self.maxLength) * leq(self.minTaps, ntaps)),

	applicable := (self, t) >> let(n:=t.params[1], ntaps:=t.params[2].domain(),
            (IsSymbolic(n) or (n >= self.minLength and n <= self.maxLength)) and
            (IsSymbolic(ntaps) or ntaps >= self.minTaps)),

        freedoms := (self, t) >> let(ntaps := t.params[2].domain(),
	    # we split into filters with 2^k taps, to constrain the search space
	    [ List( Filtered([1..20], i -> 2^(i+1) <= ntaps), i -> [2^i, ntaps mod 2^i] ) ]),

	child := (t, fr) -> let(
	    n  := t.params[1],
	    k  := t.params[2].domain(),
	    ds := t.params[4],
            d := fr[1],
	    Cond(d[2]=0, [ _setN(DSFilt(n, d[1],ds)) ],
                         [ _setN(DSFilt(n, d[1],ds)), _setN(DSFilt(n, d[2], ds)) ])),

	apply := function(t, C, Nonterms)
	    local n,k,kk,nfilts,j,core,tot_wid,bk_wid,rem_wid,taps,
	          bk_filtlen, rem_filtlen, ds;
	    taps := t.params[2];
	    n := t.params[1];
	    k := t.filtlen();
	    ds := t.params[4];
	    tot_wid := ds*(n-1)+k;
	    bk_wid  := Cols(C[1]);

	    bk_filtlen := Nonterms[1].filtlen();
	    nfilts := idiv(k, bk_filtlen);

      	    j := Ind(nfilts);
	    Nonterms[1].setData(fCompose(taps, fAdd(k, bk_filtlen, bk_filtlen * j)));

            core := ISumAcc(j, j.range, Scat(fId(n)) * C[1] *
                                Gath(fAdd(tot_wid, bk_wid, bk_filtlen * j)));

	    if Length(C) = 1 then return core;
	    else
		rem_wid     := Cols(Nonterms[2]);
		rem_filtlen := Nonterms[2].filtlen();
		Nonterms[2].setData(fCompose(taps, fAdd(k, rem_filtlen, bk_filtlen * nfilts)));
		return SUMAcc(core, Scat(fId(n)) * C[2] *
		                    Gath(fAdd(tot_wid, rem_wid, bk_filtlen * nfilts)));
	    fi;
	end
    ),

    #F Filt_OverlapSave: Block convolution using Overlap-Save method.
    #F
    #F Filter_n -> Filter_n/b
    #F
    Filt_OverlapSave := rec(
    info             := "Filt(n,h) -> Filt(b,h)",
    forTransposition := true,
    switch           := true,
    applicable       := t -> t.params[1] > 2 and t.params[2].domain() > 1,
    blockRatio       := 2,
    children      := (self,t) >> let(
        n  := t.params[1],
        k  := t.params[2].domain(),
        v  := t.params[3],
        ds := t.params[4],
        tags := t.getTags(),
        newTags := Cond(tags=[], [],
                    ObjId(tags[1])=AParSMP, Drop(tags, 1),
                tags),
        divs := List( Filtered([1..20], i -> 2^i <= self.blockRatio * k and 2^i < n),
                  i-> [i, n mod 2^i]),
        List(divs, d ->
        When( d[2] = 0,
            [Filt(2^d[1], k, v, ds).withTags(newTags)], # d|n
            [Filt(2^d[1], k, v, ds).withTags(newTags), Filt(d[2], k, v, ds).withTags(newTags)]))), # not(d|n)

    forceBB := 128, # set to 32 to force reasonable unrolling

    apply := meth(self, t, C, Nonterms)
        local n, k, s, ds, bb, tags, smp, dosmp;
        Nonterms[1].setData(t.params[2]);
        n := t.params[1];
        k := t.params[2].domain();
        s := Rows(C[1]);
        ds := t.params[4];
        tags := t.getTags();
        smp := When(tags<>[] and ObjId(tags[1])=AParSMP, s -> SMP(tags[1].p, s), s->s);
        dosmp := tags<>[] and ObjId(tags[1])=AParSMP;

        # NOTE, this allows us not to use .area() for unrolling, since BB
        # explicitly forces unrolling onto small guys
        bb := When(not dosmp and self.forceBB > 0 and k*n <= self.forceBB, BB, x->x);
        if Length(C)=1 then
        return smp(bb(RowTensor(n/s, k-ds, C[1]))); # d|n
        else
        Nonterms[2].setData(t.params[2]);
        return RowDirectSum(k-ds, smp(bb(RowTensor(QuoInt(n,s),k-ds,C[1]))), C[2]); # not(d|n)
        fi;
    end
    ),

  Filt_Vect := rec(
    forTransposition := true,
    switch           := true,
    requiredFirstTag := AVecReg,
    applicable       := t -> (t.params[1] mod t.getTags()[1].v)=0 and t.params[2].domain() >= 1,
    freedoms := t -> [],

    child := (t, freedoms) -> let(
        n  := t.params[1], k := t.params[2].domain(), v := t.params[3], ds := t.params[4],
        vlen := t.getTags()[1].v,
        [ Filt(n/vlen, k+(vlen-1)*ds, v-(vlen-1)*ds, ds*vlen)
         .withTags(Drop(t.getTags(), 1)) ]),

    apply := meth(self, t, C, Nonterms)
        local vlen;
        vlen := t.getTags()[1].v;
        Nonterms[1].setData(FData(_vecTaps(t, vlen)));
        return VTensor(C[1], vlen) * VGath_dup(fId(Cols(t)), vlen);
    end
    ),

  Filt_VectHoriz := rec(
    forTransposition := true,
    switch           := true,
    requiredFirstTag := AVecReg,
    applicable       := t -> let(vlen := t.getTags()[1].v,
        t.params[1] = vlen and
        t.params[2].domain() >= 1
    ),
    freedoms := t -> [],

    child := (t, freedoms) -> let(vlen := t.getTags()[1].v,
        vtaps := _vecTapsHoriz(t, vlen),
        [ TTag(TL(vlen^2, vlen, 1, 1), AVecReg(vlen)),
          Filt(1, vtaps[1]) ]),

    apply := meth(self, t, C, Nonterms)
        local vlen, vtaps, filts, vv;
        vlen := t.getTags()[1].v;
        vtaps := _vecTapsHoriz(t, vlen);
        filts := [];
        for vv in vtaps do
            Nonterms[1].setData(vv);
        Add(filts, Copy(C[1]));
        od;

        return VTensor(RowVec(fConst(vlen, 1)), vlen) *
               C[1] *
           VTensor(VStack(filts), vlen);
    end
    ),

    # Vectorization by converting a filter into a v-way filterbank
    #
    Filt_VectBank := rec(
 	forTransposition := true,
 	switch           := true,
 	requiredFirstTag := AVecReg,
        libApplicable    := t -> eq(imod(t.params[1], t.firstTag().v), 0),
 	applicable       := t -> let(n:=t.params[1], vlen:=t.firstTag().v,
            IsSymbolic(n) or (n mod vlen)=0),

 	freedoms := t -> [],

	child := (self,t,fr) >> let(
 	    n := t.params[1], taps := t.params[2], v := t.params[3], ds := t.params[4], vlen := t.firstTag().v,
 	    [ Filt(n/vlen, taps, v, ds*vlen).withTags(Drop(t.getTags(), 1)) ]),

 	apply := (self, t, C, Nonterms) >> let(
 	    n := t.params[1], taps := t.params[2], ds := t.params[4], vlen := t.firstTag().v,
            #GT(C[1], HH(Cols(t), Cols(C[1]), 0, [1, ds]),
            #         HH(Rows(t), Rows(C[1]), 0, [vlen, 1]), [vlen]).withTags(new_tags)
            VTensor(C[1], vlen) *
            VGath_u(H(Cols(t), Cols(C[1]), 0, 1), vlen)
                     ),
    ),

));

RulesFor(Filt, rec(
    #F Filter_OverlapSaveFreq: Block convolution using Overlap-Save method suited
    #F for Frequency-domain rules. It blocks so that we get powers of two for
    #F the columns -> leading to Circulant rule.
    #F
    #F Filter_n -> Filter_n/b
    #F
    Filt_OverlapSaveFreq := rec(
    info             := "Filt(n,h) -> Filt(b,h)",
    forTransposition := true,
    switch           := true,
    isApplicable     := P -> let(n:=P[1], k:= P[2].domain(), x := Log2Int(n+k-2), ds := P[4],
                                     (ds=1) and (n > 2) and (k > 1) and (k < 2^x)),
    allChildren      := P -> let(
        n  := P[1],
        k  := P[2].domain(),
         m := n + k - 1,
        divs := List( Filtered([1..20], i -> 2^(i-3) <= k and 2^i >= k and 2^i < n+k-1),
                  i-> [i, n mod (2^i-k+1)]),
        List(divs, d->
        When( d[2] = 0,
            [Filt(2^d[1]-k+1, P[2], P[3])], # d|n
            [Filt(2^d[1]-k+1, P[2], P[3]), Filt(d[2], P[2], P[3])]))), # not(d|n)

    rule := (P, C) -> let(
        n := P[1],
        k := Length(P[2]),
        s := C[1].dimensions[1],
        When(Length(C)=1,
        RowTensor(n/s, k-1, C[1]), # d|n
        RowDirectSum(k-1, RowTensor(QuoInt(n,s),k-1,C[1]), C[2]))# not(d|n)
    )
    ),

    #F Filt_OverlapAdd: block convolution using Overlap-Add method.
    #F
    #F Filt_n -> Conv_n/b, Toeplitz, Toeplitz
    #F
    Filt_OverlapAdd := rec(
    info              := "Filt(n,h) -> Conv(b,h)",
    forTransposition  := true,
    isApplicable      := P -> P[1] > 2 and P[2].domain() <= P[1]+1 and P[4]=1,
    switch := false,
    allChildren := function(P)
        local n,l,divs,d, coef,list1,list2;
        coef := P[2].tolist();
        l := Length(coef);
        n := P[1]-l+1;
        divs := List(Filtered([1..20], i-> 2^(i-1) < l and 2^i < P[1]),
                 i-> [i, n mod 2^i]);
        list1:= Reversed(Concat(Replicate(l-2,0),coef{[1..l-1]}));
        list2:= Reversed(Concat(coef{[2..l]},Replicate(l-2,0)));
        return List(divs, d->
        When(d[2] = 0,
             [ Conv(2^d[1], coef, P[3]), Toeplitz(list1), Toeplitz(list2) ], # d|n
             [ Conv(2^d[1], coef, P[3]), Conv(d[2],coef, P[3]),            # not(d|n)
                                     Toeplitz(list1), Toeplitz(list2) ]));
    end,

    rule := (P, C) -> let(
        k := P[2].domain(),
        n := P[1] - k + 1,
        s := Cols(C[1]),
        When(Length(C)=3,
        ColDirectSum(k-1, C[2], ColTensor(n/s, k-1, C[1]), C[3]), # d|n
        ColDirectSum(k-1, C[3], ColTensor(QuoInt(n,s), k-1, C[1]), C[2], C[4])) # not(d|n)
    )
    ),

    #F Filt_OverlapAdd: block convolution using Overlap-Add method.
    #F
    #F Filt_n -> Ext * Conv_n/b
    #F
    Filt_OverlapAdd2 := rec(
    info              := "Filt(n,h) -> Conv(b,h)",
    switch            := false,
    forTransposition  := true,
    isApplicable      := P -> P[1] > 2 and P[2].domain() > 1 and P[4]=1,
    allChildren := function(P)
        local n,k,divs,d, coef;
        k := P[2].domain();
        n := P[1]+k-1;
        coef := P[2];
        divs := List( Filtered([1..20], i -> 2^(i-1) < k and 2^i < P[1]),
                  i -> [i, n mod 2^i]);
        return List(divs, d->
        When( d[2] = 0,
            [Conv(2^d[1], coef, P[3])], # d|n
            [Conv(2^d[1], coef, P[3]), Conv(d[2], coef, P[3])])); # not(d|n)
    end,
    rule := (P, C) -> let(
        k := P[2].domain(),
        n := P[1] + k - 1,
        s := Cols(C[1]),
        ext := Ext(P[1],k-1,k-1).transpose(),
        When(Length(C)=1,
        ext * ColTensor(n/s, k-1, C[1]), # d|n
        ext * ColDirectSum(k-1, ColTensor(QuoInt(n,s), k-1, C[1]), C[2]))# not(d|n)
    )
    ),

    #F Filt_KaratsubaSimple:
    #F
    #F Filt_n(h) -> Filt_n(h0) dirsum Filt_n(h1) dirsum Filt_n(h0+h1)
    #F
    #F h0,h1 - even and odd downsampled h
    #F This is radix-2 Karatsuba - "transposed graph" version.
    #F Both the size and the length of the filter have to be divisible by 2
    #F
    Filt_KaratsubaSimple:= rec(
    info              := "Filt(n,h) -> Filt(b,h_e) dirsum Filt(b,h_o) dirsum Filt(b,h_e+h_o)",
    forTransposition  := false,
    switch := false,

    lowerLimit := 8,
    upperLimit := 128,

    isApplicable := (self, P) >> let(n := P[1], k := P[2].domain(), ds := P[4],
        (ds=1) and
        (n >= self.lowerLimit and n <= self.upperLimit) and (n mod 2 = 0) and
        (k >= self.lowerLimit and k <= self.upperLimit) and (k mod 2 = 0) and
        not IsBound(P[2].nested)),

    allChildren := P -> let(n := P[1], k := P[2].domain(), v := Int(P[3]/2),
        [[ Filt(n/2, k/2, v) ]]),

    rule := function( P, C, Nonterms )
        local n, m, l, m0, S, R, v0, v1, h0, h1, hsum, deg, coeffs, hds, f1,f2,f3;
        [n, l, deg] := [P[1], P[2].domain(), P[3]];
        coeffs := List(P[2].tolist(), x->x.ev());
        hds := DownsampleTwo([coeffs,deg]);
        [h0,v0] := hds[1];
        [h1,v1] := hds[2];

            if v0 <> v1 then Error("valuations of downsampled filters aren't the same"); fi;

        if (deg mod 2 = 1) then [h0,h1] := [h1,h0]; fi;
        hsum := ListSum(h0, h1);

        Nonterms[1].setData(FData(h0)); f1 := Copy(C[1]);
        Nonterms[1].setData(FData(h1)); f2 := Copy(C[1]);
        Nonterms[1].setData(FData(hsum)); f3 := Copy(C[1]);

        m0 := n + l - 1;
        m := (m0+1)/2 - 1;
        S := Tensor(Mat([[1,0,1],[0,1,1]]), I(n/2));
        # NOTE: use 3-band filter bank (instead of 3 separate filters) here to improve locality
        return L(n,n/2) * S *
               VStack(f1 * HStack(I(m), O(m,1), -I(m)),
                  f2 * HStack(O(m,1), I(m), -I(m)),
              f3 * Gath(fAdd(m0,m,m+1))) *
           BB(OS(m0,2));
        end
    ),

    #F Filt_Circulant
    #F
    #F Computes convolution by definition
    #F
    Filt_Circulant := rec(
    info             := "Filt -> Circulant",
    forTransposition := true,
    isApplicable     := P -> let(n:= P[1], k:= P[2].domain(), m:= n+k-1, ds:=P[4],
                               ds = 1 and m = 2^(Log2Int(m)) and n/k <= 8 and k/n <= 2),
    allChildren      := P -> let(
             n := P[1], l := P[3], coef := P[2], r := l+coef.domain()-1,
         [[ Circulant(n+PosInt(-l)+ PosInt(r), FUnk(coef.domain()), l) ]]),

    rule := (P,C,Nonterms) -> let(
        n := P[1], l := P[3], r := l + P[2].domain() - 1,
        nt := Nonterms[1].setData(P[2]),
        Ext(n, PosInt(-l), PosInt(r)).transpose() * C[1] )
    ),

    #F Filt_Blocking
    #F
    #F Linear Convolution (Filter) -> Blocks (Toeplitzes)
    #F
    #F Filter is partitioned into square blocks of sizes given by all
    #F proper divisors of the input size.
    #F
    Filt_Blocking := rec(
    info             := "Filt(n,h) -> Toeplitz(b,h')",
    forTransposition := false,
    switch           := true,
    isApplicable     := P -> let(n := P[1], k := P[2].domain(), ds := P[4],
        ds = 1 and n > 2 and n < k and ObjId(P[2]) in [FData, FUnk, FList]),

        # split into dense and half-sparse toeplitz, rest will be unrolled
    allChildren := P -> let(n := P[1], k := P[2].domain(), b := n,
        [[ Toeplitz(2*b-1),
           Toeplitz(2*b-1, 0, b-1),
           Toeplitz(2*b-1, b, 2*b-1) ]]),

    #When((k-1) mod b <>0, [ Toeplitz(2*((k-1) mod b) - 1) ], [ ]),

    rule := function(P, C, Nonterms)
        local n, k, b, Ls, Lsdata, S, j, rem, ofs, jj;
        n  := P[1];
        k  := P[2].domain();
        b := n;

            # first block major part
        Ls := List(P[2].tolist(), x->x.eval());
        Ls := Concat(Replicate(b-1,0), Ls, Replicate(2*b - k mod b,0));
        Lsdata := FData(Ls);

        j := Ind(Int((k-1)/b)-1); # we unroll first and last iteration
        jj := Ind(Int((k-1)/b)+1); # we unroll first and last iteration
        Nonterms[1].setData(fCompose(Lsdata, fAdd(Lsdata.domain(), 2*b-1, (j+1)*b), J(2*b-1)));
        Nonterms[2].setData(fCompose(Lsdata, fAdd(Lsdata.domain(), 2*b-1, jj*b), J(2*b-1)));
        Nonterms[3].setData(fCompose(Lsdata, fAdd(Lsdata.domain(), 2*b-1, jj*b), J(2*b-1)));
        Nonterms[3].setNonZero(b - (k mod b), 2*b - 1);

        S := SUMAcc(
        Data(jj, V(0),
            Scat(fId(b)) * C[2] * Gath(fTensor(fBase(j.range+2,jj), fId(b)))),
        ISumAcc(j, j.range,
            Scat(fId(b)) * C[1] * Gath(fTensor(fBase(j.range+2,j+1), fId(b)))),
        Data(jj, V(j.range+1),
            Scat(fId(b)) * C[3] * Gath(fTensor(fBase(j.range+2,jj), fId(b))))
        );

            # we have to add a small remaining triangle to the right
        if (k-1)mod b <> 0 then
        ofs := b * Int(((k-1)/b)+1) + 1;
        rem := Reversed( Concat( Ls{[ofs .. ofs + (k-1)mod b - 1]}, Replicate((k-1)mod b -1, 0)));
        S := ColDirectSum((b+k-1) mod b, S, toeplitz(rem));
        fi;

        return S;
    end
    ),

    Filt_KaratsubaFast := rec(
    info              := "Filt(n,h) -> Filt(b,h_e) dirsum Filt(b,h_o) dirsum Filt(b,h_e+h_o)",
    forTransposition  := false,
    switch := false,

    lowerLimit := 8,
    upperLimit := 128,

    isApplicable := (self, P) >> let(n := P[1], k := P[2].domain(), ds := P[4],
        (ds=1) and
        (n >= self.lowerLimit and n <= self.upperLimit) and (n mod 2 = 0) and
        (k >= self.lowerLimit and k <= self.upperLimit) and (k mod 2 = 0)),

    allChildren := P -> let(n := P[1], k := P[2].domain(), v := Int(P[3]/2),
        [[ Filt(n/2, k/2, v) ]]),

    rule := function( P, C, Nonterms )
        local m, m0, S, R, l, h0, h1, hsum, deg, n, t, v, coeffs, newcoeffs, hds, j, dirsum, gfunc, one, minone;
        n := P[1]; l := P[2].domain(); deg := P[3];
        coeffs := List(P[2].tolist(), x->x.ev());
        hds := DownsampleTwo([coeffs,deg]);
            if hds[1][2] <> hds[2][2] then Error("valuations of downsampled filters aren't the same"); fi;

        v := hds[1][2]; h0 := hds[1][1]; h1 := hds[2][1];
        if (deg mod 2 = 1) then t := h1; h1 := h0; h0 := t; fi;
        hsum := ListSum(h0, h1);

        newcoeffs := FData(Concat(h0, h1, hsum));
        j := Ind(3);
        gfunc := fCompose(newcoeffs, fTensor(fBase(j), fId(l/2)));
        Nonterms[1].setData(gfunc);
        dirsum := ISum(j, j.range,
            Scat(fTensor(fBase(j), fId(n/2))) * C[1] * Gath(fTensor(fBase(j), fId(Cols(C[1])))));
        m0 := n + l - 1;
        m := (m0+1)/2 - 1;
        S := Tensor(Mat([[1,0,1],[0,1,1]]), I(n/2));

        one := Diag(fConst(m,1));
        minone := Diag(fConst(m,-1));

        R := VStack(
            SUMAcc(one*Gath(fAdd(m0,m,0)), minone*Gath(fAdd(m0,m,m+1))),
            SUMAcc(one*Gath(fAdd(m0,m,1)), minone*Gath(fAdd(m0,m,m+1))),
        one*Gath(fAdd(m0,m,m+1))
        );

        return L(n,n/2) * S * dirsum * R * OS(m0,2);
        end
    )

));


text:
Mean := function(v)
  local n;
  n := Length(v);
  if n = 0 then
    return 0;
  else
    return Sum(v)/n;
  fi;
end;

Mean([3, 1, 4, 1, 5, 9]);
# 23/6


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Import(rewrite);

VectorStrategySum := [
    MergedRuleSet(JoinDirectSums,StretchRaderMid),
    MergedRuleSet(StandardSumsRules, RulesPropagate),
    # when to terminate Grp() has become an issue: too early -> diags go the wrong way; too late->Rader stretch breaks
	RulesSMP,
    MergedRuleSet(StandardSumsRules, RulesPropagate, StretchRaderMid, JoinDirectSums, RulesTermGrp),
    MergedRuleSet(StandardSumsRules, RulesPropagate, StretchRaderMid, JoinDirectSums, RulesTermGrp, RulesVDiag),
    MergedRuleSet(StandardSumsRules, RulesPropagate, StretchRaderMid, JoinDirectSums, RulesTermGrp, RulesSMP, RulesVec),
    RulesKickout
];

VectorStrategyTerm := [
#    MergedRuleSet(RulesConjTerm, RulesFuncSimp),
    MergedRuleSet(StretchRaderMid, JoinDirectSums),
    RulesTerm, RulesTermGrp, RulesKickout,
];

VectorStrategyTerm2 := [
#    RulesConjTerm,
    TerminateDirectSums, TerminateRaderMid, TerminateDirectSums, TerminateSymSPL, TerminateDPrm,
];

VectorStrategyVDiag := [
    RulesVDiag, RulesVec, RulesSMP, RulesPropagate,
    StandardSumsRules
#   MergedRuleSet(RulesVDiag, StandardSumsRules)
];

VectorStrategyRC := [
    MergedRuleSet(RulesSMP, RulesSplitComplex),
    RulesKickout,
    RulesVBlkInt,
    RulesVDiag,
    MergedRuleSet(RulesPropagate, RulesSMP, RulesSplitComplex, RulesVRC, RulesSums),
    MergedRuleSet(RulesPropagate, RulesVRC, RulesCxRC, RulesVRCTerm, RulesSMP, RulesVec, TerminateSymSPL, RulesRC, RulesSplitComplex),
    RulesTerm,
    RulesTermGrp,
    RulesSMP,
    RulesVec,
    RulesPropagate,
    StandardSumsRules,
    RulesKickout
];

VectorStrategyRCUnterm := [
    MergedRuleSet(RulesVRC, RulesCxRC, RulesRC),
    MergedRuleSet(RulesSMP, RulesVec, RulesCxRC, RulesPropagate)
];

# VectorStrategy := [
#     MergedRuleSet(StandardSumsRules)
#     MergedRuleSet(RulesVec, RulesPropagate, RulesDiag),
#     RulesVDiag,
#     MergedRuleSet(RulesVec, RulesPropagate),
#     MergedRuleSet(StandardSumsRules)
#     MergedRuleSet(RulesVRC, RulesCxRC, RulesRC), RulesVec,
#     RulesKickout,
#     MergedRuleSet(StandardSumsRules)
# ];

VectorStrategyRCRec := [
    MergedRuleSet(RulesVRC, RulesCxRC, RulesRC),     RulesKickout,
    MergedRuleSet(RulesSMP, RulesVec, RulesPropagate),
    StandardSumsRules
];

fix_fAdd := [
    MergedRuleSet(StandardSumsRules, Rules_fAdd),
];


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details
#
# n-qubit X transforms
#

#F qXT( n ) - X Gate non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies an n-point X Transform to the target qubits
#F Note:       qXT(1) denotes the matrix[[0, 1], [1, 0]],
#F             qXT(_) is symmetric.
#F
#F qXT( n ) -> an n-qubit X transform
Class(qXT, TaggedNonTerminal, rec(
  abbrevs   := [ n -> Checked(IsPosInt(n), [n]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> Tensor(Replicate(self.params[1], qX())),
  isReal    := self >> true,
  groups    := self >> [Replicate(self.params[1], 1)],
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1]),
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1]),
  connected := self >> true,
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Complex(64),
));

NewRulesFor(qXT, rec(

    # qXT_BinSplit rule
    # qXT_BinSplit qXT_(k) -> (qXT_(k1) tensor I) (I tensor qXT_(k2))
    # k1 + k2 = k
    qXT_BinSplit := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] > 1,
        children         := nt -> List( [1..nt.params[1] - 1],  i -> [ Tensor(qXT(i), qXT(nt.params[1]-i)).withTags(nt.getTags()) ]  ),
        apply            := (nt, c, cnt) -> c[1],
        switch := true,
    ),

    #F qXT_Base: qXT(1) = qX() SPL object
    #F Directly represent as an implementable gate
    qXT_Base := rec(
        info             := "qXT_(1) -> qX()",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1,
        apply            := (nt, c, cnt) -> qX(),
    )

));

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#the multiple number of buffers is equal as size
multiple_buffers := function(elements)
                local number, s;

                number := Length(elements);

                s := elements[1];

                return [ var("S", TArray(TArray(TReal, s.t.size), 2)) ];
end;

construct_loop := function(i, low, high, ld, st, val, br, s1, s2, s)
               local new_chain, new_loop, new_high;

               new_chain := chain(
                         SubstVars(Copy(st), rec((s2.id) := nth(s, imod(i,2)))),
                         SubstVars(Copy(ld), rec((i.id) := add(i,val), (s1.id) := nth(s, imod(i,2)))),
                         br);
               new_high := high - val;
               new_loop := swp_loop(i, [low..new_high.v], new_chain);

               return new_loop;
end;

swp_peel := function(fun, opts)
         #peeling instruction from loops
         local l, loops, ld, st, br, i, low, high, prologue, epilogue, new_chain, new_fun, s1, s2, svars, s;

         svars := [];

         loops := Collect(fun, swp_loop);

         new_fun := Copy(fun);

         for l in loops do
             i := l.var;
             low := Minimum(l.range);
             high := Maximum(l.range);

             ld := Collect(l, [loop, @1, @2, dma_load]);
             if Length(ld) = 0 then ld := Collect(l, dma_load)[1]; else ld := ld[1]; fi;
             st := Collect(l, [loop, @1, @2, dma_store]);
             if Length(st) = 0 then st := Collect(l, dma_store)[1]; else st := st[1]; fi;
             br := Collect(l, barrier_cmd)[1];

             s1 := Filtered(ld.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[1];
             s2 := Filtered(st.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[1];
             Add(svars, s1);
             Add(svars, s2);

             s := multiple_buffers(svars)[1];

             prologue := chain(
                      SubstVars(Copy(ld), rec((i.id) := low, (s1.id) := nth(s, low))),
                      br,
                      SubstVars(Copy(ld), rec((i.id) := low + 1, (s1.id) := nth(s, low + 1))),
                      br);

             epilogue := chain(
                      SubstVars(Copy(st), rec((i.id) := high - 1, (s2.id) := nth(s, bin_and(high - 1, 1)))),
                      br,
                      SubstVars(Copy(st), rec((i.id) := high, (s2.id) := nth(s, bin_and(high, 1)))));

             new_chain := When((low = 0 and high = 1), chain(prologue, epilogue), chain(prologue, construct_loop(i, low, high, ld, st, V(2), br, s1, s2, s), epilogue));
             new_chain := RulesStrengthReduce(new_chain);
             new_fun := SubstBottomUp(new_fun, @(1, swp_loop, e->e.var = i), e->new_chain);
         od;

         return new_fun;
end;

swp_merge_variables := function(fun, opts)
                   local l, loops, i, new_loop, new_fun, s1, s2, svars, s;

                   svars := [];

                   loops := Collect(fun, swp_loop);

                   new_fun := Copy(fun);

                   for l in loops do
                       i := l.var;

                       s1 := Filtered(l.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[1];
                       s2 := Filtered(l.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[2];
                       Add(svars, s1);
                       Add(svars, s2);

                       s := multiple_buffers(svars)[1];

                       new_loop := SubstVars(Copy(l), rec((s1.id) := nth(s, imod(i, 2)), (s2.id) := nth(s, imod(i, 2))));

                       new_loop := RulesStrengthReduce(new_loop);
                       new_fun := SubstBottomUp(new_fun, @(1, swp_loop, e->e.var = i), e->new_loop);
                   od;

                   return new_fun;
end;

Class(SWPDMACodegen, DefaultCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, io, t, initcode;
        [x, y] := self.initXY(x, y, opts);

        Add(x.t.qualifiers, opts.memModifier);
        Add(y.t.qualifiers, opts.memModifier);

        o := o.child(1);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);
        [icode,t] := UTimedAction(ESReduce(icode, opts)); #PrintLine("ESReduce ", t);
        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        #PrintLine("---compile--");

        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        io := When(x=y, [x], [y, x]);
        sub := Cond(IsBound(opts.subNameDMA), opts.subNameDMA, "sub_dma");
        icode := func(TVoid, sub, Concatenation(io, params), icode);

        return icode;
    end,

    LSKernel := meth(self, o, y, x, opts)
        return barrier_cmd(self.swp_var);
    end,

    DMAGath := meth(self, o, y, x, opts)
        local i, func, rfunc, ix, size;
        Add(self.loadbuffers, y);
        Add(self.membuffers, x);
        i := Ind();

        func := o.func;
        size := 1;
        if ObjId(func) = fTensor and ObjId(Last(func.children())) = fId then
            size := Last(func.children()).domain();
            func := fTensor(Concat(DropLast(func.children(), 1), [fBase(size, 0)]));
        fi;

        func := func.lambda();
        return loop(i, func.domain(), dma_load(y+(i*size), x+func.at(i), size));
    end,

    DMAScat := meth(self, o, y, x, opts)
        local i, func, rfunc, ix, size;

        Add(self.storebuffers, x);
        Add(self.membuffers, y);
        i := Ind();

        func := o.func;
        size := 1;
        if ObjId(func) = fTensor and ObjId(Last(func.children())) = fId then
            size := Last(func.children()).domain();
            func := fTensor(Concat(DropLast(func.children(), 1), [fBase(size, 0)]));
        fi;

        func := func.lambda();
        return loop(i, func.domain(), dma_store(y+func.at(i), x+(i*size), size));
    end,

    DMAFence := (self,o,y,x,opts) >> chain(self(o.child(1), y, x, opts), dma_fence()),

    swp_var := false,
    SWPSum := meth(self, o, y, x, opts)
        return swp_loop(o.var, o.domain, self(o.child(1), y, x, opts));
    end

));

Class(SWPCPUCodegen, DefaultCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, io, t, initcode;
        [x, y] := self.initXY(x, y, opts);

        o := o.child(1);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);
        [icode,t] := UTimedAction(ESReduce(icode, opts)); #PrintLine("ESReduce ", t);
        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        #PrintLine("---compile--");

        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        io := When(x=y, [x], [y, x]);
        Add(params, Filtered(icode.free(), i->not IsBound(i.init)));
        sub := Cond(IsBound(opts.subNameCompute), opts.subNameCompute, "sub_cpu");

        icode := func(TVoid, sub, params, icode);

        return icode;

    end,

    DMAGath := meth(self, o, y, x, opts)
        local dma;
        if not IsBound(x.t.qualifiers) then x.t.qualifiers := []; fi;
        if not IsBound(y.t.qualifiers) then y.t.qualifiers := []; fi;
        Add(x.t.qualifiers, opts.memModifier);
        y.t.qualifiers := [ opts.scratchModifier ];
        Add(self.loadbuffers, y);
        return nop_cmd(self.swp_var);
    end,

    DMAScat := meth(self, o, y, x, opts)
        local dma;
        if not IsBound(x.t.qualifiers) then x.t.qualifiers := []; fi;
        if not IsBound(y.t.qualifiers) then y.t.qualifiers := []; fi;
        x.t.qualifiers := [ opts.scratchModifier ];
        Add(y.t.qualifiers, opts.memModifier);
        Add(self.storebuffers, x);
        return barrier_cmd(self.swp_var);
    end,

    LSKernel := (self, o, y, x, opts) >> self(o.child(1), y, x, opts),

    DMAFence := (self,o,y,x,opts) >> chain(barrier_cmd(self.swp_var), self(o.child(1), y, x, opts)),

    swp_var := false,
    SWPSum := meth(self, o, y, x, opts)
        return swp_loop(o.var, o.domain, self(o.child(1), y, x, opts));
    end
));

Class(SWPBarrierScratchCodegen, DefaultCodegen, rec(
    CPUCodegen := SWPCPUCodegen,
    DMACodegen := SWPDMACodegen,
    MainCodegen := ScratchMainCodegen,

    Formula := meth(self, o, y, x, opts)
        local compute_time, tag, main, initfunc, cpufunc, dmafunc, prog, params, sub, initsub, memvars, scratchvars, v, io, loadvar, storevar, v, substrec, svars, initcode, datas, dvars, dv, counter;

        datas := Collect(o, FDataOfs);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");

        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            initcode := chain(List(datas, x -> SReduce(x.var.init, opts)));
            initfunc := func(TVoid, initsub, params :: Set(Collect(initcode, param)), initcode);
        else
            initfunc := func(TVoid, initsub, params, chain());
        fi;

        dvars := List(datas, x->x.var);
        for dv in dvars do
            dv.t.qualifiers := [opts.romModifier];
        od;

        self.DMACodegen.loadbuffers := Set([]);
        self.DMACodegen.storebuffers := Set([]);
        self.DMACodegen.membuffers := Set([]);
        dmafunc := self.DMACodegen.Formula(o, y, x, opts);
        self.DMACodegen.membuffers := Set(self.DMACodegen.membuffers);
        SubtractSet(self.DMACodegen.membuffers, Set([x, y]));
        for v in self.DMACodegen.membuffers do v.t.qualifiers := [opts.memModifier]; od;

        self.CPUCodegen.loadbuffers := Set([]);
        self.CPUCodegen.storebuffers := Set([]);
        cpufunc := self.CPUCodegen.Formula(o, y, x, opts);

        memvars := Set(Concat(
            Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.memModifier in j.t.qualifiers),
            self.DMACodegen.membuffers));
        scratchvars := Set(Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.scratchModifier in j.t.qualifiers));

        tag := opts.tags[1];

        loadvar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        loadvar.t.qualifiers := [opts.scratchModifier];
        storevar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        storevar.t.qualifiers := [opts.scratchModifier];

        substrec := rec();
        for v in self.CPUCodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.CPUCodegen.storebuffers do substrec.(v.id) := storevar; od;
        for v in self.DMACodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.DMACodegen.storebuffers do substrec.(v.id) := storevar; od;
        cpufunc :=  SubstVars(cpufunc, substrec);
        dmafunc :=  SubstVars(dmafunc, substrec);

        cpufunc := SubstTopDown(cpufunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));
        dmafunc := SubstTopDown(dmafunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));

        [cpufunc, dmafunc, svars] := When(opts.swp,[swp_merge_variables(cpufunc, opts), swp_peel(dmafunc, opts), multiple_buffers([loadvar, storevar])],[cpufunc,dmafunc,[loadvar,storevar]]);

        x.t.qualifiers := Set(x.t.qualifiers);
        y.t.qualifiers := Set(y.t.qualifiers);

        compute_time := true;

        counter := var("count",TInt);

        cpufunc := func(cpufunc.ret, cpufunc.id, When(compute_time, Concat(cpufunc.params, [counter]), cpufunc.params), chain(register(), cpufunc.cmd));
	    dmafunc := func(dmafunc.ret, dmafunc.id, When(compute_time, Concat(dmafunc.params, [counter]), dmafunc.params), chain(register(), dmafunc.cmd));
        initfunc := func(initfunc.ret, initfunc.id, initfunc.params, chain(initialization(), initfunc.cmd));
        #PrintLine(svars);
        #main := self.MainCodegen.genMain(opts, sub, dmafunc, cpufunc, [y,x], params);
        prog := program(
            decl(Concat(Set(memvars), svars, dvars), chain(
                initfunc,
                cpufunc,
		        dmafunc
        )));

        prog.dimensions := o.dims();
        return prog;
    end
));