!PI con cifre decimaliVersione febbraio 2008 ( linguaggio fortran 95 compilabile con www.g95.org )
module decimale
implicit none
integer,public,parameter::p=selected_real_kind(12,300)
integer,public,dimension(:),allocatable::pi,vp,vpm,vpq,rs1,rs2
integer,public::base=100000
public:: fa_pi,stampa,imposta,valezero,cambiosegno,sommo,divido
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) then
write(unit=*,fmt=*)k/8,n
end if
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)) then
exit
end if
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(len=*),intent(in),optional::str
integer::j,nu
nu=6
if(present(nw)) then
nu=nw
end if
if(present(str)) then
write(unit=nu,fmt=*)str
end if
do j=1,size(vc,1),10
if(j+9>size(vc,1)) then
exit
end if
write(unit=nu,fmt=*) stri(vc(j:j+9))
end do
if(nu==6) then
read(unit=*,fmt=*)
end if
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) then
return
end if
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) then
return
end if
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(len=*),parameter::t=char(60),x=char(62)//char(60)
integer::ncf,t1,t2
write(unit=*,fmt=*)"Calcolo delle cifre decimali di PI"
write(unit=*,fmt=*)"(durata del calcolo proporzionale al quadrato delle cifre richieste)"
write(unit=*,fmt=*)" "
write(unit=*,fmt=*)"Assegna quante cifre (consigliabile almeno un centinaio...)"
write(unit=*,fmt=*)"...con oltre 2000 cifre salva su file pi-cifre.htm ..."
read(unit=*,fmt=*) ncf
call imposta(ncf)
write(unit=*,fmt=*)"Usa sei vettori di ",size(pi,1)," cifre ossia in tutto",24*size(pi,1),&
" ottetti"
write(unit=*,fmt=*)
if(ncf>2000) then
open(file="pi-cifre.htm",unit=12,status="replace",action="write")
write(unit=12,fmt="(1a)")t//"html"//x//"head"//x//"title>pi"
write(unit=12,fmt="(1a)")t//"/title"//x//"/head"//x//"body>"
write(unit=*,fmt=*)"Scrivera' sul file pi-cifre.htm"
write(unit=*,fmt=*)"Attendere..."
call system_clock(t1)
write(unit=*,fmt=*)"t1=",t1
call fa_pi()
call system_clock(t2)
write(unit=*,fmt=*)"t2=",t2," Diff.",t2-t1
call stampa(pi,t//"h2>π"//t//"/h2"//x//"pre>",12)
write(unit=12,fmt="(1a)")t//"/pre"//x//"/body"//x//"/html>"
close(unit=12)
write(unit=*,fmt=*)"Un invio per finire..."
read(unit=*,fmt=*)
else
write(unit=*,fmt=*)"Attendere...."
call fa_pi()
call stampa(pi,"pi")
end if
stop
end program fa
!
Dopo avere istallato il g95 ed aver chiamato questo file pi-dieci.f03 l'eseguibile si ottiene col comando: g95 -std=F pi-dieci.f03 -o pi-dieci.exe