      subroutine un_init(nn) ! bernd berg may 29 2001.
c initialize coefficients to calculate the probability distribution
c of the sums of n uniform distributions.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      common/uniform_n/ a(0:nmax,0:nmax),u_df(nmax,nmax),n,nm1
      dimension b(0:nmax,0:nmax) ! recursion from b to a.
      ltest=.true.
      ltest=.false.
      if(ltest) stop "un_init: ltest."
      if(nn.lt.1 .or. nn.gt.nmax) then
        write(6,'(" un_init: nn out of range! nn,nmax =",2i8)') nn,nmax
        stop "un_init: n out of range."
      end if

      write(6,'(/," un_init: recursion up to n =",i8,/)') nn
      write(6,'(a,/)')
     &  ' n, um, u2m, sigma(u), e1sigma, e2sigma, e3sigma, e4sigma:'
      do k=0,nn
      do i=0,nn
        b(i,k)=zero
        a(i,k)=zero
      end do
      end do
      b(0,1)=one

      do n=2,nn
        nm1=n-1

        do k=1,n ! the constant-term coefficient.
          a(0,k)=zero
          do i=0,nm1
            a(0,k)=a(0,k)+b(i,k-1)/(i+1)
          end do
        end do

        do k=1,n ! the higher order x**i coefficients.
        do i=0,nm1
          a(i+1,k)=(b(i,k)-b(i,k-1))/(i+1)
        end do
        end do

        do k=1,n ! prepare next recursion step.
        do i=0,nm1
          b(i,k)=a(i,k)
        end do
        end do

        do k=1,n
          un_add=a(nm1,k)/n
          do i=nm1-1,0,-1
            un_add=un_add+a(i,k)/(i+1)
          end do
          u_df(k,n)=un_add
        end do

c calculation of the first (for check result n/2) and second moment (unknown):
        um=zero
        u2m=zero
        do k=1,n
          xkm1=(k-1)*one
          txkm1=two*xkm1
          xkm12=xkm1**2
          um=um+xkm1*a(nm1,k)/n
          um=um+a(nm1,k)/(n+1)
          u2m=u2m+xkm12*a(nm1,k)/n
          u2m=u2m+txkm1*a(nm1,k)/(n+1)
          u2m=u2m+a(nm1,k)/(n+2)
          do i=nm1-1,0,-1
            um=um+xkm1*a(i,k)/(i+1)
            um=um+a(i,k)/(i+2)
            u2m=u2m+xkm12*a(i,k)/(i+1)
            u2m=u2m+txkm1*a(i,k)/(i+2)
            u2m=u2m+a(i,k)/(i+3)
          end do
        end do
        sigma=sqrt(u2m-um**2)
        u1=half*n-sigma
        e1sigma=two*un_df(u1)
        u1=half*n-two*sigma
        e2sigma=zero
        if(u1.gt.zero) e2sigma=two*un_df(u1)
        u1=half*n-three*sigma
        e3sigma=zero
        if(u1.gt.zero) e3sigma=two*un_df(u1)
        u1=half*n-four*sigma
        e4sigma=zero
        if(u1.gt.zero) e4sigma=two*un_df(u1)
        write(6,'(i4,f7.2,f9.3,2f8.3,f9.4,f10.5,g12.3)') 
     &    n,um,u2m,sigma,e1sigma,e2sigma,e3sigma,e4sigma

      end do

      n=nn
      write(6,'(/," k,a(i,k) for n =",i8)') n
      do k=1,n
        write(6,'(/," k =",i8)') k
        write(6,'(6g12.4)') (a(i,k),i=0,nm1)
        asum=zero
        do i=0,nm1
          asum=asum+a(i,k)
        end do
        if(k.lt.n) write(6,'(" asum,a(0,k+1):",2g22.12)') asum,a(0,k+1)
      end do

      return
      end


      function un_test(u) ! bernd berg, may 29 2001.
c calculates the probability distribution of the sums of 
c n uniform distributions.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      common/uniform_n/ a(0:nmax,0:nmax),u_df(nmax,nmax),n,nm1
      if(n.gt.nmax) stop "un_test: enlarge nmax!"
      k=1+int(u)
      if(k.lt.1) k=1+int(u+em12)
      if(k.gt.n) k=1+int(u-em12)
      if(k.lt.1 .or. k.gt.n) then
        write(*,'(" un_test: k out of range!")')
        write(*,'("          u,k,n =",g16.7,2i8)') u,k,n
        stop "un_test: k out of range."
      end if
      x=u-k+1
      un_test=zero
      do i=0,nm1
        un_test=un_test+a(i,k)*x**i
      end do
      return
      end


      function un_pd(u) ! bernd berg, june 5 2001.
c calculates the probability density for the sum of n uniform distributions.
c call un_init first to set the parameters of the common block /uniform_n/.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      common/uniform_n/ a(0:nmax,0:nmax),u_df(nmax,nmax),n,nm1
      if(n.gt.nmax) stop "un: enlarge nmax!"
      uu=u
      if(u.gt.half*n) uu=one*n-u
      k=1+int(uu)
      if(k.lt.1) k=1+int(uu+em12)
      if(k.lt.1) then
        write(6,'(" un: k out of range! u,uu,k,n =",2g16.7,2i8)') u,uu,k,n
        stop "un: k out of range!"
      end if
      x=uu-k+1
      un_pd=a(nm1,k)*x
      do i=nm1-1,1,-1
        un_pd=(un_pd+a(i,k))*x
      end do
      un_pd=un_pd+a(0,k)
      return
      end


      function un_df(u) ! bernd berg, june 5 2001.
c calculates the cumulative distribution function for the sum of n uniform 
c distributions. call un_init first to set the parameters of the common block.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      common/uniform_n/ a(0:nmax,0:nmax),u_df(nmax,nmax),n,nm1
      if(n.gt.nmax) stop "un: enlarge nmax!"
      uu=u
      if(u.gt.half*n) uu=one*n-u
      k=1+int(uu)
      if(k.lt.1) k=1+int(uu+em12)
      if(k.lt.1) then
        write(6,'(" un: k out of range! u,uu,k,n =",2g16.7,2i8)') u,uu,k,n
        stop "un: k out of range!"
      end if
      x=uu-k+1
      un_df=(a(nm1,k)/n)*x
      do i=nm1-1,0,-1
        un_df=(un_df+a(i,k)/(i+1))*x
      end do
      km1=k-1
      do k=1,km1
        un_df=un_df+u_df(k,n)
      end do
      if(u.gt.half*n) un_df=one-un_df
      return
      end


      function un_dfln(u) ! bernd berg, june 5 2001.
c calculates the cumulative distribution function for the sum of n uniform 
c distributions. call un_init first to set the parameters of the common block.
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(nmax=100)
      common/uniform_n/ a(0:nmax,0:nmax),u_df(nmax,nmax),n,nm1
      ltest=.false.
      if(n.gt.nmax) stop "un: enlarge nmax!"
      if(u.ge.one*n) then
        if(u.gt.one*n) write(6,'(" warning: u =",1g22.12)') u
        un_dfln=zero
        return
      end if
      uu=u
      if(u.gt.half*n) uu=one*n-u
      k=1+int(uu)
      if(k.lt.1) k=1+int(uu+em12)
      if(k.lt.1) then
        write(6,'(" un: k out of range! u,uu,k,n =",2g16.7,2i8)') u,uu,k,n
        stop "un: k out of range!"
      end if
      x=uu-k+one
      if(x.lt.em12) then
        un_dfln=-one/em12
      else
        un_dfln=log(a(nm1,k)/n)+log(x)
        do i=nm1-1,0,-1
          a_ln=log(a(i,k)/(i+1))
          un_dfln=addln(un_dfln,a_ln)+log(x)
        end do
      end if
      if(ltest) write(11,*) "k,uu,x,un_dfln:",k,uu,x,un_dfln
      km1=k-1
      do k=1,km1
        a_ln=log(u_df(k,n))
        un_dfln=addln(un_dfln,a_ln)
      end do
      if(ltest) write(11,*) "u,un_dfln =",u,un_dfln
      return
      end
