      program calcpi
      implicit real*8 (a-h,o-z)
      parameter(nmax=50001)
      i=0
      rad2=dble(nmax+1)**2
      do 10 n=1,nmax,2
        ann=dble(n)**2
      do 20 m=1,nmax,2
        amm=dble(m)**2
        rr=ann+amm
        if(rr.lt.rad2) i=i+1
 20   continue
 10   continue
      s=dble(i)/dble((nmax+1)/2)**2
      pi=4.0d0*s
      write(6,900) nmax,pi
 900  format(' nmax=',i6,' pi=',f20.16)
      end
