docs: def qnwequi(n, a, b, kind="N", equidist_pp=None, random_state=None):
      """
      Generates equidistributed sequences with property that averages
      value of integrable function evaluated over the sequence converges
      to the integral as n goes to infinity.
  
      Parameters
      ----------
      n : int
          Number of sequence points
  
      a : scalar or array_like(float)
          A length-d iterable of lower endpoints. If a scalar is given,
          that constant is repeated d times, where d is the number of
          dimensions
  
      b : scalar or array_like(float)
          A length-d iterable of upper endpoints. If a scalar is given,
          that constant is repeated d times, where d is the number of
          dimensions
  
      kind : string, optional(default="N")
          One of the following:
  
          - N - Neiderreiter (default)
          - W - Weyl
          - H - Haber
          - R - pseudo Random
  
      equidist_pp : array_like, optional(default=None)
          TODO: I don't know what this does
  
      random_state : int or np.random.RandomState, optional
          Random seed (integer) or np.random.RandomState instance to set
          the initial state of the random number generator for
          reproducibility. If None, a randomly initialized RandomState is
          used.
  
      Returns
      -------
      nodes : np.ndarray(dtype=float)
          Quadrature nodes
  
      weights : np.ndarray(dtype=float)
          Weights for quadrature nodes
  
      Notes
      -----
      Based of original function ``qnwequi`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      random_state = check_random_state(random_state)
  
      if equidist_pp is None:
          import sympy as sym
          equidist_pp = np.sqrt(np.array(list(sym.primerange(0, 7920))))
  
      n, a, b = list(map(np.atleast_1d, list(map(np.asarray, [n, a, b]))))
  
      d = max(list(map(len, [n, a, b])))
      n = np.prod(n)
  
      if a.size == 1:
          a = np.repeat(a, d)
  
      if b.size == 1:
          b = np.repeat(b, d)
  
      i = np.arange(1, n + 1)
  
      if kind.upper() == "N":  # Neiderreiter
          j = 2.0 ** (np.arange(1, d+1) / (d+1))
          nodes = np.outer(i, j)
          nodes = (nodes - fix(nodes)).squeeze()
      elif kind.upper() == "W":  # Weyl
          j = equidist_pp[:d]
          nodes = np.outer(i, j)
          nodes = (nodes - fix(nodes)).squeeze()
      elif kind.upper() == "H":  # Haber
          j = equidist_pp[:d]
          nodes = np.outer(i * (i+1) / 2, j)
          nodes = (nodes - fix(nodes)).squeeze()
      elif kind.upper() == "R":  # pseudo-random
          nodes = random_state.rand(n, d).squeeze()
      else:
          raise ValueError("Unknown sequence requested")
  
      # compute nodes and weights
      r = b - a
      nodes = a + nodes * r
      weights = (np.prod(r) / n) * np.ones(n)
  
      return nodes, weights