module simpsons_rule integer::r2=2**30 contains real (kind = 8) function f(x) implicit none real (kind=8),intent(in)::x f=sqrt(r2-x**2.) if(isnan(f)) then f=0.0 end if return end function f real (kind=8) function simple(a,b) implicit none real (kind=8),intent(in)::a,b simple=((b-a)/6.)*(f(a)+4*f((a+b)/2.)+f(b)) return end function simple real (kind=8) function composite(a,b,accuracy) implicit none real , intent(in)::a,b integer, intent(in)::accuracy integer::n,i real (kind=8)::h,integral integral=0.0 n=accuracy if(mod(n,2)==0) then n=n+1 end if h=(b-a)/n do i=0,n-1 integral=integral+simple(a+i*h,a+(i+1)*h) end do composite=integral return end function composite end module simpsons_rule program pi use simpsons_rule implicit none integer::max_t=18 integer::t,f_e_r_t,c=3 real::r r=floor(sqrt(real(r2))) write(*,*)"Approximations of pi:" do t=1,max_t,2 f_e_r_t=3**t write(*,*)"n=",f_e_r_t,":",composite((-1)*r,r,f_e_r_t)*(2./r2) end do write(*,*)"n=",2**30,":",composite((-1)*r,r,2**23)*(2./r2) end program pi