subroutine glrule(n,x,w) 

! grule calculates the nodes x and weights w for the n-point gauss rule.
! DJI's version of grule from Davis & Rabinowitz, "Methods of Numerical Integration". 
! This routine is more transparent but requires (3JJ+2)M + 2A extra operations, where
! JJ is the number of j iterations. This is not a problem, since the weights and
! nodes need only be calculated once.
! Nodes/weights should be accurate to 15/13 significant figures.

double precision pkm1,pk,t1,pkp1,u,v,h,fx,x0
double precision den,d1,dpn,d2pn,d3pn,d4pn,p,dp,d2p,d3p
double precision x(n),w(n)

m=(n+1)/2
e1=n*(n+1)

do i=1,m
   t=(4*i-1)*3.141592653589793/(4*n+2)
   x0=(1.-(1.-1./n)/(8.*n*n))*cos(t)

   do j=1,2   
      ! Calculate Pn(x0) from recurrence relation 1:
      pkm1=1.
      pk=x0
      
      do k=2,n
         t1=x0*pk
         pkp1=t1-pkm1-(t1-pkm1)/k+t1
         pkm1=pk
         pk=pkp1
      enddo

      ! Calculate first derivative of Pn at x0 from recurrence relation 2
      ! and second, third, fourth derivatives of Pn at x0 from d.e. and 
      ! its derivatives:
      den=1.-x0*x0
      d1=n*(pkm1-x0*pk)
      dpn=d1/den
      d2pn=(2.*x0*dpn-e1*pk)/den
      d3pn=(4.*x0*d2pn+(2.-e1)*dpn)/den
      d4pn=(6.*x0*d3pn+(6.-e1)*d2pn)/den

      ! First iteration of root-finder:
      u=pk/dpn
      v=d2pn/dpn
      h=-u*(1.+.5*u*v)
      x0=x0+h
      den=1.-x0*x0

      ! Second iteration of root-finder: calculating Pn and its derivatives 
      ! about new iterate from 5-term Taylor expansion (avoids having to
      ! run through recurrence relation 1 all over again):
      p=pk+h*(dpn+.5*h*(d2pn+h/3.*(d3pn+.25*h*d4pn)))
      dp=dpn+h*(d2pn+.5*h*(d3pn+h/3.*d4pn))
      d2p=d2pn+h*(d3pn+.5*h*d4pn)
      u=p/dp
      v=d2p/dp
      h=-u*(1.+.5*u*v)
      x0=x0+h 
    
      ! Must use recurrence relation 1 here since accuracy of Taylor expansion
      ! exhausted. Can trade off iteration of recurrence relation 1 with more 
      ! terms in Taylor expansion - worthwhile for very large n:
   enddo   

   x(i)=x0
   x(n-i+1)=-x(i) 

   ! Taylor expansion to evaluate (1-x^2)dPn/dx about previous iterate.
   ! The d.e. can be used to evaluate the derivatives in terms of Pn and
   ! its derivatives:
 
   d3p=d3pn+h*d4pn 
   fx=den*dp-h*e1*(p+.5*h*(dp+h/3.*(d2p+.25*h*(d3p+.2*h*d4pn))))
   w(i)=2.*(1.-x(i)*x(i))/(fx*fx) 

   ! Alternatively (and less accurately), calculate Pn-1 using recurrence 
   ! relation 1:
   !pkm1=1.
   !pk=x0
   !do k=2,n
      !t1=x0*pk
      !pkp1=t1-pkm1-(t1-pkm1)/k+t1
      !pkm1=pk
      !pk=pkp1
   !enddo
   !w(i)=2.*(1.-x(i)*x(i))/(n*pkm1)**2

   w(n-i+1)=w(i)
enddo

if(m+m.gt.n) x(m)=0.

end subroutine glrule