Commit 24959b26 authored by ulrich_y's avatar ulrich_y

Merge branch 'ieps'

parents 3ec2beb9 3215f5fb
MODE=DEBUG
HAVE_GINAC=1
FC=gfortran
AR=ar rcs
CC=gcc
MCC=mcc
ifeq ($(HAVE_GINAC),1)
LD=g++
LFLAGS=-lgfortran
else
LD=gfortran
endif
FFLAGS=-fdefault-real-8 -cpp -pedantic-errors -std=f2008
FFLAGS+= -Werror -Wall -Wno-maybe-uninitialized -Wno-uninitialized
......@@ -17,6 +23,11 @@ else
FFLAGS += -ffpe-trap=invalid,overflow -fdump-core -fbacktrace
endif
ifeq ($(HAVE_GINAC),1)
FFLAGS += -DHAVE_GINAC
endif
files=globals.o ieps.o utils.o shuffle.o maths_functions.o mpl_module.o gpl_module.o
objects = $(addprefix build/,$(files))
......@@ -31,6 +42,9 @@ build/%.o: src/%.f90
@echo "F90 $@"
@$(FC) $(FFLAGS) -c $< -J build -o $@
build/%.o: src/%.cpp
@echo "C++ $@"
@$(CC) -c $< -o $@
# Mathlink related
......@@ -61,11 +75,18 @@ gpl: build/gpl.o libgpl.a build/mcc.internals
eval: libgpl.a build/eval.o
@echo "LD $@"
@$(LD) -o $@ build/eval.o -L. -lgpl
@$(LD) -o $@ build/eval.o -L. -lgpl $(LFLAGS)
ifeq ($(HAVE_GINAC),1)
test: $(objects) build/ginac.o build/test.o
@echo "LD $@"
@$(LD) -o $@ $^ -lcln -lginac $(LFLAGS)
else
test: $(objects) build/test.o
@echo "LD $@"
@$(LD) -o $@ $^ $(LFLAGS)
$(LD) -o $@ $^ $(LFLAGS)
endif
check: test
./$<
......
#include <ginac/ginac.h>
using namespace GiNaC;
#include <cln/cln.h>
#include <stdio.h>
#include <iostream>
typedef struct {double r,i;} complex_t;
extern "C"{
complex_t geval_(complex_t * z, int* n);
};
complex_t geval_(complex_t * z, int* n) {
cln::cl_inhibit_floating_point_underflow = true;
lst w;
for(long i=0;i<(*n)-1;i++)
{
w.append((z->r)+(z->i)*I);
z++;
}
ex ans = G(w,z->r).evalf();
return {
.r = ex_to<numeric>(evalf(real_part(ans))).to_double(),
.i = ex_to<numeric>(evalf(imag_part(ans))).to_double()
};
}
......@@ -4,8 +4,6 @@ MODULE globals
integer, parameter :: prec = selected_real_kind(15,32)
real, parameter :: epsilon = 1e-20 ! used for the small imaginary part
real, parameter :: zero = 1e-15 ! values smaller than this count as zero
real, parameter :: pi = 3.14159265358979323846
......
:Evaluate:
gpl`args2r[a_]:=Re[N[a/.SubPlus|SubMinus->Identity]];
gpl`args2i[a_]:=Im[N[a/.SubPlus|SubMinus->Identity]];
gpl`args2e[a_]:=Switch[Head[#], SubPlus, 1, SubMinus, -1, _, 1]& /@ a;
:Begin:
:Function: gpl
:Pattern: G[a__]/;And @@ (NumberQ /@ ({a} /. SubPlus | SubMinus -> Identity))
:Arguments: {Re@N[{a}], Im@N[{a}]}
:ArgumentTypes: {RealList,RealList}
:Pattern: gG[a__]/;And @@ (NumberQ /@ ({a} /. SubPlus | SubMinus -> Identity))
:Arguments: { gpl`args2r[{a}], gpl`args2i[{a}], gpl`args2e[{a}] }
:ArgumentTypes: {RealList,RealList,IntegerList}
:ReturnType: Manual
:End:
......@@ -13,17 +17,21 @@
#include <assert.h>
typedef struct {double r,i;} complex;
typedef struct {complex c; signed char i0;} inum;
extern complex __gpl_module_MOD_g_superflatn(complex*,long*);
extern complex __gpl_module_MOD_g_superflatn(inum*,long*);
void gpl(double * re, long nr, double * im, long ni)
void gpl(double * re, long nr, double * im, long ni, int*ieps, long ne)
{
assert(nr==ni);
complex input[nr], ans;
assert(nr==ne);
inum input[nr];
complex ans;
for(long i=0;i<nr;i++)
{
input[i].r = *(re+i);
input[i].i = *(im+i);
input[i].c.r = *(re+i);
input[i].c.i = *(im+i);
input[i].i0 = *(ieps+i);
}
ans = __gpl_module_MOD_g_superflatn(&input[0],&nr);
......
This diff is collapsed.
......@@ -11,59 +11,27 @@ MODULE ieps
type(inum), parameter :: izero=inum( 0.,di0)
type(inum), parameter :: imone=inum(-1.,di0)
type(inum), parameter :: ione=inum(+1.,di0)
type(inum), parameter :: marker=inum(0.,5)
interface operator (*)
module procedure multinum
end interface operator (*)
interface operator (+)
module procedure addinum
end interface operator (+)
interface operator (-)
module procedure subinum
end interface operator (-)
interface operator (**)
module procedure powinum
end interface operator (**)
interface operator (/)
module procedure divint, divinum
end interface operator (/)
interface abs
module procedure absinum, absinumv
end interface abs
interface log
module procedure loginum
end interface log
interface toinum
module procedure toinum_cmplx, toinum_real, toinum_int
end interface toinum
interface tocmplx
module procedure tocmplxv, tocmplxs
end interface tocmplx
interface real
module procedure realis, realiv
end interface real
interface aimag
module procedure imags, imagv
end interface aimag
CONTAINS
FUNCTION MULTINUM(n1, n2)
implicit none
type(inum), intent(in) :: n1, n2
type(inum) :: multinum
multinum = inum( n1%c*n2%c, int(sign(1._prec,real(n1%c)*n2%i0 + real(n2%c)*n1%i0)) )
END FUNCTION MULTINUM
FUNCTION ADDINUM(n1, n2)
implicit none
type(inum), intent(in) :: n1, n2
type(inum) :: addinum
!TODO: what *is* the sum?
addinum = inum(n1%c + n2%c, n1%i0 )
END FUNCTION ADDINUM
FUNCTION SUBINUM(n1, n2)
implicit none
type(inum), intent(in) :: n1, n2
type(inum) :: subinum
!TODO: what *is* the sum?
subinum = inum(n1%c - n2%c, n1%i0 )
END FUNCTION SUBINUM
FUNCTION ABSINUM(n1)
implicit none
type(inum), intent(in) :: n1
......@@ -78,42 +46,6 @@ CONTAINS
absinumv = abs(n1%c)
END FUNCTION ABSINUMV
FUNCTION POWINUM(n1, m)
implicit none
type(inum), intent(in) :: n1
integer, intent(in) :: m
type(inum) :: powinum
if (aimag(n1%c)<zero) then
powinum = inum( cmplx(real(n1%c)**m,0.), int(sign(1._prec,real(n1%c)**m)) )
else
powinum = inum( n1%c**m, n1%i0 )
endif
END FUNCTION POWINUM
FUNCTION DIVINT(n1, m)
implicit none
type(inum), intent(in) :: n1
integer, intent(in) :: m
type(inum) :: divint
divint = inum( n1%c/m, n1%i0*sign(1,m))
END FUNCTION DIVINT
FUNCTION DIVINUM(n1, n2)
implicit none
type(inum), intent(in) :: n1, n2
type(inum) :: divinum
divinum = inum( n1%c/n2%c, int(sign(1., real(n2%c)*n1%i0 - real(n1%c)*n2%i0)))
END FUNCTION DIVINUM
FUNCTION LOGINUM(n1)
implicit none
type(inum), intent(in) :: n1
type(inum) :: loginum
loginum = inum( log(n1%c), n1%i0 * int(sign(1._prec, real(n1%c))) )
END FUNCTION LOGINUM
FUNCTION TOINUM_cmplx(z, s)
complex(kind=prec) :: z(:)
type(inum) :: toinum_cmplx(size(z))
......@@ -129,7 +61,7 @@ CONTAINS
toinum_cmplx(i) = inum(z(i), ss)
enddo
END FUNCTION TOINUM_cmplx
FUNCTION TOINUM_real(z, s)
real(kind=prec) :: z(:)
type(inum) :: toinum_real(size(z))
......@@ -163,5 +95,38 @@ CONTAINS
enddo
END FUNCTION TOINUM_int
FUNCTION TOCMPLXv(z)
type(inum) :: z(:)
complex(kind=prec) tocmplxv(size(z))
tocmplxv = z%c
END FUNCTION
FUNCTION TOCMPLXs(z)
type(inum) :: z
complex(kind=prec) tocmplxs
tocmplxs = z%c
END FUNCTION
FUNCTION REALIV(z)
type(inum) :: z(:)
real(kind=prec) realiv(size(z))
realiv = real(z%c)
END FUNCTION
FUNCTION REALIS(z)
type(inum) :: z
real(kind=prec) realis
realis = real(z%c)
END FUNCTION
FUNCTION IMAGV(z)
type(inum) :: z(:)
real(kind=prec) imagv(size(z))
imagv = aimag(z%c)
END FUNCTION
FUNCTION IMAGS(z)
type(inum) :: z
real(kind=prec) imags
imags = aimag(z%c)
END FUNCTION
END MODULE IEPS
......@@ -3,6 +3,9 @@ MODULE maths_functions
use globals
use utils
implicit none
interface polylog
module procedure polylog1, polylog2
end interface polylog
CONTAINS
FUNCTION zeta(n)
......@@ -120,19 +123,15 @@ CONTAINS
END FUNCTION Li2
RECURSIVE FUNCTION dilog(x) result(res)
! evaluates dilog for any argument
! evaluates dilog for any argument |x|<1
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
else
res = naive_polylog(2,x)
endif
if(abs(aimag(x)) < zero ) then
res = Li2(real(x))
else
res = -dilog(1/x) - (pi**2) /6 - log(add_ieps(-x))**2 / 2
end if
res = naive_polylog(2,x)
endif
END FUNCTION dilog
FUNCTION Li3(x)
......@@ -240,18 +239,14 @@ CONTAINS
END FUNCTION Li3
FUNCTION trilog(x) result(res)
! evaluates trilog for any argument
! evaluates trilog for any argument |x|<1
complex(kind=prec) :: res
complex(kind=prec) :: x
if(abs(x) <= 1.0) then
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
else
res = naive_polylog(3,x)
endif
if(abs(aimag(x)) < zero ) then
res = Li3(real(x))
else
res = naive_polylog(3,sub_ieps(x)**(-1)) - (log(-sub_ieps(x)))**3/6 - pi**2/6 * log(-sub_ieps(x))
end if
res = naive_polylog(3,x)
endif
END FUNCTION trilog
FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res)
......@@ -296,30 +291,56 @@ CONTAINS
END FUNCTION
FUNCTION polylog(m,x) result(res)
RECURSIVE FUNCTION polylog1(m,x) result(res)
! computes the polylog
integer :: m
complex(kind=prec) :: x,res
type(inum) :: x, inv
complex(kind=prec) :: res
if(verb >= 70) print*, 'called polylog(',m,',',x,')'
if ((m.le.9).and.(abs(x-1.).lt.zero)) then
if(verb >= 70) print*, 'called polylog(',m,',',x%c,x%i0,')'
if ((m.le.9).and.(abs(x%c-1.).lt.zero)) then
res = zeta(m)
else if ((m.le.9).and.(abs(x+1.).lt.zero)) then
return
else if ((m.le.9).and.(abs(x%c+1.).lt.zero)) then
res = -(1. - 2.**(1-m))*zeta(m)
else if(m == 2) then
res = dilog(x)
return
else if (abs(x) .gt. 1) then
inv = inum(1./x%c, x%i0)
res = (-1)**(m-1)*polylog(m,inv) &
- cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*conjg(log(-x%c))/2/pi) / factorial(m)
return
endif
if(m == 2) then
res = dilog(x%c)
else if(m == 3) then
res = trilog(x)
res = trilog(x%c)
else
if (abs(x).gt.1) then
res = (-1)**(m-1)*naive_polylog(m,1./x) &
- cmplx(0,2*pi)**m * bernoulli_polynomial(m, 0.5-cmplx(0.,1.)*clog(-x)/2/pi) / factorial(m)
else
res = naive_polylog(m,x)
endif
res = naive_polylog(m,x%c)
end if
END FUNCTION polylog
END FUNCTION polylog1
RECURSIVE FUNCTION polylog2(m,x,y) result(res)
type(inum) :: x, y
integer m
complex(kind=prec) :: res
!TODO!!
res=polylog1(m,inum(x%c/y%c,di0))
END FUNCTION POLYLOG2
FUNCTION PLOG1(a,b)
! calculates log(1-a/b)
implicit none
type(inum) :: a,b
complex(kind=prec) plog1
!TODO!!
plog1 = log(1.-a%c/b%c)
END FUNCTION
END MODULE maths_functions
......
This diff is collapsed.
MODULE GLOBAL_DEF
implicit none
!! ----------
!! parameters
!! ----------
integer, parameter :: prec = selected_real_kind(15,32)
real (kind=prec), parameter :: cw = 0.876613
real (kind=prec), parameter :: sw = 0.481196
real (kind=prec), parameter :: pi = 3.14159265358979323846_prec
real (kind=prec), parameter :: z3 = 1.20205690315959428540_prec
real (kind=prec), parameter :: log2 = 0.693147180559945309417_prec
real (kind=prec), parameter :: conv = 3.893850E+8 ! convert GeV to pb
real (kind=prec), parameter :: xsc = 0._prec ! FDH=>1 vs HV=>0
real (kind=prec), parameter :: Nc = 3._prec
real (kind=prec), parameter :: Tf = 0.5_prec
real (kind=prec), parameter :: Cf = (Nc**2-1)/(2*Nc)
real (kind=prec), parameter :: Nh = 1._prec
real (kind=prec), parameter :: Nf = 5._prec
complex (kind=prec), parameter :: imag = (0.0_prec,1.0_prec)
real (kind=prec), parameter :: zero = 1.0E-50_prec
real (kind=prec), parameter :: alpha_ew = 0.03394_prec
! real (kind=prec), parameter :: alpha = 1/127.9_prec
! real (kind=prec), parameter :: alpha = 1./137.0359997_prec
! real (kind=prec), parameter :: GF = 1.16637E-11_prec ! MeV^-2
real (kind=prec), parameter :: GF = 1._prec
real (kind=prec), parameter :: alpha = 1._prec
real (kind=prec), parameter :: Mmu = 105.658372_prec ! MeV
real (kind=prec), parameter :: Mel = 0.51099893_prec ! MeV
! real (kind=prec), parameter :: Mel = 10._prec ! MeV
real (kind=prec), parameter :: Mtau = 1776.82_prec ! MeV
real (kind=prec), parameter :: xi_sep = 1.0E-10_prec
real (kind=prec), parameter :: del_sep = 1.0E-10_prec
! character (len=3), parameter :: cgamma = "exp"
character (len=3), parameter :: cgamma = "gam"
integer print_ok, throw_away
!! ---------
!! variables
!! ---------
integer :: ran_seed = 1
real (kind=prec) :: p1(4),p2(4),p3(4),p4(4),p5(4),p6(4),p7(4), &
p8(4),p9(4), pol1(4)
real (kind=prec) :: mu, musq, delcut, xinormcut
real (kind=prec) :: xinormcut1, xinormcut2
character (len=8) :: flavour
real (kind=prec) :: Mm ! MeV
real (kind=prec) :: Me ! MeV
contains
SUBROUTINE CRASH(function_name)
character(len=*) :: function_name
write(6,*) "Program crashes because of a call to the function ", &
function_name
stop
END SUBROUTINE CRASH
END MODULE GLOBAL_DEF
! An implementation of the shuffle algebra
! in accordance with 1904.07279v1, polylogs for the masses, p.7-8
! This implementation defines words as strings of characters and shuffles them
! into sums of words.
PROGRAM shuffle_algebra
implicit none
! Currently words can be no longer than the following values
! Might need to be adjusted
integer, parameter :: max_word_size = 6
integer, parameter :: max_word_sum_size = 1000
type word
character(len=max_word_size) :: letters
integer :: length
endtype word
type word_sum
type(word), dimension(max_word_sum_size) :: words
integer :: length
end type word_sum
type(word) :: v1 = word("abc",3)
type(word) :: v2 = word("123",3)
type(word) :: w1
type(word_sum) :: ws, ws1, ws2
ws = shuffle_product(v1,v2)
print*, ws%words
CONTAINS
RECURSIVE FUNCTION shuffle_product(v1, v2) result(res)
! takes two words and returns shuffle product as a word sum
type(word_sum) :: res, p1, p2, s1, s2
type(word) :: v1,v2,w1,w2
type(character) :: alpha, beta
! print*, '----------------------'
! print*, 'v1 = ', v1
! print*, 'v2 = ', v2
if(v1%length == 0) then
res = word_sum((/ v2 /),1)
else if(v2%length == 0) then
res = word_sum((/ v1 /),1)
else
alpha = v1%letters(1:1)
beta = v2%letters(1:1)
w1 = word(v1%letters(2:),v1%length-1)
w2 = word(v2%letters(2:),v2%length-1)
p1 = shuffle_product(w1,v2)
p2 = shuffle_product(v1,w2)
s1 = times(alpha, p1)
s2 = times(beta, p2)
res = combined_word_sum(s1,s2)
end if
END FUNCTION shuffle_product
FUNCTION combined_word_sum(s1, s2)
! combines word sums s1 and s2 into s
type(word_sum) :: s1, s2, combined_word_sum
type(word), dimension(s1%length + s2%length) :: combined_words
integer :: length
length = s1%length + s2%length
combined_words(1:s1%length) = s1%words(1:s1%length)
combined_words(s1%length+1:length) = s2%words(1:s2%length)
combined_word_sum = word_sum(combined_words,length)
END FUNCTION combined_word_sum
FUNCTION times(l, ws)
! computes word sum from letter times word sum, e.g. a(bc + cd) = abc + acd
character :: l
type(word_sum) :: ws
type(word_sum) :: times
integer :: i
times%length = ws%length
do i=1,ws%length
times%words(i)%letters = l // ws%words(i)%letters
end do
END FUNCTION times
SUBROUTINE print_word_sum(ws)
! prints a sum of word with plusses for easy readibility
integer :: i
type(word_sum) :: ws
do i = 1,ws%length
write(*, fmt="(1xai0)", advance="no") ws%words(i)%letters
if(i /= ws%length) then
write(*, fmt="(1xai0)", advance="no") " + "
end if
end do
END SUBROUTINE print_word_sum
END PROGRAM shuffle_algebra
......@@ -8,9 +8,9 @@ CONTAINS
FUNCTION append_to_each_row(a, m) result(res)
! appends element a to each row of m
complex(kind=prec) :: a, m(:,:)
type(inum) :: a, m(:,:)
integer :: i
complex(kind=prec) :: res(size(m,1),size(m,2)+1)
type(inum) :: res(size(m,1),size(m,2)+1)
do i=1,size(m,1)
res(i,:) = [a,m(i,:)]
end do
......@@ -18,21 +18,21 @@ CONTAINS
FUNCTION stack_matrices_vertically(m1, m2) result(res)
! appends to matrix m1 the rows of matrix m2
complex(kind=prec) :: m1(:,:), m2(:,:)
complex(kind=prec) :: res(size(m1,1)+size(m2,1), size(m1,2))
type(inum) :: m1(:,:), m2(:,:)
type(inum) :: res(size(m1,1)+size(m2,1), size(m1,2))
res(1:size(m1,1), :) = m1
res(size(m1,1)+1:size(res,1),:) = m2
END FUNCTION stack_matrices_vertically
RECURSIVE FUNCTION shuffle_product(v1, v2) result(res)
complex(kind=prec) :: v1(:), v2(:)
type(inum) :: v1(:), v2(:)
integer :: i
complex(kind=prec) :: res(product((/(i,i=1,size(v1)+size(v2))/))/ &