docs: 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:
          # TODO algorithm?
          import numpy as np
          return float(np.linalg.det(matrix))