!PI con cifre decimali
module decimale
implicit none
integer,parameter::p=selected_real_kind(12,300)
integer,dimension(:),allocatable::pi,vp,vpm,vpq,rs1,rs2
integer::base=100000
contains !
!
! Calcola π con la formula di David Bailey, Peter Borwein e Simon Plouffe usata
! tradizionalmente ossia senza sfruttare la possibilità di calcolo veloce in caso
! di uso della base esadecimale.
!
subroutine fa_pi()
integer::k,kc,n
n=size(pi,1)
pi(1:)=0
vp(1:)=0
vp(1)=4*base/10
kc=max(64,8*min(n/23,100))
k=0
do
if(modulo(k,kc)+8==kc)write(*,*)k/8,n
call divido(rs1,vp,k+1)
call sommo (rs2,rs1,pi)
call divido(vpm,vp,2)
call divido(rs1,vpm,k+4)
call divido(vpq,vpm,2)
call divido(vp,vpq,k+5)
call sommo (vpm,rs1,vp)
call divido(rs1,vpq,k+6)
call sommo (vp,rs1,vpm)
call cambiosegno(rs1,vp)
call sommo(pi,rs2,rs1)
call divido(vp,vpq,4)
if(valezero(vp))exit
k=k+8
end do
return
end subroutine fa_pi
!
subroutine imposta(ncf)
integer,intent(in)::ncf
integer::n
n=50*(abs(ncf)/50)+75
allocate(pi(n),vp(n),vpm(n),vpq(n),rs1(n),rs2(n))
pi(:)=0
vp=pi
vpm=pi
vpq=pi
rs1=pi
rs2=pi
return
end subroutine imposta
!
subroutine stampa(vc,str,nw)
integer,dimension(:),intent(in)::vc
integer,intent(in),optional::nw
character(*),intent(in),optional::str
integer::j,nu
nu=6
if(present(nw))nu=nw
if(present(str))write(nu,*)str
do j=1,size(vc,1),10
if(j+9>size(vc,1)) exit
write(nu,*) stri(vc(j:j+9))
end do
if(nu==6) read(*,*)
return
contains
function stri(nn)result(s)
integer,intent(in),dimension(:)::nn
character(len=5*size(nn,1))::s
integer::k,kk,m
do k=1,size(nn,1)
m=nn(k)
kk=m/10000
s(5*k-4:5*k-4)=char(48+kk)
m=m-10000*kk
kk=m/1000
s(5*k-3:5*k-3)=char(48+kk)
m=m-1000*kk
kk=m/100
s(5*k-2:5*k-2)=char(48+kk)
m=m-100*kk
kk=m/10
s(5*k-1:5*k-1)=char(48+kk)
m=m-10*kk
s(5*k:5*k)=char(48+m)
end do
return
end function stri
end subroutine stampa
!
function valezero(va)result(z)
integer,dimension(:),intent(in)::va
logical::z
integer::k
do k=1,size(va,1)
z=va(k)==0
if(.not.z)return
end do
return
end function valezero
!
subroutine cambiosegno(wc,va)
integer,dimension(:),intent(in)::va
integer,dimension(:),intent(out)::wc
integer::j,jj,ns
ns=min(size(va,1),size(wc,1))
jj=0
do j=ns,1,-1
if(va(j)==0) then
wc(j)=0
else
wc(j)=base-va(j)
jj=j-1
exit
end if
end do
do j=jj,1,-1
wc(j)=base-va(j)-1
end do
return
end subroutine cambiosegno
!
subroutine sommo(wc,va,vb)
integer,dimension(:),intent(in)::va,vb
integer,dimension(:),intent(out)::wc
integer::j,ns,unm
ns=min(size(va,1),size(vb,1),size(wc,1))
unm=0
do j=ns,1,-1
unm=unm+va(j)+vb(j)
wc(j)=modulo(unm,base)
unm=unm/base
end do
return
end subroutine sommo
!
subroutine divido(wc,vc,nm)
integer,intent(in)::nm
integer,dimension(:),intent(in)::vc
integer,dimension(:),intent(out)::wc
integer:: ns,j,jj
real(kind=p)::unm,nmr
ns=min(size(vc,1),size(wc,1))
jj=0
do j=1,ns
if(vc(j)==0)then
wc(j)=vc(j)
else
jj=j
exit
end if
end do
if(jj==0) return
unm=0.0_p
nmr=nm
do j=jj,ns
unm=base*unm+vc(j)
wc(j)=unm/nm
unm=unm-wc(j)*nmr
end do
return
end subroutine divido
!
end module decimale
!
program fa
use decimale
implicit none
character(*),parameter::t=char(60),x=char(62)//char(60)
integer::ncf,t1,t2
write(*,*)"Calcolo delle cifre decimali di PI"
write(*,*)"(durata del calcolo proporzionale al quadrato delle cifre richieste)"
write(*,*)" "
write(*,*)"Assegna quante cifre (consigliabile almeno un centinaio...)"
read(*,*) ncf
call imposta(ncf)
write(*,*)"Usa sei vettori di ",size(pi,1)," cifre ossia in tutto",24*size(pi,1),&
" ottetti"
write(*,*)
if(ncf>2000) then
open(file="pi-cifre.htm",unit=12,status="replace",action="write")
write(12,"(1a)")t//"html"//x//"head"//x//"title>pi"
write(12,"(1a)")t//"/title"//x//"/head"//x//"body>"
write(*,*)"Scrivera' sul file pi-cifre.htm"
write(*,*)"Attendere..."
call system_clock(t1)
write(*,*)"t1=",t1
call fa_pi()
call system_clock(t2)
write(*,*)"t2=",t2," Diff.",t2-t1
call stampa(pi,t//"h2>π"//t//"/h2"//x//"pre>",12)
write(12,"(1a)")t//"/pre"//x//"/body"//x//"/html>"
close(12)
write(*,*)"Un invio per finire..."
read(*,*)
else
write(*,*)"Attendere...."
call fa_pi()
call stampa(pi,"pi")
end if
stop
end program fa
!