I have a system of equation like bellow:

"L3 + L4 + S5 + S12 + L1 + D4 + L8 + S3 + L7 + D8 + D5 + L5 == 1",
"L4 + D9 + S5 + L1 + D16 + L8 + L6 + S8 + L7 + D8 == 1",
"L4 + L1 + D16 + S60 + L2 == 1",
"L3 + D12 + L1 + S9 + S3 + D5 + S105 + L2 + L7 + D28 + L5 == 1",
"S11 + L3 + S72 + D10 + D72 + D9 + S5 + D16 + S9 + S60 + L6 + S105 + L2 + L8 + D5 == 1",
"L3 + S60 + L2 + L4 == 1",
"D72 + L6 + S105 + L7 + D28 + L5 == 1",
"S72 + D72 + L8 + L6 + L5 == 1",
"D4 + S12 + S11 + D10 == 1",
"D12 + D10 + S9 + S8 + D8 + S12 == 1",
"S11 + D12 + S72 + D9 + D4 + S3 + S8 + D28 == 1",
"S11 + D12 + D10 + D72 + D9 + D28 + S72 + S5 + S12 + D4 + D16 + S9 + S3 + S60 + S8 + S105 + D8 + D5 == 1"

and i need all possible solutions where solution are in [0 or 1] And solution would be like : S1 = [L3,D12, S9] S2 = [L4, D72, S8] solution should take on element from L, D and S. I need to make this code on python, and I was trying to model the equation system as matrix : M_L*L + M_D*D + M_S*S = Ones but dont have a hint for coding that. Can somebody help please ?

I tried loops and it was not easy to implement. I tried as well the library sympy: here is my code the issue is that i dont have combination of solution and should add constraint

from sympy import symbols, Eq, solve

L1, L2, L3, L4, L5, L6, L7, L8, S3, S5, S8, S9, S11, S12, S60, S72, S105, D4, D5, D8, D9, D10, D12, D16, D28, D72 = symbols('L1 L2 L3 L4 L5 L6 L7 L8 S3 S5 S8 S9 S11 S12 S60 S72 S105 D4 D5 D8 D9 D10 D12 D16 D28 D72', boolean=True)


eq1 = Eq(L3 + L4 + S5 + S12 + L1 + D4 + L8 + S3 + L7 + D8 + D5 + L5, 1)
eq2 = Eq(L4 + D9 + S5 + L1 + D16 + L8 + L6 + S8 + L7 + D8, 1)
eq3 = Eq(L4 + L1 + D16 + S60 + L2, 1)
eq4 = Eq(L3 + D12 + L1 + S9 + S3 + D5 + S105 + L2 + L7 + D28 + L5, 1)
eq5 = Eq(S11 + L3 + S72 + D10 + D72 + D9 + S5 + D16 + S9 + S60 + L6 + S105 + L2 + L8 + D5, 1)
eq6 = Eq(L3 + S60 + L2 + L4, 1)
eq7 = Eq(D72 + L6 + S105 + L7 + D28 + L5, 1)
eq8 = Eq(S72 + D72 + L8 + L6 + L5, 1)
eq9 = Eq(D4 + S12 + S11 + D10, 1)
eq10 = Eq(D12 + D10 + S9 + S8 + D8 + S12, 1)
eq11 = Eq(S11 + D12 + S72 + D9 + D4 + S3 + S8 + D28, 1)
eq12 = Eq(S11 + D12 + D10 + D72 + D9 + D28 + S72 + S5 + S12 + D4 + D16 + S9 + S3 + S60 + S8 + S105 + D8 + D5, 1)

solution = solve((eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12))

print(solution)

You could try the z3 smt solver. A BitVec of length 1 can hold the values 0 and 1.

Bit vectors of length 1

The following example code uses bitvectors of length 1. If I understand correctly, the sums will be done modulo 2.

from z3 import Solver, BitVec, sat, Or

s = Solver()

param = [BitVec(name, 1) for name in
         'L1 L2 L3 L4 L5 L6 L7 L8 S3 S5 S8 S9 S11 S12 S60 S72 S105 D4 D5 D8 D9 D10 D12 D16 D28 D72'.split()]
L1, L2, L3, L4, L5, L6, L7, L8, S3, S5, S8, S9, S11, S12, S60, S72, S105, D4, D5, D8, D9, D10, D12, D16, D28, D72 = param

s.add(L3 + L4 + S5 + S12 + L1 + D4 + L8 + S3 + L7 + D8 + D5 + L5 == 1)
s.add(L4 + D9 + S5 + L1 + D16 + L8 + L6 + S8 + L7 + D8 == 1)
s.add(L4 + L1 + D16 + S60 + L2 == 1)
s.add(L3 + D12 + L1 + S9 + S3 + D5 + S105 + L2 + L7 + D28 + L5 == 1)
s.add(S11 + L3 + S72 + D10 + D72 + D9 + S5 + D16 + S9 + S60 + L6 + S105 + L2 + L8 + D5 == 1)
s.add(L3 + S60 + L2 + L4 == 1)
s.add(D72 + L6 + S105 + L7 + D28 + L5 == 1)
s.add(S72 + D72 + L8 + L6 + L5 == 1)
s.add(D4 + S12 + S11 + D10 == 1)
s.add(D12 + D10 + S9 + S8 + D8 + S12 == 1)
s.add(S11 + D12 + S72 + D9 + D4 + S3 + S8 + D28 == 1)
s.add(S11 + D12 + D10 + D72 + D9 + D28 + S72 + S5 + S12 + D4 + D16 + S9 + S3 + S60 + S8 + S105 + D8 + D5 == 1)
res = s.check()
count = 0
while res == sat:
    count += 1
    M = s.model()
    sol = [M[p] for p in param]
    # print(count, sol)
    print(count, ", ".join([f"{p}:{M[p]}" for p in param]))
    s.add(Or([p != M[p] for p in param]))
    res = s.check()

The output starts with

1 L1:0, L2:0, L3:0, L4:0, L5:0, L6:0, L7:0, L8:1, S3:0, S5:1, S8:1, S9:1, S11:0, S12:1, S60:1, S72:1, S105:1, D4:0, D5:0, D8:0, D9:0, D10:0, D12:0, D16:0, D28:1, D72:1
2 L1:1, L2:1, L3:1, L4:1, L5:1, L6:1, L7:1, L8:1, S3:1, S5:0, S8:0, S9:0, S11:1, S12:0, S60:0, S72:0, S105:0, D4:1, D5:1, D8:0, D9:0, D10:1, D12:0, D16:0, D28:0, D72:0
3 L1:1, L2:1, L3:1, L4:1, L5:1, L6:1, L7:0, L8:0, S3:1, S5:0, S8:0, S9:0, S11:1, S12:0, S60:0, S72:1, S105:0, D4:1, D5:1, D8:0, D9:0, D10:1, D12:0, D16:0, D28:1, D72:0
4 L1:1, L2:0, L3:1, L4:0, L5:0, L6:0, L7:0, L8:0, S3:1, S5:0, S8:0, S9:0, S11:1, S12:0, S60:0, S72:1, S105:0, D4:1, D5:1, D8:0, D9:0, D10:1, D12:0, D16:0, D28:1, D72:0
5 L1:1, L2:0, L3:1, L4:0, L5:0, L6:0, L7:0, L8:0, S3:0, S5:1, S8:0, S9:0, S11:1, S12:0, S60:0, S72:0, S105:0, D4:1, D5:1, D8:0, D9:1, D10:1, D12:0, D16:0, D28:0, D72:1
...

Bit vectors of length 8 with 0 or 1 constraint

If you change it to bit vectors of length 8, and add a constraint that each value should be either 0 or 1, no solutions are found. If you'd limit the equations to e.g. the first 8, 3 solutions are found.

from z3 import Solver, BitVec, sat, Or, And

s = Solver()

param = [BitVec(name, 8) for name in
         'L1 L2 L3 L4 L5 L6 L7 L8 S3 S5 S8 S9 S11 S12 S60 S72 S105 D4 D5 D8 D9 D10 D12 D16 D28 D72'.split()]
L1, L2, L3, L4, L5, L6, L7, L8, S3, S5, S8, S9, S11, S12, S60, S72, S105, D4, D5, D8, D9, D10, D12, D16, D28, D72 = param
L = [p for p in param if str(p).startswith('L')]
S = [p for p in param if str(p).startswith('S')]
D = [p for p in param if str(p).startswith('D')]

s.add(And([And(p >= 0, p <= 1) for p in param]))
eqs = [L3 + L4 + S5 + S12 + L1 + D4 + L8 + S3 + L7 + D8 + D5 + L5,
       L3 + D12 + L1 + S9 + S3 + D5 + S105 + L2 + L7 + D28 + L5,
       L4 + D9 + S5 + L1 + D16 + L8 + L6 + S8 + L7 + D8,
       L4 + L1 + D16 + S60 + L2,
       S11 + L3 + S72 + D10 + D72 + D9 + S5 + D16 + S9 + S60 + L6 + S105 + L2 + L8 + D5,
       S72 + D72 + L8 + L6 + L5,
       L3 + S60 + L2 + L4,
       D72 + L6 + S105 + L7 + D28 + L5,
       D4 + S12 + S11 + D10,
       S11 + D12 + D10 + D72 + D9 + D28 + S72 + S5 + S12 + D4 + D16 + S9 + S3 + S60 + S8 + S105 + D8 + D5,
       D12 + D10 + S9 + S8 + D8 + S12,
       S11 + D12 + S72 + D9 + D4 + S3 + S8 + D28]
s.add(And([eq == 1 for eq in eqs]))
res = s.check()
count = 0
while res == sat:
    count += 1
    M = s.model()
    sol = [M[p] for p in param]
    if count % 100 == 0 or count <= 10:
        print(count, ", ".join([f"{p}:{M[p]}" for p in param]))
    s.add(Or([p != M[p] for p in param]))
    res = s.check()
print(f"{count} solutions found")

Output:

0 Solutions found