c-----*----------------------------------------------------6---------7--
c     This subroutine is to calculate the m-th numerical derivative of 
c       equally spaced data by use of an (n+1)-point formula. 
c       Here the maximum of m is 8. 
c       You can easily compile a subroutine for higher derivative than 8 through 
c           modifying the subroutine root(n,m,x,a).
c     input: nn,n,m,h,f
c       nn: (nn+1) is the number of sampling data, from 0 to nn
c       n: (n+1)-point formula,from 0 to n, n< or =nn
c       m: m-th derivative, m is less than or equal to n.
c       h: stepsize or sampling period
c       f(0:nn): function values at (nn+1) points
c     output: df
c       df(0:nn): approximations of the m-th derivative at (nn+1) sampling points 
c     work array: d
c       d(0:n,0:n): stores (n+1) coefficients of the (n+1)-point formula. 
c          In d(j,i), the index i=0,n denote (n+1) points xi (i=0,n), 
c          and the index j=0,n represent (n+1) coefficients at the 
c          conference point xi (i=0,n). 
c     subroutine: eqcoeff is to obtain the the coefficients d of equally
c       spaced data.
cccccc
c     Reference: 
c        Li, J., 2005: General explicit difference formulas for numerical differentiation. 
c                      J. Comput. Appl. Math., 183(1), 29-52, doi: 10.1016/j.cam.2004.12.026
cccccc
c     By Dr. Jianping Li, April 2, 2003.
      subroutine Ediffer(nn,n,m,h,f,df)
      dimension f(0:nn),df(0:nn)
      real d(0:n,0:n) !work array
      call eqcoeff(n,m,d)
      hm=h**m
      do i=0,nn
        df(i)=0.
      enddo
      k=int(n/2.)
      do i=0,k-1
        do j=0,n
          df(i)=df(i)+f(j)*d(j,i)
        enddo
      enddo
      do i=k,nn-k-1
        ik=i-k
        do j=0,n
          df(i)=df(i)+f(ik+j)*d(j,k)
        enddo
      enddo
      nk=nn-n
      do i=nn-k,nn
        ii=n-(nn-i)
        do j=0,n
          df(i)=df(i)+f(nk+j)*d(j,ii)
        enddo
      enddo
      do i=0,nn
        df(i)=df(i)/hm
      enddo
      return
      end

c-----*----------------------------------------------------6---------7--
c     This subroutine is used to calculate coefficients of an (n+1)-point 
c       formula for the m-th numerical derivative of equally spaced data.
c       Here the maximum of m is 8. 
c       One can easily calculate those coefficients for larger m
c       through inputting the common denominator array and 
c       writing the corresponding code in the subroutine ROOT.  
c     input: n,m
c       n: (n+1)-point formula
c       m: m-th derivative, m is less than or equal to n.
c     output: d
c       d(0:n,0:n): coefficient array,
c          In d(j,i), the index i=0,n denote (n+1) points xi (i=0,n), 
c          and the index j=0,n represent (n+1) coefficients of m-th numerical
c          derivative at the conference point xi (i=0,n). 
c     work array: x
c     work function: ifact
c       ifact(i):  external function, factorial of i 
c     By Dr. Jianping Li, April 2, 2003.
      subroutine eqcoeff(n,m,d)
      dimension d(0:n,0:n)
      dimension x(n)  !work array
      real ifact
	If(m.gt.n)then
	  write(*,*)'m must be less than or equal to n'
	  stop
	endif
      fm=ifact(m)
      do i=0,n
        do j=0,n
          d(j,i)=0.
        enddo
        do j=0,n
          if(j.ne.i)then
            kn=0
            do k=0,n
              if(k.ne.i.and.k.ne.j)then
                kn=kn+1
                x(kn)=k-i
              endif
            enddo
            call root(kn,m-1,x,a)

            d1=(-1)**(m-j)*fm*a
            f3=ifact(j)
            f4=ifact(n-j)
            d2=f3*f4
            d(j,i)=d1/d2
          endif
        enddo
        di=0.
        do j=0,n
          di=di+d(j,i) 
        enddo
        d(i,i)=-di
      enddo
      return
      end 
c-----*----------------------------------------------------6---------7--
c     Computing the factorial n! of a nutural number n.
      real function ifact(n)
      if(n.lt.0)then
        print*,'Illegal factorial',n
        stop
      else
        ifact=1.0
        do i=2,n
          ifact=ifact*float(i)
        enddo
      endif
      return
      end  
c-----*----------------------------------------------------6---------7--
c     This subroutine is to calculate the root coefficient of x1,...,xn 
c        (m)
c       An =Am(x1,x2,...,xn).
c     input: n,m,x
c       n: number of x
c       m: superscript index of An(m)
c       x(n): n numbers 
c     output: a
c       a: =An(m)=Am(x1,x2,...,xn)
c     By Dr. Jianping Li, April 2, 2003.
      subroutine root(n,m,x,a)
      dimension x(n)
      if(m.eq.0)then
        a=1.
        do i=1,n
          a=a*x(i)
        enddo
        return
      endif
      if(m.eq.1)then
        a=0.
        do j=1,n
          b=1.
          do i=1,n
            if(i.ne.j)then
              b=b*x(i)
            endif
          enddo
          a=a+b
        enddo
        return
      endif
      if(m.eq.2)then
        a=0.
        do k=1,n-1
          do j=k+1,n
            b=1.
            do i=1,n
              if(i.ne.j.and.i.ne.k)then
                b=b*x(i)
              endif
            enddo
            a=a+b
          enddo
        enddo
        return
      endif
      if(m.eq.3)then
        a=0.
        do l=1,n-2
          do k=l+1,n-1
            do j=k+1,n
              b=1.
              do i=1,n
                if(i.ne.j.and.i.ne.k.and.i.ne.l)then
                  b=b*x(i)
                  endif
              enddo
              a=a+b
            enddo
          enddo
        enddo
        return
      endif
      if(m.eq.4)then
        a=0.
        do 40 k1=   1,n-3
        do 40 k2=k1+1,n-2
        do 40 k3=k2+1,n-1
        do 40 k4=k3+1,n
          b=1.
          do i=1,n
            if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4)then
              b=b*x(i)
            endif
          enddo
          a=a+b
 40     continue
        return
      endif
      if(m.eq.5)then
        a=0.
        do 50 k1=   1,n-4
        do 50 k2=k1+1,n-3
        do 50 k3=k2+1,n-2
        do 50 k4=k3+1,n-1
        do 50 k5=k4+1,n
          b=1.
          do i=1,n
            if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and.
     *         i.ne.k5)then
              b=b*x(i)
            endif
          enddo
          a=a+b
 50     continue
        return
      endif
      if(m.eq.6)then
        a=0.
        do 60 k1=   1,n-5
        do 60 k2=k1+1,n-4
        do 60 k3=k2+1,n-3
        do 60 k4=k3+1,n-2
        do 60 k5=k4+1,n-1
        do 60 k6=k5+1,n
          b=1.
          do i=1,n
            if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and.
     *         i.ne.k5.and.i.ne.k6)then
              b=b*x(i)
            endif
          enddo
          a=a+b
 60     continue
        return
      endif
      if(m.eq.7)then
        a=0.
        do 70 k1=   1,n-6
        do 70 k2=k1+1,n-5
        do 70 k3=k2+1,n-4
        do 70 k4=k3+1,n-3
        do 70 k5=k4+1,n-2
        do 70 k6=k5+1,n-1
        do 70 k7=k6+1,n
          b=1.
          do i=1,n
            if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and.
     *         i.ne.k5.and.i.ne.k6.and.i.ne.k7)then
              b=b*x(i)
            endif
          enddo
          a=a+b
 70     continue
        return
      endif
      if(m.eq.8)then
        a=0.
        do 80 k1=   1,n-7
        do 80 k2=k1+1,n-6
        do 80 k3=k2+1,n-5
        do 80 k4=k3+1,n-4
        do 80 k5=k4+1,n-3
        do 80 k6=k5+1,n-2
        do 80 k7=k6+1,n-1
        do 80 k8=k7+1,n
          b=1.
          do i=1,n
            if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and.
     *         i.ne.k5.and.i.ne.k6.and.i.ne.k7.and.i.ne.k8)then
              b=b*x(i)
            endif
          enddo
          a=a+b
 80     continue
        return
      endif
      return
      end