Commit e2de68e8 authored by Luca's avatar Luca

utils modul für eigene printfunktionen etc, funktionen in functions.f90

parent 74156e8d
MODULE functions
implicit none
CONTAINS
FUNCTION dilog(x,n)
! Computes the dilog Li_2(x) using the series representation up to order n
integer :: n
complex :: x, dilog
integer :: i
integer :: j(n)
j = (/(i, i=1,n,1)/)
dilog = sum(x**j / j**2)
END FUNCTION dilog
FUNCTION polylog(m,x,n)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m,n
complex :: x, polylog
integer :: i
integer :: j(n)
j = (/(i, i=1,n,1)/)
polylog = sum(x**j / j**m)
END FUNCTION polylog
recursive FUNCTION multiple_polylog(m, x, n) result(res)
! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
integer :: m(:)
complex :: x(:)
complex :: res
integer :: n, i
! print*, 'm = ', m
! print*, 'x = ', x
! print*, 'n = ', n
if(size(m) /= size(x)) then
print*, 'Error: m and x must have the same length'
end if
if(size(m) == 1) then
! base case
! print*, 'landed in base case'
res = polylog(m(1),x(1), n)
else
! recursion step
! print*, 'landed in step'
res = 0
do i = 1, n
res = res + x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i - 1)
end do
! a nicer way to do it would be but problem is i
! i = (/(j, j=1,n, 1)/)
! res = sum( x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i(1) - 1) )
end if
END FUNCTION multiple_polylog
END MODULE functions
PROGRAM test
! Used to test the procedures defined in this module
use functions
use utils
implicit none
complex :: result
! print*, dilog((0.8,0),100) ! should be 1.07479 + 0i
! print*, dilog((0.2,0.5),100) ! should be 0.133909 + 0.537628i
! print*, polylog(2,(0.2,0.5),100) ! should be 0.133909 + 0.537628i
! print*, polylog(5,(0.2,0.5),100) ! should be 0.192872345 + 0.505898833i
! result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /), 100)
! print*, 'result = ', result
! result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 100)
! print*, 'result = ', result
END PROGRAM test
#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;
int main()
{
symbol x("x"), y("y");
ex poly;
for (int i=0; i<3; ++i)
poly += factorial(i+16)*pow(x,i)*pow(y,2-i);
cout << poly << endl;
return 0;
}
! Contains some functions that might be useful later
! Write your own print function with ability to suppress print
! Muss immer alle prints und warnings ausschalten können
! Test Programm schreiben mit exit codes -> gfortran 'test.f90' und dann 'echo $?'
! Define GPL infinity
! Mach n optional
! Kommentar schreiben zu anderer Notation
MODULE utils
implicit none
logical :: print_enabled = .true.
logical :: warnings_enabled = .true.
CONTAINS
FUNCTION dilog(x,n)
! Computes the dilog Li_2(x) using the series representation up to order n
integer :: n
complex :: x, dilog
integer :: i
integer :: j(n)
j = (/(i, i=1,n,1)/)
dilog = sum(x**j / j**2)
END FUNCTION dilog
FUNCTION polylog(m,x,n)
! Computes the classical polylogarithm Li_m(x) using series representation up to order n
integer :: m,n
complex :: x, polylog
integer :: i
integer :: j(n)
j = (/(i, i=1,n,1)/)
polylog = sum(x**j / j**m)
END FUNCTION polylog
recursive FUNCTION multiple_polylog(m, x, n) result(res)
! Computes the multiple polylogarithm Li_{m1,...,mk} (x1,...,xk) up to order n
integer :: m(:)
complex :: x(:)
complex :: res
integer :: n, i
! print*, 'm = ', m
! print*, 'x = ', x
! print*, 'n = ', n
if(size(m) /= size(x)) then
print*, 'Error: m and x must have the same length'
subroutine print(s1,s2,s3,s4,s5)
character(len = *), intent(in), optional :: s1, s2, s3, s4, s5
if(print_enabled) then
print*, s1, s2, s3, s4, s5
end if
if(size(m) == 1) then
! base case
! print*, 'landed in base case'
res = polylog(m(1),x(1), n)
else
! recursion step
! print*, 'landed in step'
res = 0
do i = 1, n
res = res + x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i - 1)
end do
end subroutine print
! a nicer way to do it would be but problem is i
! i = (/(j, j=1,n, 1)/)
! res = sum( x(1)**i / i**m(1) * multiple_polylog(m(2:), x(2:), i(1) - 1) )
subroutine warn(s1,s2,s3,s4,s5)
character(len = *), intent(in), optional :: s1, s2, s3, s4, s5
if(warnings_enabled) then
print*, 'Warning: ', s1, s2, s3, s4, s5
end if
END FUNCTION multiple_polylog
end subroutine warn
END MODULE utils
PROGRAM test
! Used to test the procedures defined in this module
use utils
implicit none
complex :: result
! print*, dilog((0.8,0),100) ! should be 1.07479 + 0i
! print*, dilog((0.2,0.5),100) ! should be 0.133909 + 0.537628i
! print*, polylog(2,(0.2,0.5),100) ! should be 0.133909 + 0.537628i
! print*, polylog(5,(0.2,0.5),100) ! should be 0.192872345 + 0.505898833i
result = multiple_polylog((/ 5 /),(/ (0.2,0.5) /), 100)
print*, 'result = ', result
result = multiple_polylog((/ 5, 5 /),(/ (0.8,0),(0.3,0.5) /), 100)
print*, 'result = ', result
END PROGRAM test
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment