<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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