module simpsons_rule contains real (kind = 8) function f(x) implicit none real (kind=8),intent(in)::x f=sqrt(16.-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::i write(*,*)"Approximations of pi:" do i=50000,10000000,1000000 write(*,*)"n=",i+1,composite(-4.0,4.0,i+1)/8. end do write(*,*)"n=",20000000,composite(-4.,4.,20000000)/8; end program pi