docs: def cy_bispev(tx, ty, c, kx, ky, x, y):
      '''Possible optimization: Do not evaluate derivatives, ever.
      '''
      nx = len(tx)
      ny = len(ty)
      mx = len(x)
      my = len(y)
  
      kx1 = kx + 1
      ky1 = ky + 1
  
      nkx1 = nx - kx1
      nky1 = ny - ky1
      
      wx = [[0.0]*kx1]*mx
      wy = [[0.0]*ky1]*my
      lx = [0]*mx
      ly = [0]*my
  
      size_z = mx*my
      
      z = [0.0]*size_z
      init_w(tx, kx, x, lx, wx)
      init_w(ty, ky, y, ly, wy)
      
      for j in range(my):
          for i in range(mx):
              sp = 0.0
              err = 0.0
              for i1 in range(kx1):
                  for j1 in range(ky1):
                      l2 = lx[i]*nky1 + ly[j] + i1*nky1 + j1
                      a = c[l2]*wx[i][i1]*wy[j][j1] - err
                      tmp = sp + a
                      err = (tmp - sp) - a
                      sp = tmp
              z[j*mx + i] += sp
      return z