module quadruple ! module, to be called from Python for quadruple precision
! calculations
! Universal definition for high precision falling back to double if not
! supported by the compiler:
!integer, parameter :: qk = &
! (abs(min(max(selected_real_kind(18), -4), 0) + 2)/2 ) * selected_real_kind(18) + &
! (abs(min(max(selected_real_kind(14), -4), 0) + 2)/2 - &
! abs(min(max(selected_real_kind(18), -4), 0) + 2)/2 ) * selected_real_kind(14)
!integer, parameter :: rk = selected_real_kind(15,307) ! double
integer, parameter :: qk=16
integer, parameter :: rk=8
real(kind=qk),dimension(:),allocatable,save :: fi,fi_prim
real(kind=qk),dimension(:,:),allocatable,save :: PHI
real(kind=qk),save :: c_star=0.0_qk,ny,lambda_0=0.0_qk,lambda_1=0.0_qk
integer,save :: r_0,r_1,n,m,m_0,m_1
real(kind=qk),dimension(0:4),save :: alphap
real(kind=qk),dimension(:,:),allocatable,save :: beta_ij ! quadruple precision array
real(kind=qk),dimension(:,:),allocatable,save :: a_ij ! quadruple precision array
contains
subroutine calculate_f_fii(r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
integer :: nf
real(kind=qk),dimension(:),allocatable :: ffff ! long double
real(kind=qk), parameter :: pi = 3.14159265358979323846264338327950_qk
nf=size(fi)
allocate(ffff(0:nf-1))
ffff=0.0_qk
ffff=1.-pi/2.0_qk-2.*fi**0.5_qk-2.0_qk*(1.0_qk-fi)**0.5_qk &
-fi*log(1.0_qk+(1.0_qk-fi)**0.5_qk)-(1.0_qk-fi)*log(1.0_qk+fi**0.5_qk) &
+0.5_qk*fi*log(fi)+0.5_qk*(1.0_qk-fi)*log(1.0_qk-fi)
!write(*,*)'ffff=',ffff
r1(1:nf)=ffff(0:nf-1)
end subroutine calculate_f_fii
subroutine calculate_known_sol_fii(r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
integer :: nf
real(kind=qk),dimension(:),allocatable :: ffff ! long double
nf=size(fi)
allocate(ffff(0:nf-1))
ffff=0.0_qk
ffff=1.0_qk+fi**0.5_qk+(1.0_qk-fi)**0.5_qk
r1(1:nf)=ffff(0:nf-1)
end subroutine calculate_known_sol_fii
subroutine calculate_c_star(r_0_in,r_1_in,n_in,m_in,m_0_in,m_1_in,ny_in,c_star_out)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: r_0_in,r_1_in,n_in,m_in,m_0_in,m_1_in
real(kind=rk) :: ny_in
real(kind=rk) :: c_star_out(*)
real(kind=qk) :: c,sgn
integer :: k
!write(*,*)'before c_star=',c_star
r_0=r_0_in
r_1=r_1_in
n=n_in
m=m_in
m_0=m_0_in
m_1=m_1_in
ny=ny_in
c=0.0_qk
sgn=-1.0_qk
do k=0,r_1-1
sgn=-sgn
c=c+1.*sgn/(r_0+k)*comb(r_1-1,k)
enddo
c_star_out(1)=c
c_star=c
!write(*,*)'c_star=',c_star
end subroutine calculate_c_star
subroutine calculate_alpha_p()
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
if (m==4) then
alphap(0)=1.0_qk/36.0_qk !-2,4
alphap(1)=-5.0_qk/18.0_qk !-1,4
alphap(2)=3.0_qk/2.0_qk ! 0,4
alphap(3)=alphap(1) ! 1,4
alphap(4)=alphap(0) ! 2,4
else
print *,"m!=4 not considered yet..."
endif
end subroutine calculate_alpha_p
function comb(mm,i) result(res)
integer :: mm,i,lug,nim,res,k
lug=1
do k=mm-i+1,mm
lug=lug*k
enddo
nim=1
do k=1,i
nim=nim*k
enddo
res=lug/nim
end function comb
subroutine calculate_fi(nt,t_in,r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: nt
real(kind=rk),dimension(*) :: t_in
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
real(kind=qk) :: sgn
real(kind=qk),dimension(:),allocatable :: tt,ff,t
integer :: k
if (allocated(fi)) deallocate(fi)
allocate(fi(0:nt-1))
allocate(tt(0:nt-1))
allocate(ff(0:nt-1))
allocate(t(0:nt-1))
t(0:nt-1)=t_in(1:nt)
sgn=-1.0_qk
ff=0.0_qk
tt=1.0_qk
do k=0,r_1-1
sgn=-sgn
ff=ff+1.0_qk*sgn/(r_0+k)*comb(r_1-1,k)*tt
if (k<r_1-1) then
tt=tt*t
endif
enddo
fi=ff/c_star*t**r_0
!write(*,*)'fi=',fi
deallocate(t)
deallocate(ff)
deallocate(tt)
r1(1:nt)=fi(0:nt-1)
end subroutine calculate_fi
subroutine calculate_fi_prim(nt,t_in)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: nt
real(kind=rk),dimension(*) :: t_in
real(kind=qk),dimension(:),allocatable :: t
if (allocated(fi_prim)) deallocate(fi_prim)
allocate(fi_prim(0:nt-1))
allocate(t(0:nt-1))
t(0:nt-1)=t_in(1:nt)
fi_prim=t**(r_0-1)*(1-t)**(r_1-1)/c_star
!write(*,*)'fi_prim=',fi_prim
deallocate(t)
end subroutine calculate_fi_prim
subroutine calculate_PHI(nt,ns,t_in,s_in,r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: nt,ns
real(kind=rk),dimension(*) :: t_in,s_in
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
real(kind=qk),dimension(:),allocatable :: t,s
integer :: i,j,k
if (allocated(PHI)) deallocate(PHI)
allocate(PHI(0:nt-1,0:ns-1))
allocate(t(0:nt-1))
allocate(s(0:ns-1))
t(0:nt-1)=t_in(1:nt)
s(0:ns-1)=s_in(1:ns)
do i=0,nt-1
do j=0,ns-1
if (i==j) then
PHI(i,j)=fi_prim(i)
else
PHI(i,j)=(fi(i)-fi(j))/(t(i)-s(j))
if (PHI(i,j)<0.0_qk) then ! IMPORTANT! to avoid nan-s in Acurly calculations...
!print "FII<0: FII,fii_i, fii_j, t_i, s_j:",r[i,j],fii[i],fii[j],t[i],s[j]
PHI(i,j)=0.0_qk
endif
endif
enddo
enddo
!write(*,*)'PHI=',PHI
deallocate(s)
deallocate(t)
k=0
do i=0,nt-1
do j=0,ns-1
k=k+1
r1(k)=PHI(i,j)
enddo
enddo
end subroutine calculate_PHI
subroutine calculate_Acurly(nt,ns,t_in,s_in,r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: nt,ns
real(kind=rk),dimension(*) :: t_in,s_in
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
real(kind=qk),dimension(:),allocatable :: t,s
integer :: i,j,k
allocate(a_ij(0:nt-1,0:ns-1))
allocate(t(0:nt-1))
allocate(s(0:ns-1))
t(0:nt-1)=t_in(1:nt)
s(0:ns-1)=s_in(1:ns)
a_ij=0.0_qk
do i=0,nt-1
if (t(i)>=0.0_qk.and.t(i)<=1.0_qk) then
do j=0,ns-1
if (s(j)>=0.0_qk.and.s(j)<=1.0_qk) then
if (PHI(i,j)/=0.0_qk) then
a_ij(i,j)=PHI(i,j)**(-ny)*fi(j)**(-lambda_0)* &
(1.0_qk-fi(j))**(-lambda_1)*fi_prim(j)
endif
endif
enddo
endif
enddo
deallocate(s)
deallocate(t)
k=0
do i=0,nt-1
do j=0,ns-1
k=k+1
r1(k)=a_ij(i,j)
enddo
enddo
!deallocate(a_ij)
end subroutine calculate_Acurly
subroutine calculate_bet_ij(idm,jdim,ii_in,jjm_in,m,n,ny,g1) !: # (54)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: idm,jdim
real(kind=rk),dimension(*) :: ii_in,jjm_in
real(kind=rk),dimension(:),allocatable :: ii,jjm
integer :: m,n,ln
real(kind=rk) :: ny
real(kind=rk),dimension(*) :: g1 ! result as 1-dimensional array
integer :: k,j,i
real(kind=qk),dimension(0:m-1) :: fakt1,fakt4,mult
real(kind=qk) :: m1m,ff,h,h1
allocate(beta_ij(0:idm-1,0:jdim-1))
allocate(ii(0:idm-1))
allocate(jjm(0:jdim-1))
ii(0:idm-1)=ii_in(1:idm)
jjm(0:jdim-1)=jjm_in(1:jdim)
beta_ij=0.0_qk
m1m=(-1._qk)**m
! Prepare the multiplier for the row 1 and 4:
fakt1=0.0_qk
fakt4=0.0_qk
mult=1.0_qk
fakt1(0)=m1m
fakt4(0)=1.0_qk
mult(m-1)=1.0_qk-ny
do k=1,m-1
fakt1(k)=-fakt1(k-1)*k
fakt4(k)=fakt4(k-1)*k
mult(m-k-1)=mult(m-k)*(k+1.0_qk-ny)
enddo
!write (*,*)'fakt1=',fakt1
!write (*,*)'fakt4=',fakt4
!write (*,*)'mult=',mult
do j=0,jdim-1
if (jjm(j)<0.0_rk) then ! : # 1st row
do k=0,m-1
beta_ij(:,j)=beta_ij(:,j)+(ii(:)+m/2.0_qk)**(m-k-ny)/fakt1(k)*jjm(j)**k/mult(k)
enddo
elseif (jjm(j)>n) then ! 4th row
do k=0,m-1
beta_ij(:,j)=beta_ij(:,j)+(n-ii(:)-m/2.0_qk)**(m-k-ny)/fakt4(k)*(jjm(j)-n)**k/mult(k)
enddo
else
do i=0,idm-1
if (jjm(j)<ii(i)+m/2) then !: # 2nd row
ff=m1m
do k=1,m ! # 1,2,...,m
ff=ff/(k-ny)
enddo
beta_ij(i,j)=ff*(ii(i)-jjm(j)+m/2.0_qk)**(m-ny)
else !: # 3rd row
ff=1.0_qk
do k=1,m ! # 1,2,...,m
ff=ff/(k-ny)
enddo
beta_ij(i,j)=ff*(jjm(j)-ii(i)-m/2.0_qk)**(m-ny)
endif
enddo
endif
enddo
ln=jdim
do k=1,m
ln=ln-1
do j=0,ln-1
beta_ij(:,j)=beta_ij(:,j+1)-beta_ij(:,j)
enddo
enddo
h=1.0_qk/n
h1=h**(1.0_qk-ny)
beta_ij=h1*beta_ij
deallocate(jjm)
deallocate(ii)
!deallocate(beta_ij)
k=0
do i=0,idm-1
do j=0,jdim-m-1
k=k+1
g1(k)=beta_ij(i,j)
enddo
enddo
end subroutine calculate_bet_ij
subroutine calculate_tau(ni,nj,beta_j0_in,b_ij_in,r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: ni,nj
real(kind=rk),dimension(*) :: beta_j0_in,b_ij_in
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
real(kind=qk),dimension(:),allocatable :: beta_j0
real(kind=rk),dimension(:,:),allocatable :: b_ij
real(kind=qk),dimension(:,:),allocatable :: tau
integer :: k,j,i
allocate(b_ij(0:nj-1,0:ni-1))
!b_ij=reshape(b_ij_in(1:ni*nj),(/nj,ni/) ) ! we get this in trasposed form
k=0
do i=0,ni-1
do j=0,nj-1
k=k+1
b_ij(i,j)=b_ij_in(k)
enddo
enddo
allocate(beta_j0(0:ni-1))
beta_j0(0:ni-1)=beta_j0_in(1:ni)
allocate(tau(0:ni-1,0:nj-1))
tau=0.0_qk
do k=0,ni-1
do j=k,k+2*(m_1-1)
tau(:,k)=tau(:,k)+alphap(j-k)*(beta_ij(:,j)*a_ij(:,k)+beta_j0(k)*b_ij(k,:))
!if (k==0) then
! print *,'tau(0,0)=',tau(0,0),beta_ij(0,j),a_ij(0,0),beta_j0(0),b_ij(0,0)
!endif
enddo
enddo
deallocate(beta_j0)
deallocate(b_ij)
deallocate(a_ij)
deallocate(beta_ij)
k=0
do i=0,ni-1
do j=0,nj-1
k=k+1
r1(k)=tau(i,j)
enddo
enddo
deallocate(tau)
end subroutine calculate_tau
subroutine calculate_delta_m_gamjt(r1)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
real(kind=rk),dimension(*) :: r1 ! result as 1-dimensional array
real(kind=qk),dimension(:),allocatable :: tt,ff,t
real(kind=qk) :: m1m
integer,dimension(:),allocatable :: jjj
integer :: i,k,len1,len2
len1=n
allocate(jjj(0:len1))
!jjj=linspace(0.,n,n+1)
jjj=(/ (i,i=0,len1) /)
allocate(ff(0:len1))
ff=1.0_qk
m1m=(-1)**m
ff(0:m/2)=m1m*(m/2.0_qk-jjj(0:m/2))**(m-ny)
ff(m/2+1:len1)=(jjj(m/2+1:len1)-m/2.0_qk)**(m-ny)
do i = 1,m
do k=0,len1-i
ff(k)=ff(k+1)-ff(k)
enddo
enddo
!write(*,*)'ff=',ff
len2=len1-m+1
r1(1:len2)=ff(0:len2-1)
deallocate(ff)
deallocate(jjj)
end subroutine calculate_delta_m_gamjt
!subroutine Calculate_D2_m(w,m)
! implicit none
! integer, parameter :: qk=16
! real(kind=qk),dimension(:,:),intent(inout) :: w
! integer :: m
! integer :: ln,k
! ln=size(w,dim=2)+1
! do k=1,m
! ln=ln-1
! call Calculate_D2(w,ln)
! enddo
!end subroutine Calculate_D2_m
!
!subroutine Calculate_D2(w,ln)
! implicit none
! integer, parameter :: qk=16
! real(kind=qk),dimension(0:ln-1,0:ln-1),intent(inout) :: w
! integer :: ln,j
! do j=0,ln-1
! w(:,j)=w(:,j)-w(:,j+1)
! enddo
!end subroutine Calculate_D2
subroutine test(n,ii,gg)
implicit none
integer, parameter :: qk=16
integer, parameter :: rk=8
integer :: n
real(kind=rk), dimension(*) :: ii
!!double precision, dimension(:), intent(inout) :: ii ! töötab aga size 0
!double precision, dimension(:), pointer :: ii ! töötab aga size 0
!double precision, dimension(*) :: ii ! on OK!
!real(kind=rk),dimension(size(ii)) :: gg ! double precision array
!!double precision,dimension(size(ii)) :: gg ! double precision array
double precision,dimension(*) :: gg ! double precision array
integer :: i,j
double precision,dimension(0:n-1,0:n-1) :: ri,rg
ri= reshape(ii(1:n*n),(/n,n/) )
do j=0,n-1
do i=0,n-1
rg(i,j)=n*(j)+1.0D0*i+1
enddo
enddo
!write(*,*)'size(ii) is:',size(ii)
write(*,*)'ri is:',ri(0:n-1,0:n-1)
gg(1:n*n)=reshape(rg,(/n*n/))
return
end subroutine test
end module quadruple