code: def bilinear(x, W, y, input_size, seq_len, batch_size, num_outputs=1, bias_x=False, bias_y=False):
      
      if bias_x:
          x = nd.concat(x, nd.ones((1, seq_len, batch_size)), dim=0)
      if bias_y:
          y = nd.concat(y, nd.ones((1, seq_len, batch_size)), dim=0)
  
      nx, ny = input_size + bias_x, input_size + bias_y
      # W: (num_outputs x ny) x nx
      lin = nd.dot(W, x)
      if num_outputs > 1:
          lin = reshape_fortran(lin, (ny, num_outputs * seq_len, batch_size))
      y = y.transpose([2, 1, 0])  # May cause performance issues
      lin = lin.transpose([2, 1, 0])
      blin = nd.batch_dot(lin, y, transpose_b=True)
      blin = blin.transpose([2, 1, 0])
      if num_outputs > 1:
          blin = reshape_fortran(blin, (seq_len, num_outputs, seq_len, batch_size))
      return blin
docs: Do xWy
  
      Parameters
      ----------
      x : NDArray
          (input_size x seq_len) x batch_size
      W : NDArray
          (num_outputs x ny) x nx
      y : NDArray
          (input_size x seq_len) x batch_size
      input_size : int
          input dimension
      seq_len : int
          sequence length
      batch_size : int
          batch size
      num_outputs : int
          number of outputs
      bias_x : bool
          whether concat bias vector to input x
      bias_y : bool
          whether concat bias vector to input y
  
      Returns
      -------
      output : NDArray
          [seq_len_y x seq_len_x if output_size == 1 else seq_len_y x num_outputs x seq_len_x] x batch_size