!pi plouffeπSorgente in Fortran 95 compilabile con www.g95.orgpi ( versione 20080214 )
module pi_pezzi
!
! Applica il metodo di estrazione delle cifre, recentemente scoperto
! per calcolare il gruppo di cifre di pi_greco che si desidera.
! Si possono trovare 2**24 ossia circa 16 milioni di cifre di pi_greco,
! un pezzettino alla volta e senza dover calcolare le precedenti....
! Notare che trova l'espressione di pi_greco in cifre esadecimali
! e quindi in cifre binarie. Non e' attualmente nota una formula
! che calcoli le cifre decimali di pi_greco e quindi, per trovare
! la, mettiamo, milionesima cifra decimale di pi_greco occorre tuttora
! conoscere tutte le cifre prima. In esadecimale invece no!
! La milionesima cifra dopo la virgola e' 2 ed e' sia preceduta
! che seguita dalla cifra 6.
!
implicit none
integer,public,parameter:: ntp=25
integer,public,parameter::pr=selected_real_kind(12,300)
real(kind=pr),public,dimension(ntp) :: tp = 0.0_pr
real(kind=pr),public,parameter:: eps = 1.0e-17_pr
character(len=*),public,parameter,dimension(0:15):: hx = (/ "0", "1", &
"2", "3", "4", "5", "6", "7", "8", "9", "A", "B", "C", &
"D", "E", "F" /)
private:: xexpm
public:: xseries,ihex
!
contains
!
subroutine ihex (x, nhx, chx)
!
! Il risultato e' la stringa chx che contiene nhx cifre
! esadecimali della parte frazionaria di x.
!
real(kind=pr),intent(in) :: x
real(kind=pr) :: y
integer,intent(in) :: nhx
integer :: i
character(len=*),dimension(1:),intent(out):: chx
!
y = abs (x)
do i = 1, nhx
y = 16.0_pr * (y - aint (y))
chx(i) = hx(int(y))
end do
return
end subroutine ihex
!
subroutine xseries ( series, m, ic)
!
! Questa funzione calcola la serie : sommatoria, indice k,
! di 16**(ic-k)/(8*k+m) usando la tecnica di esponentazione
! modulare.
!
real(kind=pr) :: ak, p, s, t
real(kind=pr),intent(out)::series
integer,intent(in) :: ic,m
integer :: k
!
s = 0.0_pr
!
! Somma le serie fino a ic.
!
do k = 0, ic - 1
ak = 8 * k + m
p = ic - k
call xexpm (t, p, ak)
s = modulo (s + t / ak, 1.0_pr)
end do
!
! Calcola un po' di termini dove k.ge.ic.
!
do k = ic, ic + 100
ak = 8 * k + m
t = 16.0_pr ** (ic - k) / ak
if (eps > t ) then
exit
end if
s = modulo (s + t, 1.0_pr)
end do
!
series = s
return
end subroutine xseries
!
subroutine xexpm(expm, p, ak)
!
! expm = mod( 16**p, ak). Questa funzione usa lo schema di
! esponentazione binario da sinistra a destra. E' valida per
! ogni ak non superiore a 2**24.
!
real(kind=pr),intent(in) :: p,ak
integer :: i
real(kind=pr),intent(out):: expm
real(kind=pr) :: p1, pt, r
integer::jexpm,jak !!!!
!
! Se questa e' la prima chiamata di expm, riempie la tabella
! di potenze di due tp.
!
if(0.0_pr >= tp(1)) then
tp(1) = 1.0_pr
do i = 2, ntp
tp(i) = 2.0_pr * tp(i-1)
end do
end if
!
if (0.0_pr >= abs(ak-1.0_pr) ) then
expm = 0.0_pr
return
end if
!
! Trova la piu' grande potenza di due minore o uguale a p.
!
do i = 1, ntp
if (tp(i) > p) then
exit
end if
end do
!
pt = tp(i-1)
p1 = p
r = 1.0_pr
!
! Attua l'algoritmo di esponenziazione binario modulo ak.
!
do
if (p1 >= pt) then
r = modulo(16.0_pr * r, ak)
p1 = p1 - pt
end if
pt = 0.5_pr * pt
if( 1.0_pr>pt) then
exit
end if
r = modulo(r * r, ak)
end do
!
expm = modulo (r, ak)
!
! Controllo che il risultato sia un intero:
!
jexpm=expm
jak=ak
if(abs(ak-jak)>0.0_pr) then
write(unit=*,fmt=*)"Residuo frazionario ak: ",ak,jak
read(unit=*,fmt=*)
end if
if(abs(expm-jexpm)>0.0_pr) then
write(unit=*,fmt=*)"Residuo frazionario expm: ",expm,jexpm
read(unit=*,fmt=*)
end if
return
end subroutine xexpm
!
end module pi_pezzi
!
! ------------------------------------------------------------------
!
program piqp
!
use pi_pezzi
implicit none
real(kind=pr) :: pid, s1, s2, s3, s4
integer :: ic,k
integer,parameter :: nhx = 12
character(len=1),dimension(nhx):: chx
!
! ic rappresenta la posizione della ultima cifra esadecimale considerata
! nota. Il risultato inizia alla cifra esadecimale ic+1_esima.
! Se si assegna ad es 10 si intende che si vuole dalla undicesima
! cifra inclusa in poi, dopo la virgola ossia si deve ottenere A308D...
!
write(unit=*,fmt=*)" Algoritmo di D.Bailey, P.Borwein e Simon Plouffe "
write(unit=*,fmt=*)"Trova le cifre esadecimali di Pi_greco senza dover"
write(unit=*,fmt=*)"conoscere quali sono quelle prima del punto di inizio."
write(unit=*,fmt=*)
do
assegnazione: do
write(unit=*,fmt=*) "Assegna da quale cifra ic (-1 per finire) "
read(unit=*,fmt=*) ic
if(2**24>ic) then
exit assegnazione
end if
end do assegnazione
if(0>ic) then
exit
end if
if(ic>100) then
write(unit=*,fmt=*) &
"Attendi, se hai chiesto un ic grande..."
end if
!
call xseries (s1, 1, ic)
call xseries (s2, 4, ic)
call xseries (s3, 5, ic)
call xseries (s4, 6, ic)
pid = modulo(4.0_pr * s1 - 2.0_pr * s2 - s3 - s4, 1.0_pr) + 1.0_pr
call ihex (pid, nhx, chx)
!
write(unit=*,fmt=*)" "
write(unit=*,fmt=*)"--Calcolo delle cifre esadecimali di Pi_greco--"
write(unit=*,fmt=*)" Le prime cifre esadecimali dovrebbero essere queste..."
write(unit=*,fmt=*)" 3.243F6 A8885 A308D 31319 8A2E0 37073 44A40 93822"
write(unit=*,fmt=*)"Posizione a partire da ",ic," + 1"
! write(unit=*,fmt=*)" ",pid
write(unit=*,fmt="(a)",advance="no")" "
do k=1,nhx
if(modulo(ic+k,5)==1) then
write(unit=*,fmt="(1a)",advance="no") " "
end if
write(unit=*,fmt="(1a)",advance="no") chx(k)
end do
write(unit=*,fmt=*)" "
end do
write(unit=*,fmt=*)"Amen..."
stop
end program piqp !
La compilazione si ottiene, avendo istallato il g95, e avendo dato a questo file il nome pi-plouffe.f03 , col comando:
g95 -std=F pi-plouffe.f03 -o pi-plouffe.exe