docs: def _qnwgamma1(n, a=1.0, b=1.0, tol=3e-14):
      """
      1d quadrature weights and nodes for Gamma distributed random variable
  
      Parameters
      ----------
      n : scalar : int
          The number of quadrature points
  
      a : scalar : float, optional(default=1.0)
          Shape parameter of the gamma distribution parameter. Must be positive
  
      b : scalar : float, optional(default=1.0)
          Scale parameter of the gamma distribution parameter. Must be positive
  
      tol : scalar : float, optional(default=3e-14)
          Tolerance parameter for newton iterations for each node
  
      Returns
      -------
      nodes : np.ndarray(dtype=float, ndim=1)
          The quadrature points
  
      weights : np.ndarray(dtype=float, ndim=1)
          The quadrature weights that correspond to nodes
  
      Notes
      -----
      Based of original function ``qnwgamma1`` in CompEcon toolbox by
      Miranda and Fackler
  
      References
      ----------
      Miranda, Mario J, and Paul L Fackler. Applied Computational
      Economics and Finance, MIT Press, 2002.
  
      """
      a -= 1
  
      maxit = 25
  
      factor = -math.exp(gammaln(a+n) - gammaln(n) - gammaln(a+1))
      nodes = np.zeros(n)
      weights = np.zeros(n)
  
      # Create nodes
      for i in range(n):
          # Reasonable starting values
          if i == 0:
              z = (1+a) * (3+0.92*a) / (1 + 2.4*n + 1.8*a)
          elif i == 1:
              z = z + (15 + 6.25*a) / (1 + 0.9*a + 2.5*n)
          else:
              j = i-1
              z = z + ((1 + 2.55*j) / (1.9*j) + 1.26*j*a / (1 + 3.5*j)) * \
                  (z - nodes[j-1]) / (1 + 0.3*a)
  
          # root finding iterations
          its = 0
          z1 = -10000
          while abs(z - z1) > tol and its < maxit:
              p1 = 1.0
              p2 = 0.0
              for j in range(1, n+1):
                  # Recurrance relation for Laguerre polynomials
                  p3 = p2
                  p2 = p1
                  p1 = ((2*j - 1 + a - z)*p2 - (j - 1 + a)*p3) / j
  
              pp = (n*p1 - (n+a)*p2) / z
              z1 = z
              z = z1 - p1/pp
              its += 1
  
          if its == maxit:
              raise ValueError('Failure to converge')
  
          nodes[i] = z
          weights[i] = factor / (pp*n*p2)
  
      return nodes*b, weights