Silverfrost Forums

Welcome to our forums

Riemann Zeta Function

20 Dec 2023 8:15 #30879

Hello,

I would like to convert the following Scilab code that I generated into Fortran 95, but not entirely sure of how to deal with it. Two issues, one how to deal with an equivalent to linspace and secondly the recursive function in zeta_3.

// Riemann Zeta Function (all domains)

// Re > 1 function s=zeta_1(z,n) n=linspace(1,n,n); //s=0.0; if z == 0 s = -0.5 elseif z == 1 s = %inf else s=sum(n.^-z); end

endfunction

// Re < 0 function zfn = zeta_2(s,n) // Riemann's functional equation // Analytic contiuation for negative values zfn = 2.s .* %pi.(s - 1) .* sin(%pi.s./2) . gamma(1 - s) .* zeta_1((1 - s),n) endfunction

// 0 < Re < 1 function zs1 = zeta_3(s,n) // Vectorised version
zs1=0 k=linspace(1,n,n); zs1 = sum((-1).^(k+ 1)./(k.s )); zs1 = 1./(1 - 2.( 1-s )).*zs1; endfunction

Any suggestions would be most helpful. Thanks

Lester

20 Dec 2023 10:31 #30880

There are no zeta functions here, if n is finite.

Try this for your zeta_1:

program tstz1
   print *, zta_1(0.5,10),zta_1(-0.4,20)
contains
   function zta_1(z,n) result(s)
   integer,allocatable :: iv(:)
   real s, z
   integer i,n
   allocate (iv(n))
   iv = (/(i,i=1,n)/)
   s = sum(z**iv)
   end function zta_1
end program

There is no recursion in your zeta_2. A function name behaves as an ordinary variable when no argument list follows it.

20 Dec 2023 1:31 #30881

Thanks for the solution, I have coded and tested the three functions and all looks fine:

program Riemann_Zeta_test

   print *, zeta_1(2.0,1000), zeta_2(-2.0,1000), zeta_2(-3.0,1000), zeta_3(0.5,1000), zeta_3(0.6,1000)
   
contains
   
! Re > 1.0 (if z=0 then s=-0.5 and if z=1 then s=infinity)
   function zeta_1(z,n) result(s)
   integer,allocatable :: nn(:)
   real s, z
   integer i,n
   allocate (nn(n))
   nn = (/(i,i=1,n)/)
   s = sum(nn**(-1*z))
   end function zeta_1

! Re < 0.0
   function zeta_2(s,n) result(zfn)
   integer,allocatable :: nn(:)
   integer i,n
   real s
   real :: pi= 4.0*ATAN(1.0)
   allocate (nn(n))
   nn = (/(i,i=1,n)/)
   zfn=2**s*pi**(s-1)*sin(pi*s/2)*gamma(1-s)*zeta_1((1-s),n)
   end function zeta_2

! 0.0 < Re < 1.0
   function zeta_3(s,n) result(zs1)
   integer,allocatable :: nn(:)
   integer i,n
   real s
   allocate (nn(n))
   nn = (/(i,i=1,n)/)   
   zs1 = sum((-1)**(nn + 1)/(nn**s ))
   zs1 = 1/(1 - 2**( 1-s ))*zs1  
   end function zeta_3

end program Riemann_Zeta_test

The values at 0.0 and 1.0 are already defined as -0.5 and infinity respectively.

Quoted from mecej4 There are no zeta functions here, if n is finite.

Try this for your zeta_1: [code]program tstz1 print *, zta_1(0.5,10),zta_1(-0.4,20) contains function zta_1(z,n) result(s) integer,allocatable :: iv(:) real s, z integer i,n allocate (iv(n)) iv = (/(i,i=1,n)/) s = sum(z**iv) end function zta_1 end program

There is no recursion in your zeta_2. A function name behaves as an ordinary variable when no argument list follows it.

20 Dec 2023 5:17 #30882

Given the function Zeta_3

   function zeta_3(s,n) result(zs1)
   integer,allocatable :: nn(:)
   integer i,n
   real s
   allocate (nn(n))
   nn = (/(i,i=1,n)/)   
   zs1 = sum((-1)**(nn + 1)/(nn**s ))
   zs1 = 1/(1 - 2**( 1-s ))*zs1  
   end function zeta_3

How can the variable (s) be passed as a complex number of the form s = 0.5 + it , where t would be a linear vector range and i is the imaginary number?

Thanks

21 Dec 2023 2:41 #30883

It is not clear to me what you want. If t is a real array, s = 0.5 + i.t is a complex array, not a complex number.

Will the following do, perhaps?

program tz3
  integer, parameter :: k = 7
  integer j
  complex s(k)
  do j = 1, k
     s(j) = cmplx(0.5,j*0.1)
  end do
  print *,zeta_3(s,5)
contains
elemental function zeta_3(s,n) result(zs1)
   integer,allocatable :: nn(:)
   integer, intent(in) :: n
   complex zs1
   complex, intent(in) :: s
   allocate (nn(n))
   nn = (/(i,i=1,n)/)
   zs1 = sum((-1)**(nn + 1)/(nn**s ))
   zs1 = 1/(1 - 2**( 1-s ))*zs1
end function zeta_3
end program
21 Dec 2023 12:42 #30884

That worked fine, many thanks for the guidance.

Had to stick the write statement in a do-loop to get the columns right:

do i = 1, k write(25,'(f15.12,1x,f15.12)') real_part(i), imaginary_part(i) ! write data to file end do

Quoted from mecej4 It is not clear to me what you want. If t is a real array, s = 0.5 + i.t is a complex array, not a complex number.

Will the following do, perhaps?

program tz3
  integer, parameter :: k = 7
  integer j
  complex s(k)
  do j = 1, k
     s(j) = cmplx(0.5,j*0.1)
  end do
  print *,zeta_3(s,5)
contains
elemental function zeta_3(s,n) result(zs1)
   integer,allocatable :: nn(:)
   integer, intent(in) :: n
   complex zs1
   complex, intent(in) :: s
   allocate (nn(n))
   nn = (/(i,i=1,n)/)
   zs1 = sum((-1)**(nn + 1)/(nn**s ))
   zs1 = 1/(1 - 2**( 1-s ))*zs1
end function zeta_3
end program
21 Dec 2023 3:42 #30885

You could replace the explicit DO loop

do i = 1, k
write(25,'(f15.12,1x,f15.12)') real_part(i), imaginary_part(i) ! write data to file
end do 

with the implied DO loop (available since Fortran 77)

write(25,'(f15.12,1x,f15.12)') (real_part(i), imaginary_part(i), i=1,k)
Please login to reply.