︠4f58aad8-aa03-4417-90a5-ffaec91c96da︠
### matrix
︡dc8821bb-9e64-4287-9ece-9cdb7c24831f︡{"done":true}︡
︠860348a5-7c41-4936-a49f-29c3fcf67997s︠
### matrix is not a standard object in Python
### but is standard in Sage.

### Use matrix( list of list ) to construct a matrix.
M = matrix([
[1,2,3],
[4,5,6],
[7,8,9]
])
print M
show(M)
︡daf6e6d7-dd42-477e-8fd2-a88e73ed4c9f︡{"stdout":"[1 2 3]\n[4 5 6]\n[7 8 9]\n"}︡{"html":"<div align='center'>$\\displaystyle \\left(\\begin{array}{rrr}\n1 &amp; 2 &amp; 3 \\\\\n4 &amp; 5 &amp; 6 \\\\\n7 &amp; 8 &amp; 9\n\\end{array}\\right)$</div>"}︡{"done":true}︡
︠334ac53b-c2c4-4c28-83b9-0728882dc978s︠
### Use M[i,j] to obtain the i,j-entry of M

M[1,0]
︡e4e65372-b84a-4752-a18d-5f981c07c948︡{"stdout":"4\n"}︡{"done":true}︡
︠cc3c6be6-4006-4a18-bb72-d5a72c8323ffs︠
### You may assign value to an entry by =

M[1,0] = -1
M
︡4540574b-61fe-486e-91bc-6e86ca367190︡{"stdout":"[ 1  2  3]\n[-1  5  6]\n[ 7  8 10]\n"}︡{"done":true}︡
︠4d0504f5-a149-4624-9f24-5698895ded0fs︠
### Alternatively, create a zero matrix and change the values of entries

M = zero_matrix(3,3)
for i in range(3):
    for j in range(3):
        M[i,j] = i*j
M
︡d752d3df-c6e6-4f92-bc63-34d9d4fbaf73︡{"stdout":"[0 0 0]\n[0 1 2]\n[0 2 4]\n"}︡{"done":true}︡
︠ae4adacb-33dc-4fe2-8004-7ba467901fdds︠
### Or input the number of rows and a list.

M = matrix(3,range(9))
M
︡a87f9f5b-8a53-4177-bf32-06bf8b5b07d8︡{"stdout":"[0 1 2]\n[3 4 5]\n[6 7 8]\n"}︡{"done":true}︡
︠6d0cd56d-5e9c-4a55-8c2b-4b48f48eea22︠
### invertible or not

M.is_invertible()
︡27f0172e-ee1e-47ba-bfd5-ad664fade489︡{"stdout":"False\n"}︡{"done":true}︡
︠76c2092a-55e3-4c7a-a32e-02dc138f2859s︠
### If invertible, M.inverse() gives the inverse.

M = matrix([
[1,2,3],
[4,5,6],
[7,8,10]
])

M.inverse()
︡17f13e59-d931-4bd9-a456-fca1cdf29127︡{"stdout":"[-2/3 -4/3    1]\n[-2/3 11/3   -2]\n[   1   -2    1]\n"}︡{"done":true}︡
︠8fbfd124-f814-461f-9282-3a0d5176090bs︠
### Taking a submatrix

### When rows == [0,1,2], rows[0:2] == [0,1]
### When cols == [0,1,2], cols[1:] == [1,2]

M[0:2,1:]

︡97d26b50-3da5-4517-bb3e-22434b956818︡{"stdout":"[2 3]\n[5 6]\n"}︡{"done":true}︡
︠3803f5a1-1cac-4bff-a5d1-90c03fc402d5s︠
### Create a matrix for the poset Dn

n = 50
Dn = [k for k in range(1,51) if 50 % k == 0]
size = len(Dn)
print Dn
︡81466c22-538e-4a88-94c4-2c0f10420ffd︡{"stdout":"[1, 2, 5, 10, 25, 50]\n"}︡{"done":true}︡
︠4f12417b-33c5-47d4-b8ec-3c95aff79668s︠
### Create the mappings between Dn and index

Dn_to_ind = {Dn[i]: i for i in range(size)}
print Dn_to_ind
ind_to_Dn = {i: Dn[i] for i in range(size)}
print ind_to_Dn
︡d26520d5-19e4-4659-9ee6-c4631843d3fc︡{"stdout":"{1: 0, 2: 1, 5: 2, 10: 3, 50: 5, 25: 4}\n"}︡{"stdout":"{0: 1, 1: 2, 2: 5, 3: 10, 4: 25, 5: 50}\n"}︡{"done":true}︡
︠9a8151a1-c22b-4abd-bff9-4d06407b5fe0s︠
### Create the matrix for the cover function

M = zero_matrix(size,size)
for i in range(size):
    for j in range(size):
        a = ind_to_Dn[i]
        b = ind_to_Dn[j]
        if b % a == 0 and Integer(b/a).is_prime(): ### This means b covers a
            M[i,j] = 1

M
︡df2429e3-6d59-433c-ae5e-80bf840ddd53︡{"stdout":"[0 1 1 0 0 0]\n[0 0 0 1 0 0]\n[0 0 0 1 1 0]\n[0 0 0 0 0 1]\n[0 0 0 0 0 1]\n[0 0 0 0 0 0]\n"}︡{"done":true}︡
︠5253af92-1fa9-409c-8066-1290676183d3︠
︡b0ec7fcb-1f57-4f71-9ebc-d48d28ff9588︡
︠9f20d4a5-6626-47fd-9308-d8da5f58503di︠
%md

## Project 1

Define two functions `zeta_func(n)` and `moebius_func(n)`.
`zeta_func(n)` returns the matrix of the zeta function on $D_n$.
`moebius_func(n)` returns the matrix of the Moebius function on $D_n$.
You may try if the two matrices are the inverse of each other.

︡1491331a-ff0e-4231-ae95-9c18acf6be1f︡{"done":true,"md":"\n## Project 1\n\nDefine two functions `zeta_func(n)` and `moebius_func(n)`.  \n`zeta_func(n)` returns the matrix of the zeta function on $D_n$.  \n`moebius_func(n)` returns the matrix of the Moebius function on $D_n$.  \nYou may try if the two matrices are the inverse of each other."}
︠49a8da51-288f-4181-b23f-c90972e0f3f6︠
def zeta_func(n):
    ### Your answer below
    Dn = [k for k in range(1,51) if 50 % k == 0]
    size = len(Dn)
    Dn_to_ind = {Dn[i]: i for i in range(size)}
    ind_to_Dn = {i: Dn[i] for i in range(size)}
    M = zero_matrix(size,size)
    for i in range(size):
        for j in range(size):
            a = ind_to_Dn[i]
            b = ind_to_Dn[j]
            if b % a == 0:
                M[i,j] = 1
    return M

def moebius_func(n):
    ### Your answer below
    ### Bonus:  Build the matrix by the definition of the Moebius function
    ###         ( instead of returnning zeta_func(n).inverse() )
    Dn = [k for k in range(1,51) if 50 % k == 0]
    size = len(Dn)
    Dn_to_ind = {Dn[i]: i for i in range(size)}
    ind_to_Dn = {i: Dn[i] for i in range(size)}
    M = zero_matrix(size,size)
    for i in range(size):
        for j in range(size):
            if i == j:
                M[i,j] = 1
            else:
                a = ind_to_Dn[i]
                b = ind_to_Dn[j]
                if b % a == 0:
                    mu_ij = -sum([M[i,Dn_to_ind[c]] for c in Dn if c % a == 0 and b % c == 0 and c != b])
                    M[i,j] = mu_ij
    return M

n = 8
print "n =", n
print "zeta"
print zeta_func(n)
print "moebius"
print moebius_func(n)

print "product"
print zeta_func(n) * moebius_func(n)
︡ce55ac9a-2957-4651-93f8-0237009ab68b︡{"stdout":"n = 8\n"}︡{"stdout":"zeta\n"}︡{"stdout":"[1 1 1 1 1 1]\n[0 1 0 1 0 1]\n[0 0 1 1 1 1]\n[0 0 0 1 0 1]\n[0 0 0 0 1 1]\n[0 0 0 0 0 1]\n"}︡{"stdout":"moebius\n"}︡{"stdout":"[ 1 -1 -1  1  0  0]\n[ 0  1  0 -1  0  0]\n[ 0  0  1 -1 -1  1]\n[ 0  0  0  1  0 -1]\n[ 0  0  0  0  1 -1]\n[ 0  0  0  0  0  1]\n"}︡{"stdout":"product\n"}︡{"stdout":"[1 0 0 0 0 0]\n[0 1 0 0 0 0]\n[0 0 1 0 0 0]\n[0 0 0 1 0 0]\n[0 0 0 0 1 0]\n[0 0 0 0 0 1]\n"}︡{"done":true}︡
︠1635e1c3-02e1-45a0-9728-3837b837ee80︠