!Microlibreria di algebra modulareMicrolibreria di algebra modulare programmata in Fortran 95
Attenzione: questo è un sorgente Fortran compilabile senza modifiche con il gratuito Fortran Silverfrost o con l'altrettanto gratuito G95 ma anche non il non gratuito ma molto professionale compilatore INTEL Visual Fortran
Vedere: http://www.elegio.it/doc/tn/ ed in particolare la pagina http://www.elegio.it/doc/tn/altaprecisione_modulare.html
module algebramodulare
implicit none
!
integer,parameter::cifre=4 ! == elementi del vettore usato
integer,parameter::cifrx=9 ! == potenza del 10 della base (il max e' 15)
!
integer,parameter::q=selected_real_kind(15,300)
real(kind=q),parameter::basecifre=10.0_q**min(cifrx,15)
type,public::zahlen
!private
real(kind=q),dimension(0:cifre-1)::z
end type zahlen
type(zahlen):: zahlenmod
!
interface assignment(=)
module procedure zahlcost,zahlvet,zahlivet,stringazahlen,numerone
end interface
interface operator(+)
module procedure fisomcomod,fi2somcomod,frsomcomod,fr2somcomod
module procedure fsommamod
end interface
interface operator(-)
module procedure fidifcomod,fi2difcomod,frdifcomod,fr2difcomod
module procedure fdiffemod
end interface
interface operator(*)
module procedure prodotto
end interface
interface operator(**)
module procedure potenza
end interface
!
private:: zahlcost,zahlvet,zahlivet,stringazahlen,numerone, &
fisomcomod,fi2somcomod,frsomcomod,fr2somcomod, &
fsommamod, &
fidifcomod,fi2difcomod,frdifcomod,fr2difcomod, &
fdiffemod,prodotto,potenza
!
public:: inizahlen,somcomod,sommamod,difcomod,diffemod,nonzero,nonuno, &
menuno,dimezzo,numeropari,trascrivo
!
contains
subroutine inizahlen()
zahlenmod%z=0.0_q
zahlenmod%z(cifre-1)=1.0_q
return
end subroutine inizahlen
!
subroutine zahlcost(z,j)
integer,intent(in)::j
type(zahlen),intent(out)::z
z%z(0:)=0.0_q
z%z(0)=j
z%z(0)=modulo(z%z(0),basecifre)
return
end subroutine zahlcost
!
subroutine zahlvet(z,v)
real(kind=q),dimension(0:),intent(in)::v
type(zahlen),intent(out)::z
integer::j
z%z(0:)=0.0_q
do j=0,size(v,1)-1
z%z(j)=modulo(v(j),basecifre)
end do
return
end subroutine zahlvet
!
subroutine zahlivet(z,iv)
integer,dimension(0:),intent(in)::iv
type(zahlen),intent(out)::z
integer::j
z%z(0:)=0.0_q
do j=0,size(iv,1)-1
z%z(j)=iv(j)
z%z(j)=modulo(z%z(j),basecifre)
end do
return
end subroutine zahlivet
!
subroutine numerone(z,str)
character(len=*),intent(in)::str
type(zahlen),intent(out)::z
real(kind=q)::qe
integer::j,m,n,kk
character(len=len(str))::s
do j=0,cifre-1
z%z(j)=0
end do
n=0
s=" "
do j=1,len(str)
m=ichar(str(j:j))
if(m>32) then
n=n+1
s(n:n)=char(m)
end if
end do
if(n==0) return
qe=1.0_q
kk=0
do j=n,1,-1
z%z(kk)=z%z(kk)+modulo(ichar(s(j:j))+2,10)*qe
qe=qe*10.0_q
if(qe>=basecifre) then
qe=1.0_q
kk=min(kk+1,cifre-1)
end if
end do
return
end subroutine numerone
!
subroutine stringazahlen(str,z)
type(zahlen),intent(in)::z
character(len=*),intent(out)::str
character(len=cifre*cifrx)::testo
real(kind=q)::zz
integer::j,k,n,p,np
p=0
do j=0,cifre-1
zz=z%z(j)
do k=1,cifrx
p=p+1
n=modulo(zz,10.0_q)
testo(p:p)=char(48+n)
zz=(zz-n)/10.0_q
if(n>0) np=p
end do
end do
k=min(len(str),np)
str=" "
do j=1,k
str(j:j)=testo(k+1-j:k+1-j)
end do
if(1>len_trim(str)) str="0"
return
end subroutine stringazahlen
!
function fisomcomod(z,c)result(s)
type(zahlen),intent(in)::z
integer,intent(in)::c
type(zahlen)::s
real(kind=q)::zz
zz=c
call somcomod(s%z,z%z,zz)
return
end function fisomcomod
!
function fi2somcomod(c,z)result(s)
type(zahlen),intent(in)::z
integer,intent(in)::c
type(zahlen)::s
real(kind=q)::zz
zz=c
call somcomod(s%z,z%z,zz)
return
end function fi2somcomod
!
function frsomcomod(z,zz)result(s)
type(zahlen),intent(in)::z
real(kind=q),intent(in)::zz
type(zahlen)::s
call somcomod(s%z,z%z,zz)
return
end function frsomcomod
!
function fr2somcomod(zz,z)result(s)
type(zahlen),intent(in)::z
real(kind=q),intent(in)::zz
type(zahlen)::s
call somcomod(s%z,z%z,zz)
return
end function fr2somcomod
!
subroutine somcomod(sz,vz,a)
real(kind=q),dimension(0:),intent(in)::vz
real(kind=q),intent(in)::a
real(kind=q),dimension(0:),intent(out)::sz
real(kind=q)::rcs
integer::j
rcs=min(abs(a),basecifre-1.0_q)
do j=0,cifre-1
sz(j)=rcs+vz(j)-zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else if(0>sz(j)) then
sz(j)=sz(j)+basecifre
rcs=-1.0_q
else
rcs=0.0_q
end if
end do
if(0>rcs) then
rcs=0.0_q
do j=0,cifre-1
sz(j)=sz(j)+rcs+zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else
rcs=0.0_q
end if
end do
end if
return
end subroutine somcomod
!
function fsommamod(z,zz)result(s)
type(zahlen),intent(in)::z,zz
type(zahlen)::s
call sommamod(s%z,z%z,zz%z)
return
end function fsommamod
!
subroutine sommamod(sz,vz,wz)
real(kind=q),dimension(0:),intent(in)::vz,wz
real(kind=q),dimension(0:),intent(out)::sz
real(kind=q)::rcs
integer::j
rcs=0.0_q
do j=0,cifre-1
sz(j)=rcs+vz(j)+wz(j)-zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else if(0>sz(j)) then
sz(j)=sz(j)+basecifre
rcs=-1.0_q
else
rcs=0.0_q
end if
end do
if(0>rcs) then
rcs=0.0_q
do j=0,cifre-1
sz(j)=sz(j)+rcs+zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else
rcs=0.0_q
end if
end do
end if
return
end subroutine sommamod
!
function fidifcomod(z,c)result(s)
type(zahlen),intent(in)::z
integer,intent(in)::c
type(zahlen)::s
real(kind=q)::zz
zz=c
call difcomod(s%z,z%z,zz)
return
end function fidifcomod
!
function fi2difcomod(c,z)result(s)
type(zahlen),intent(in)::z
integer,intent(in)::c
type(zahlen)::ss,s
real(kind=q)::zz
call diffemod(ss%z,zahlenmod%z,z%z)
zz=c
call somcomod(s%z,ss%z,zz)
return
end function fi2difcomod
!
function frdifcomod(z,zz)result(s)
type(zahlen),intent(in)::z
real(kind=q),intent(in)::zz
type(zahlen)::s
call difcomod(s%z,z%z,zz)
return
end function frdifcomod
!
function fr2difcomod(zz,z)result(s)
type(zahlen),intent(in)::z
real(kind=q),intent(in)::zz
type(zahlen)::s,ss
call diffemod(ss%z,zahlenmod%z,z%z)
call somcomod(s%z,ss%z,zz)
return
end function fr2difcomod
!
subroutine difcomod(sz,vz,a)
real(kind=q),dimension(0:),intent(in)::vz
real(kind=q),intent(in)::a
real(kind=q),dimension(0:),intent(out)::sz
real(kind=q)::rcs
integer::j
rcs=min(-abs(a),basecifre-1.0_q)
do j=0,cifre-1
sz(j)=rcs+vz(j)
if(0>sz(j)) then
sz(j)=sz(j)+basecifre
rcs=-1.0_q
else
rcs=0.0_q
end if
end do
if(0>rcs) then
rcs=0.0_q
do j=0,cifre-1
sz(j)=sz(j)+rcs+zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else
rcs=0.0_q
end if
end do
end if
return
end subroutine difcomod
!
function fdiffemod(z,zz)result(s)
type(zahlen),intent(in)::z,zz
type(zahlen)::s
call diffemod(s%z,z%z,zz%z)
return
end function fdiffemod
!
subroutine diffemod(sz,vz,wz)
real(kind=q),dimension(0:),intent(in)::vz,wz
real(kind=q),dimension(0:),intent(out)::sz
real(kind=q)::rcs
integer::j
rcs=0.0_q
do j=0,cifre-1
sz(j)=rcs+vz(j)-wz(j)
if(0>sz(j)) then
sz(j)=sz(j)+basecifre
rcs=-1.0_q
else
rcs=0.0_q
end if
end do
if(0>rcs) then
rcs=0.0_q
do j=0,cifre-1
sz(j)=sz(j)+rcs+zahlenmod%z(j)
if(sz(j)>=basecifre) then
sz(j)=sz(j)-basecifre
rcs=1.0_q
else
rcs=0.0_q
end if
end do
end if
return
end subroutine diffemod
!
function nonzero(z)result(f)
type(zahlen),intent(in)::z
logical::f
integer::j
do j=cifre-1,0,-1
if(abs(z%z(j))>0.0_q) then
f=.true.
return
end if
end do
f=.false.
return
end function nonzero
!
function nonuno(z)result(f)
type(zahlen),intent(in)::z
logical::f
integer::j
do j=cifre-1,1,-1
if(z%z(j)>0.0_q) then
f=.true.
end if
end do
f=(z%z(0)-1.0_q)>0.0_q
return
end function nonuno
!
function menuno(z)result(f)
type(zahlen),intent(in)::z
logical :: f
integer::j
f=.true.
do j=cifre-1,1,-1
if(abs(z%z(j)-zahlenmod%z(j))>0.0_q) then
f=.false.
return
end if
end do
if(abs(z%z(0)-zahlenmod%z(0)+1.0_q)>0.0_q ) then
f=.false.
end if
return
end function menuno
!
subroutine dimezzo(pari,z)
type(zahlen),intent(in out)::z
logical,intent(out)::pari
real(kind=q)::res
integer::j
res=0.0_q
do j=cifre-1,0,-1
z%z(j)=z%z(j)+res
if(modulo(z%z(j),2.0_q)>0.0_q) then
z%z(j)=(z%z(j)-1.0_q)/2.0_q
res=basecifre
else
z%z(j)=z%z(j)/2.0_q
res=0.0_q
end if
end do
pari=0.0_q>=res
return
end subroutine dimezzo
!
function numeropari(z)result(f)
type(zahlen),intent(in)::z
logical::f
f=0.0_q>=modulo(z%z(0),2.0_q)
return
end function numeropari
!
subroutine trascrivo(z,a)
type(zahlen),intent(in)::a
type(zahlen),intent(out)::z
integer::j
do j=cifre-1,0,-1
z%z(j)=modulo(a%z(j),basecifre)
end do
return
end subroutine trascrivo
!
function prodotto(aa,bb)result(vs)
type(zahlen),intent(in):: aa,bb
type(zahlen)::vs
type(zahlen)::as,bs
logical::pari
call trascrivo(as,aa)
call trascrivo(bs,bb)
vs=0
do
call dimezzo(pari,as)
if(.not.pari) then
vs=vs+bs
end if
bs=bs+bs
if(.not.nonzero(as)) exit
end do
end function prodotto
!
function potenza(aa,n)result(vp)
type(zahlen),intent(in):: aa,n
type(zahlen)::vp
type(zahlen)::bp,mp
logical::pari
call trascrivo(bp,aa)
call trascrivo(mp,n)
vp=1
do
call dimezzo(pari,mp)
if(.not.pari) then
vp=bp*vp
end if
bp=bp*bp
if(.not.nonzero(mp)) exit
end do
end function potenza
!
function pseudoprimoforte(a,n)result(f)
type(zahlen),intent(in)::a,n
logical::f,pari
type(zahlen)::sb,vp,nn
integer::j
call trascrivo(nn,zahlenmod)
call trascrivo(zahlenmod,n)
call trascrivo(sb,n)
sb=n-1
j=0
do
j=j+1
vp=sb
call dimezzo(pari,sb)
if(.not.pari) then
sb=vp
exit
end if
end do
vp=a**sb
f=.not.nonuno(vp)
if(f) then
call trascrivo(zahlenmod,nn)
return
end if
do
f=menuno(vp)
if(f) exit
sb=sb+sb
vp=vp*vp
j=j-1
if(0>j) exit
end do
call trascrivo(zahlenmod,nn)
return
end function pseudoprimoforte
!
subroutine controllo_pseudoprimoforte()
type(zahlen)::a,mm
character(79):: ilnumero
integer::k
write(*,*)"Assegna il possibile numero primo"
read(*,"(1a)") ilnumero
if(1>len_trim(ilnumero)) return
mm=ilnumero
ilnumero=mm
write(*,*)"Analizza il numero:"
write(*,*) ilnumero
ilnumero=" "
k=2
do
a=k
if(pseudoprimoforte(a,mm) ) then
write(*,*)"Con la base",k," il numero sembra primo"
else
write(*,*)"Con la base",k," numero NON PRIMO"
end if
k=2*((k+1)/2)+1
write(*,*)"Digita 0 per finire"
read(*,"(1a)") ilnumero
if(ilnumero(1:1)=="0") exit
end do
return
end subroutine controllo_pseudoprimoforte
!
subroutine cercagermain()
type(zahlen)::a,uno,due,tre,cinque,mm,poi
character(79):: ilnumero
integer::k,nl,tot
logical::siono
character(len=*),parameter::fs="almod-germain.txt"
zahlenmod%z(:)=0.0_q
zahlenmod%z(cifre-1)=basecifre/10.0_q
!
tot=0
uno=1
due=2
tre=3
cinque=5
!
write(*,*)"Catene di primi di Sophie Germain"
write(*,*)" "
write(*,*)"Numeri notevoli trovati con questo programma"
write(*,*)"Livello 6: 89,63419,127139,405269,810809,...,181774409"
write(*,*)"Livello 7: 1122659,...,183561929"
write(*,*)"Livello 8: 19099919"
write(*,*)"Livello 9: 85864769"
write(*,*)" "
write(*,*)"Assegna il numero dispari di partenza:"
read(*,"(1a)") ilnumero
mm=ilnumero
ilnumero=mm
write(*,*)"Inizia dal numero:"
write(*,*) ilnumero
write(*,*)"Quanti livelli di Sophie Germain:"
read(*,*) nl
if(1>nl) return
nl=max(1,nl)
open(file=fs,unit=12,status="replace",action="write")
write(12,*)"Ricerca fino al livello",nl
close(unit=12)
write(*,*)"Ora indaga e si arresta quando lo trova: invia per iniziare"
read(*,*)
do
poi=mm
controllo:do k=1,nl
! ilnumero=poi
! write(*,*)"Controllo passo:",k
! write(*,*) ilnumero
siono=pseudoprimoforte(due,poi)
if(.not.siono) exit controllo
siono=pseudoprimoforte(tre,poi)
if(.not.siono) exit controllo
siono=pseudoprimoforte(cinque,poi)
if(.not.siono) exit controllo
poi=poi*due
poi=poi+uno
if(k>max(2,nl-6)) then
ilnumero=mm
open(file=fs,unit=12,status="replace",action="write",&
position="append")
write(12,*)"Livello ",k
write(12,*) ilnumero
close(12)
end if
end do controllo
if(siono) then
ilnumero=mm
write(*,*)"Pseudoprimo di SF di livello ",nl
write(*,*) ilnumero
open(file=fs,unit=12,status="replace",action="write",&
position="append")
write(12,*)"Pseudoprimo di SF di livello ",nl
write(12,*) ilnumero
close(12)
write(*,*)"Digita 1 per proseguire, altrimenti smetto..."
read(*,*) k
if(k/=1) exit
end if
mm=mm+2
tot=tot+1
if(modulo(tot,100000)==0) then
ilnumero=mm
write(*,*) ilnumero
tot=1
end if
end do
return
end subroutine cercagermain
!
end module algebramodulare
!
program almod
use algebramodulare
implicit none
character(len=78)::stringona
type(zahlen)::a,b,c
call inizahlen()
!
zahlenmod="100000"
stringona=zahlenmod
do
call controllo_pseudoprimoforte()
write(*,*)"....... digita 0 per finire ....."
read(*,"(1a)") stringona
if(stringona(1:1)=="0") exit
end do
call cercagermain()
write(*,*)"Amen"
read(*,*)
stop
end program almod
!