code: def construct_zernike_lookuptable(zernike_indexes):
      
      n_max = np.max(zernike_indexes[:,0])
      factorial = np.ones((1 + n_max,), dtype=float)
      factorial[1:] = np.cumproduct(np.arange(1, 1 + n_max, dtype=float))
      width = int(n_max//2 + 1)
      lut = np.zeros((zernike_indexes.shape[0],width), dtype=float)
      for idx, (n, m) in enumerate(zernike_indexes):
          alt = 1
          npmh = (n+m)//2
          nmmh = (n-m)//2
          for k in range(0,nmmh+1):
              lut[idx,k] = \
                  (alt * factorial[n-k] /
                   (factorial[k]*factorial[npmh-k]*factorial[nmmh-k]))
              alt = -alt
      return lut
docs: Return a lookup table of the sum-of-factorial part of the radial
      polynomial of the zernike indexes passed
      
      zernike_indexes - an Nx2 array of the Zernike polynomials to be
                        computed.