-
Notifications
You must be signed in to change notification settings - Fork 2
/
vijkl2.jl
127 lines (104 loc) · 4.98 KB
/
vijkl2.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
function recurringA0(M2::Int64,m1::Int64,m2::Int64,m3::Int64,m4::Int64)
# return (M-m1-m2-m3-m4+4)/(2M-2m1-2m2-2m3-2m4+8)/(2M-2m1-2m2-2m3-2m4+7)
# return BigFloat(1)/BigFloat(M2-2m1-2m2-2m3-2m4+7) #BigFloat(1/(M2-2m1-2m2-2m3-2m4+7))
return 1//BigInt(M2-2m1-2m2-2m3-2m4+7)
end
function recurringA1(n1::Int64,n2::Int64,n3::Int64,n4::Int64,m1::Int64,m2::Int64,m3::Int64,m4::Int64)
# M = Int64((n1 + n2 + n3 + n4)/2) # integer
if m1 == m2 == m3 == 0
# coeff = (-1)*BigFloat(n4-2m4+2)*BigFloat(n4-2m4+1)/BigFloat(m4)*recurringA0(n1+n2+n3+n4,1,1,1,m4)
coeff = (-1)*BigInt(n4-2m4+2)*BigInt(n4-2m4+1)//BigInt(m4)*recurringA0(n1+n2+n3+n4,1,1,1,m4)
elseif m1 == m2 == 0
# coeff = (-1)*BigFloat(n3-2m3+2)*BigFloat(n3-2m3+1)/BigFloat(m3)*recurringA0(n1+n2+n3+n4,1,1,m3,m4)
coeff = (-1)*BigInt(n3-2m3+2)*BigInt(n3-2m3+1)//BigInt(m3)*recurringA0(n1+n2+n3+n4,1,1,m3,m4)
elseif m1 == 0
# coeff = (-1)*BigFloat(n2-2m2+2)*BigFloat(n2-2m2+1)/BigFloat(m2)*recurringA0(n1+n2+n3+n4,1,m2,m3,m4)
coeff = (-1)*BigInt(n2-2m2+2)*BigInt(n2-2m2+1)//BigInt(m2)*recurringA0(n1+n2+n3+n4,1,m2,m3,m4)
else
# coeff = (-1)*BigFloat(n1-2m1+2)*BigFloat(n1-2m1+1)/BigFloat(m1)*recurringA0(n1+n2+n3+n4,m1,m2,m3,m4)
coeff = (-1)*BigInt(n1-2m1+2)*BigInt(n1-2m1+1)//BigInt(m1)*recurringA0(n1+n2+n3+n4,m1,m2,m3,m4)
end
return coeff
end
function Vijkl2(n1::Int64,n2::Int64,n3::Int64,n4::Int64)
if isodd(n1 + n2 + n3 + n4)
return 0.0
end
M = Int64((n1 + n2 + n3 + n4)/2)
maxm1 = floor(Int64,n1/2)+1
maxm2 = floor(Int64,n2/2)+1
maxm3 = floor(Int64,n3/2)+1
maxm4 = floor(Int64,n4/2)+1
matA = zeros(Rational{BigInt},2,maxm2,maxm3,maxm4) #zeros(BigFloat,2,maxm2,maxm3,maxm4) #zeros(Float64,2,maxm2,maxm3,maxm4)
# m1=m2=m3=m4=1
matA[1,1,1,1] = Rational(BigInt(1)) #BigFloat(1.0)
# matA[1,1,1,1] = exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))
# m1=m2=m3=1
for m4 = 1:maxm4-1
# matA[1,1,1,m4+1] = matA[1,1,1,m4]*(-1)*2*(n4-2m4+2)*(n4-2m4+1)/m4*(M-m4+1)/(2M-2m4+2)/(2M-2m4+1)
matA[1,1,1,m4+1] = matA[1,1,1,m4]*recurringA1(n1,n2,n3,n4,0,0,0,m4)
end
# m1=m2=1
for m3 = 1:maxm3-1
# Threads.@threads for m4 = 1:maxm4 # parfor
for m4 = 1:maxm4
# matA[1,1,m3+1,m4] = matA[1,1,m3,m4]*(-1)*2*(n3-2m3+2)*(n3-2m3+1)/m3*(M-m3-m4+2)/(2M-2m3-2m4+4)/(2M-2m3-2m4+3)
matA[1,1,m3+1,m4] = matA[1,1,m3,m4]*recurringA1(n1,n2,n3,n4,0,0,m3,m4)
end
end
# m1=1
for m2 = 1:maxm2-1
# Threads.@threads for m34 = 1:maxm3*maxm4 # parfor
for m34 = 1:maxm3*maxm4
m3 = div(m34,maxm4)+1
m4 = mod(m34,maxm4) # m34-div(m34,maxm4)*maxm4
if m4 == 0
m3 = m3-1 # div(m34,maxm4)
m4 = maxm4 # m34-(div(m34,maxm4)-1)*maxm4 = m34-(m3-1)*maxm4 = maxm4
end
# matA[1,m2+1,m3,m4] = matA[1,m2,m3,m4]*(-1)*2*(n2-2m2+2)*(n2-2m2+1)/m2*(M-m2-m3-m4+3)/(2M-2m2-2m3-2m4+6)/(2M-2m2-2m3-2m4+5)
matA[1,m2+1,m3,m4] = matA[1,m2,m3,m4]*recurringA1(n1,n2,n3,n4,0,m2,m3,m4)
end
end
sumA = sum(matA[1,:,:,:])
for m1 = 1:maxm1-1
# Threads.@threads for m234 = 1:maxm2*maxm3*maxm4 # parfor
for m234 = 1:maxm2*maxm3*maxm4
m2 = div(m234,maxm3*maxm4)+1
m3 = div(m234 - (m2-1)*maxm3*maxm4,maxm4)+1
m4 = m234 - (m2-1)*maxm3*maxm4 - (m3-1)*maxm4
if mod(m234,maxm3*maxm4) == 0
m2 = m2-1
m3 = maxm3
m4 = maxm4
elseif mod(m234 - (m2-1)*maxm3*maxm4,maxm4) == 0
m3 = m3-1
m4 = maxm4
end
# matA[2,m2,m3,m4] = matA[1,m2,m3,m4]*(-1)*2*(n1-2m1+2)*(n1-2m1+1)/m1*(M-m1-m2-m3-m4+4)/(2M-2m1-2m2-2m3-2m4+8)/(2M-2m1-2m2-2m3-2m4+7)
matA[2,m2,m3,m4] = matA[1,m2,m3,m4]*recurringA1(n1,n2,n3,n4,m1,m2,m3,m4)
end
sumA = sumA + sum(matA[2,:,:,:])
matA[1,:,:,:] .= matA[2,:,:,:]
matA[2,:,:,:] .= Rational(BigInt(0)) #0
end
# sumA = sumA*exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))
sumA = Float64(sumA)*exp(-2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4))))
# if sumA < 0
# sumA = -sumA
# sumA = log(sumA) -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
# sumA = exp(sumA)
# sumA = -sumA
# else
# sumA = log(sumA) -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
# sumA = exp(sumA)
# end
sumA = sumA/sqrt(2*pi)/2
return sumA
# return matA
end
function A0000(n1::Int64,n2::Int64,n3::Int64,n4::Int64)
M = Int64((n1 + n2 + n3 + n4)/2)
A0 = -2M*log(2)+sum(log.(M+1:2*M))-1/2*(sum(log.(2:n1))+sum(log.(2:n2))+sum(log.(2:n3))+sum(log.(2:n4)))
return exp(A0)/sqrt(2*pi)/2
end