function p(x)   ! 1 <= x < infinity ! October 31/2004/Dubna-Sofia
       implicit real*8(a-h,o-z)
       common/wn/wn
        if(x  <=  2.d0) then
           p = x+1.d0; return
        endif
        if(x > 2.d0 .and. x <= wn) then
           ix=floor(x+0.5d0); xx=dfloat(ix); p= sq(ix,x,xx); return
        endif
        if(x > wn) then
           p= r(x); return
        endif
        contains
        function sq(ix,x,xx)  ! left-right quadric pair
        common/protiarithmi/q(10000000)
        if(x <= xx) then
           sq=-2.d0*(q(ix)-q(ix-1)-1.d0)*(x-xx)**2+x-xx+q(ix)
        else
           sq=2.d0*(q(ix+1)-q(ix)-1.d0)*(x-xx-0.5d0)**2+            &
        (2.d0*(q(ix+1)-q(ix))-1.d0)*(x-xx-0.5d0)+(q(ix)+q(ix+1))/2.d0
        endif
        end function sq
        end function p
!
!
       function dp(x) ! 1 <= x < infinity ; derivative of p(x) 
       implicit real*8(a-h,o-z)
       common/wn/wn
        if(x  <=  2.d0) then
           dp = 1.d0; return
        endif
        if(x > 2.d0 .and. x <= wn) then
           ix=floor(x+0.5d0); xx=dfloat(ix); dp= dsq(ix,x,xx); return
        endif
        if(x > wn) then
           dp=dr(x); return
        endif
        contains
        function dsq(ix,x,xx) ! derivative of left-right quadrics
        common/protiarithmi/q(10000000)
        if(x <= xx) then
           dsq=-4.d0*(q(ix)-q(ix-1)-1.d0)*(x-xx)+1.d0
        else
           dsq=4.d0*(q(ix+1)-q(ix)-1.d0)*(x-xx-0.5d0)+               &
               2.d0*(q(ix+1)-q(ix))-1.d0
        endif
        end function dsq
        end function dp
!
!
       function p_(x)   ! 2 <= x < infinity 
       implicit real*8(a-h,o-z)
       common/protiarithmi/q(10000000)/wn/wn/ck/vkoch/nn/nn
1      if(x >= 3.d0 .and. x <= q(nn)) then
         if(x <= 1.d3) ixx=floor(x/dlog(x))
         if(x >  1.d3) ixx=floor(dlogintegral(x)-vkoch*dsqrt(x)*dlog(x))
         if(q(ixx) == 0.d0) qq=r(dfloat(ixx))
         if(q(ixx) /= 0.d0) qq=q(ixx)
         if(qq <= x) then
            is=-1; si=-1.d0
         else
            is= 1; si= 1.d0
         endif
       else
         if(x < 3.d0) then
            p_=x-1.d0; return
         endif
       endif
       do i=1, nn
          ii=ixx-is*i
          if(si*(q(ii)-x) <= 0.d0) then
          ix=ii; goto 3; endif
       enddo
3      if(dabs(q(ix)-x) >= dabs(q(ix+is)-x)) ix=ix+is
       xx=q(ix); xi=dfloat(ix); p_= y(ix,xi,x,xx)
       if(vkoch /= 0.d0) goto 1
       contains
       function y(ix,xi,x,xx)   ! left-right inverse quadric pair
       implicit real*8(a-h,o-z)
       common/protiarithmi/q(10000000)/ck/vkoch
              vkoch=0.d0
        if(x <= xx) then
           a=q(ix)-q(ix-1)-1.d0; b=8.d0*a*(q(ix)-x)+1.d0
           if(b <= 0.d0) then
              vkoch=1.d0
           else
              if(a == 0.d0) then
                 y=x-1.d0; else; y=xi+(1.d0-dsqrt(b))/(4.d0*a); endif
           endif
        else
           a=q(ix+1)-q(ix)-1.d0; b=8.d0*a*(x-q(ix))+1.d0
           if(b < 0.d0) then
              vkoch=1.d0
           else
              y=xi+(dsqrt(b)-1.d0)/(4.d0*a)
           endif
        endif
       end function y
       end function p_
!
!
       function dp_(x)    ! 2 <= x < infinity; inverse of dp(x) 
       implicit real*8(a-h,o-z)
       common/protiarithmi/q(10000000)/wn/wn/ck/vkoch/nn/nn
1      if(x >= 3.d0 .and. x <= q(nn)) then
         if(x <= 1.d3) ixx=floor(x/dlog(x))
         if(x >  1.d3) ixx=floor(dlogintegral(x)-vkoch*dsqrt(x)*dlog(x))
         if(q(ixx) == 0.d0) qq=r(dfloat(ixx))
         if(q(ixx) /= 0.d0) qq=q(ixx)
          if(qq <= x) then
             is=-1; si=-1.d0
          else
             is= 1; si= 1.d0
          endif
       else
         if(x < 3.d0) then
            dp_=1.d0
            return
         endif
       endif
       do i=1, nn
          ii=ixx-is*i
          if(si*(q(ii)-x) <= 0.d0) then
          ix=ii; goto 3; endif
       enddo
3      if(dabs(q(ix)-x) >= dabs(q(ix+is)-x)) ix=ix+is
       xx=q(ix); xi=dfloat(ix); dp_= dy(ix,x,xx)
       if(vkoch /= 0.d0) goto 1
       contains
       function dy(ix,x,xx)   ! left-right inverse quadric pair
       common/protiarithmi/q(10000000)/ck/vkoch
        vkoch=0.d0
        if(x <= xx) then
           a=q(ix)-q(ix-1)-1.d0; b=1.d0+8.d0*a*(q(ix)-x)
           if(b <= 0.d0) then
              vkoch=1.d0
           else
              dy=1.d0/dsqrt(b)
           endif
        else
           a=q(ix+1)-q(ix)-1.d0; b=1.d0+8.d0*a*(x-q(ix))
           if(b <= 0.d0) then
              vkoch=1.d0
           else
              dy=1.d0/dsqrt(b)
           endif
        endif
       end function dy
       end function dp_
!
!
         function r(t)!asymptote of the function p(t) (M. Cipolla, 1902)
         implicit real*8(a-h,o-z)
         common/rcorr/rc,drc
       r=t*(dlog(t)+dlog(dlog(t))+(dlog(dlog(t))-2.d0)/              &
         dlog(t)-((dlog(dlog(t)))**2-6.d0*dlog(dlog(t))+11.d0)/      &
         (2.d0*dlog(t)**2)-1.d0)+rc
       end function r
!
!
         function dr(x)  ! derivative of the asymptote r(t)
         implicit real*8(a-h,o-z)
         common/rcorr/rc,drc
      t1 = dlog(x); t2 = dlog(t1); t3 = t1**2; t4 = t3*t1
      t7 = t1-2.d0; t8 = dlog(t7); t18 = t2**2; t21 = t3**2
      t35 = -56.d0-4.d0*t2*t4+4.d0*t8*t1+50.d0*t1-13.d0*t3+          &
             2.d0*t4-6.d0*t8*t3-26.d0*t2*t1+4.d0*t18*t1+             &
             2.d0*t2*t21+2.d0*t8*t4-1.d0*t18*t3+6.d0*t2*t3+          &
             28.d0*t2-4.d0*t18-4.d0*t21+2.d0*t21*t1
      dr = 0.5d0*t35/t4/t7+drc
      end function dr
!
!
      function dlogintegral(x)
      implicit real*8(a-h,o-z)
      parameter(itol=100, small=1.d-10)
      dli=0.d0; i=0
1     i=i+1; dliw=dli; dm=(dlog(x))**i
      if(dm > small) then
         dli=dli+dfloat(factorial(i-1))/dm
         if(dabs(dliw-dli) > small.and.i <= itol) goto 1
      endif
      dlogintegral=x*dli
      contains
      integer recursive function factorial(l) result(lf)
      integer (4) l
      lf=1
      if(l > 0) then
         lf=l*factorial(l-1); return
      endif
      end function factorial
      end function dlogintegral
!
!
       subroutine eratosthenes(k)
       implicit real*8(a-h,o-z)
       common/protiarithmi/q(10000000)/rcorr/rc,drc/wn/wn/nn/nn
       common/ck/vkoch
       rc=0.d0; drc=0.d0; np=0; vkoch=0.d0
       do n=2,k
          id=1
1         id=id+1
          if(n == 2) np=1; q(1)=2.d0
          if(mod(n,id) /= 0) then
             if(id*id < n) goto 1
             np=np+1; q(np)=dfloat(n)
          endif
       enddo
       wn=dfloat(np); nn=np; rc=p(wn)-r(wn); drc=dp(wn)-dr(wn)
       end subroutine eratosthenes