def det(matrix):
    """Seem sto work fine.
    
    >> from sympy import *
    >> from sympy.abc import *
    >> Matrix([[a, b], [c, d]]).det()
    a*d - b*c
    >> Matrix([[a, b, c], [d, e, f], [g, h, i]]).det()
    a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g
    
    A few terms can be slightly factored out of the 3x dim.
    
    >> Matrix([[a, b, c, d], [e, f, g, h], [i, j, k, l], [m, n, o, p]]).det()
    a*f*k*p - a*f*l*o - a*g*j*p + a*g*l*n + a*h*j*o - a*h*k*n - b*e*k*p + b*e*l*o + b*g*i*p - b*g*l*m - b*h*i*o + b*h*k*m + c*e*j*p - c*e*l*n - c*f*i*p + c*f*l*m + c*h*i*n - c*h*j*m - d*e*j*o + d*e*k*n + d*f*i*o - d*f*k*m - d*g*i*n + d*g*j*m

    72 mult vs ~48 in cse'd version'
    
    Commented out - takes a few seconds
    >> #Matrix([[a, b, c, d, e], [f, g, h, i, j], [k, l, m, n, o], [p, q, r, s, t], [u, v, w, x, y]]).det()
    
    260 multiplies with cse; 480 without it.
    """
    size = len(matrix)
    if size == 1:
        return matrix[0]
    elif size == 2:
        (a, b), (c, d) = matrix
        return a * d - c * b
    elif size == 3:
        (a, b, c), (d, e, f), (g, h, i) = matrix
        return a * (e * i - h * f) - d * (b * i - h * c) + g * (b * f - e * c)
    elif size == 4:
        (a, b, c, d), (e, f, g, h), (i, j, k, l), (m, n, o, p) = matrix
        return (a * f * k * p - a * f * l * o - a * g * j * p + a * g * l *
            n + a * h * j * o - a * h * k * n - b * e * k * p + b * e * l *
            o + b * g * i * p - b * g * l * m - b * h * i * o + b * h * k *
            m + c * e * j * p - c * e * l * n - c * f * i * p + c * f * l *
            m + c * h * i * n - c * h * j * m - d * e * j * o + d * e * k *
            n + d * f * i * o - d * f * k * m - d * g * i * n + d * g * j * m)
    elif size == 5:
        (a, b, c, d, e), (f, g, h, i, j), (k, l, m, n, o), (p, q, r, s, t), (u,
            v, w, x, y) = matrix
        x0 = s * y
        x1 = a * g * m
        x2 = t * w
        x3 = a * g * n
        x4 = r * x
        x5 = a * g * o
        x6 = t * x
        x7 = a * h * l
        x8 = q * y
        x9 = a * h * n
        x10 = s * v
        x11 = a * h * o
        x12 = r * y
        x13 = a * i * l
        x14 = t * v
        x15 = a * i * m
        x16 = q * w
        x17 = a * i * o
        x18 = s * w
        x19 = a * j * l
        x20 = q * x
        x21 = a * j * m
        x22 = r * v
        x23 = a * j * n
        x24 = b * f * m
        x25 = b * f * n
        x26 = b * f * o
        x27 = b * h * k
        x28 = t * u
        x29 = b * h * n
        x30 = p * x
        x31 = b * h * o
        x32 = b * i * k
        x33 = p * y
        x34 = b * i * m
        x35 = r * u
        x36 = b * i * o
        x37 = b * j * k
        x38 = s * u
        x39 = b * j * m
        x40 = p * w
        x41 = b * j * n
        x42 = c * f * l
        x43 = c * f * n
        x44 = c * f * o
        x45 = c * g * k
        x46 = c * g * n
        x47 = c * g * o
        x48 = c * i * k
        x49 = c * i * l
        x50 = p * v
        x51 = c * i * o
        x52 = c * j * k
        x53 = c * j * l
        x54 = q * u
        x55 = c * j * n
        x56 = d * f * l
        x57 = d * f * m
        x58 = d * f * o
        x59 = d * g * k
        x60 = d * g * m
        x61 = d * g * o
        x62 = d * h * k
        x63 = d * h * l
        x64 = d * h * o
        x65 = d * j * k
        x66 = d * j * l
        x67 = d * j * m
        x68 = e * f * l
        x69 = e * f * m
        x70 = e * f * n
        x71 = e * g * k
        x72 = e * g * m
        x73 = e * g * n
        x74 = e * h * k
        x75 = e * h * l
        x76 = e * h * n
        x77 = e * i * k
        x78 = e * i * l
        x79 = e * i * m
        return (x0 * x1 - x0 * x24 + x0 * x27 + x0 * x42 - x0 * x45 - x0 *
            x7 - x1 * x6 + x10 * x11 - x10 * x21 - x10 * x44 + x10 * x52 + 
            x10 * x69 - x10 * x74 - x11 * x20 + x12 * x13 + x12 * x25 - x12 *
            x3 - x12 * x32 - x12 * x56 + x12 * x59 - x13 * x2 + x14 * x15 +
            x14 * x43 - x14 * x48 - x14 * x57 + x14 * x62 - x14 * x9 - x15 *
            x8 + x16 * x17 - x16 * x23 - x16 * x58 + x16 * x65 + x16 * x70 -
            x16 * x77 - x17 * x22 + x18 * x19 + x18 * x26 - x18 * x37 - x18 *
            x5 - x18 * x68 + x18 * x71 - x19 * x4 - x2 * x25 + x2 * x3 + x2 *
            x32 + x2 * x56 - x2 * x59 + x20 * x21 + x20 * x44 - x20 * x52 -
            x20 * x69 + x20 * x74 + x22 * x23 + x22 * x58 - x22 * x65 - x22 *
            x70 + x22 * x77 + x24 * x6 - x26 * x4 - x27 * x6 + x28 * x29 - 
            x28 * x34 - x28 * x46 + x28 * x49 + x28 * x60 - x28 * x63 - x29 *
            x33 + x30 * x31 - x30 * x39 - x30 * x47 + x30 * x53 + x30 * x72 -
            x30 * x75 - x31 * x38 + x33 * x34 + x33 * x46 - x33 * x49 - x33 *
            x60 + x33 * x63 + x35 * x36 - x35 * x41 - x35 * x61 + x35 * x66 +
            x35 * x73 - x35 * x78 - x36 * x40 + x37 * x4 + x38 * x39 + x38 *
            x47 - x38 * x53 - x38 * x72 + x38 * x75 + x4 * x5 + x4 * x68 - 
            x4 * x71 + x40 * x41 + x40 * x61 - x40 * x66 - x40 * x73 + x40 *
            x78 - x42 * x6 - x43 * x8 + x45 * x6 + x48 * x8 + x50 * x51 - 
            x50 * x55 - x50 * x64 + x50 * x67 + x50 * x76 - x50 * x79 - x51 *
            x54 + x54 * x55 + x54 * x64 - x54 * x67 - x54 * x76 + x54 * x79 +
            x57 * x8 + x6 * x7 - x62 * x8 + x8 * x9)
    else:
        import numpy as np
        return float(np.linalg.det(matrix))
