Merge branch 'simple_relative'

helper.jl

@@ -1,7 +1,7 @@
 using LinearAlgebra, DelimitedFiles
 
 "Sum over array while minimizing catastrophic cancellation as much as possible"
-function better_sum(arr::Array{Float64})
+function better_sum(arr::Array{T}) where T<:Real
     pos_arr = arr[arr .> 0]
     neg_arr = arr[arr .< 0]
 

ho_basis.jl

@@ -1,21 +1,8 @@
 using SparseArrays
-using SpecialFunctions
 using QuadGK
 using LRUCache
 include("helper.jl")
-
-# Gaussian potentials in HO space
-inv_factorial(n) = Iterators.prod(inv.(1:n))
-sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
-N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
-Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
-
-# numerical evaluation of V matrix elements
-sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
-sqrt_sqrt_pi = sqrt(sqrt(pi))
-laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
-ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
+include("math.jl")
 
 function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
     integrand(r) = sqrt(μω_gen) * ho_basis(l, n1, sqrt(μω_gen) * r) * ho_basis(l, n2, sqrt(μω_gen) * r) * V_of_r(r)
@@ -23,8 +10,6 @@ function V_numerical(V_of_r, l, n1, n2; μω_gen=1.0, atol=0, maxevals=10^7)
     return integral
 end
 
-##############################################################
-
 function get_sp_basis(E_max)
     Es = Int[]
     ns = Int[]
@@ -175,3 +160,47 @@ function get_jacobi_V2_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxeva
 
     return V2
 end
+
+function get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μ1ω1, μ2ω2; dtype=Float64)
+    mat = zeros(dtype, length(n1s), length(n1s))
+    Threads.@threads for idx in CartesianIndices(mat)
+        (i, j) = Tuple(idx)
+        val = racahs_reduction_formula(n1s[i], l1s[i], n2s[i], l2s[i], n1s[j], l1s[j], n2s[j], l2s[j], Λ, μ1ω1, μ2ω2)
+        if !(val ≈ 0); mat[idx] = val; end
+    end
+    return sparse(mat)
+end
+
+function get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global; atol=10^-6, maxevals=10^5)
+    _, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
+    l_max = max(maximum(l1s), maximum(l2s))
+    n_max = max(maximum(n1s), maximum(n2s))
+    mask1 = (n2s .== n2s') .&& (l2s .== l2s')
+    mask2 = (n1s .== n1s') .&& (l1s .== l1s')
+    
+    V_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω, atol=atol, maxevals=maxevals)
+    V_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
+
+    V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_cache)
+    V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2, dtype=ComplexF64, cache=V_cache)
+
+    V12 =  get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=atol, maxevals=maxevals)
+
+    return V1 + V2 + V12
+end
+
+function get_src_V12_matrix(V_of_r, E_max, Λ, μω_global; atol=10^-6, maxevals=10^5)
+    Es, n1s, l1s, n2s, l2s = get_2p_basis(E_max, Λ)
+    l_max = max(maximum(l1s), maximum(l2s))
+    n_max = max(maximum(n1s), maximum(n2s))
+    mask1 = (n2s .== n2s') .&& (l2s .== l2s')
+    
+    V_relative_elem(l, n1, n2) = V_numerical(V_of_r, l, n1, n2; μω_gen=μω_global, atol=atol, maxevals=maxevals)
+    V_relative_cache = LRU{Tuple{UInt8, UInt8, UInt8}, ComplexF64}(maxsize=(1+l_max)*(1+n_max)^2)
+
+    V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1, dtype=ComplexF64, cache=V_relative_cache)
+    U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
+    V12 =  U' * V_relative * U
+
+    return V12
+end

ho_basis_3body.jl

@@ -7,13 +7,11 @@ Va = -2
 Ra = 2
 
 E_max = 40
-μω_global = 0.5
+μω_global = 0.3
 
-# due to Jacobi coordinates
-μ1ω1 = μω_global * 1/2
-μ2ω2 = μω_global * 2
-μ1 = m * 1/2
-μ2 = m * 2/3
+# due to simple relative coordinates
+μω = μω_global * 2
+μ = m/2
 
 println("No of threads = ", Threads.nthreads())
 
@@ -26,19 +24,21 @@ end
 println("Basis size = ", length(Es))
 
 println("Constructing KE matrices")
-@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
-@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
+@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
+@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
+@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω) ./ (2*μ)
 
 println("Constructing PE matrices")
-V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μ1ω1)
+V_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω)
 V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_global)
-@time "V1" V1 = get_sp_V_matrix(V1_elem, n1s, l1s; mask=mask1)
-@time "V relative" V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1) + get_sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2)
+@time "V1" V1 = get_sp_V_matrix(V_elem, n1s, l1s; mask=mask1)
+@time "V2" V2 = get_sp_V_matrix(V_elem, n2s, l2s; mask=mask2)
+@time "V relative" V_relative = get_sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1)
 @time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
-@time "V2" V2 =  U' * V_relative * U
+@time "V12" V12 =  U' * V_relative * U
 
 println("Calculating spectrum")
-@time "H" H = T1 + T2 + V1 + V2
+@time "H" H = T1 + T2 + T_cross + V1 + V2 + V12
 @time "Eigenvalues" evals, _ = eigs(H, nev=3, ncv=30, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
 
 display(evals)

ho_basis_3body_resonance.jl

@@ -8,11 +8,9 @@ V_of_r(r) = 2 * exp(-(r-3)^2 / (1.5)^2)
 E_max = 30
 μω_global = 0.5 * exp(-2im * pi / 9)
 
-# due to Jacobi coordinates
-μ1ω1 = μω_global * 1/2
-μ2ω2 = μω_global * 2
-μ1 = m * 1/2
-μ2 = m * 2/3
+# due to simple relative coordinates
+μω = μω_global * 2
+μ = m/2
 
 println("No of threads = ", Threads.nthreads())
 
@@ -28,15 +26,16 @@ println("Basis size = ", length(Es))
 
 println("Constructing KE matrices")
 
-@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μ1ω1, μ=μ1)
-@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μ2ω2, μ=μ2)
+@time "T1" T1 = get_sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω, μ=μ)
+@time "T2" T2 = get_sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω, μ=μ)
+@time "T_cross" T_cross = get_2p_p1p2_matrix(n1s, l1s, n2s, l2s, Λ, μω, μω; dtype=ComplexF64) ./ (2*μ)
 
 println("Constructing PE matrices")
 
-@time "V" V = get_jacobi_V_matrix(V_of_r, E_max, Λ, μ1ω1, μω_global)
+@time "V" V = get_src_V_matrix(V_of_r, E_max, Λ, μω, μω_global)
 
 println("Calculating spectrum")
-@time "H" H = T1 + T2 + V
+@time "H" H = T1 + T2 + T_cross + V
 @time "Eigenvalues" evals, _ = eigs(H, nev=5, ncv=50, which=:LI, maxiter=5000, tol=1e-5, ritzvec=false, check=1)
 
 display(evals)

math.jl

@@ -0,0 +1,38 @@
+using SpecialFunctions, WignerSymbols
+include("helper.jl")
+
+# Gaussian potentials in HO space
+inv_factorial(n) = Iterators.prod(inv.(1:n))
+sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
+N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2)
+Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
+V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+
+# for numerical evaluation of V matrix elements
+sqrt_double_factorial(n) = Iterators.prod(sqrt.(n:-2:1))
+sqrt_sqrt_pi = sqrt(sqrt(pi))
+laguerre(l, n, x) = gamma(n + l + 3/2) * better_sum([(-x * x)^k / gamma(k + l + 3/2) * inv_factorial(n - k) * inv_factorial(k)  for k in 0:n])
+ho_basis(l, n, x) = (-1)^n / sqrt_sqrt_pi * 2^((n + l + 2) / 2) * sqrt_factorial(n) / sqrt_double_factorial(2*n + 2*l + 1) * x^(l + 1) * exp(-x^2 / 2) * laguerre(l, n, x)
+
+# for implementation of simple relative coordinates
+double_factorial(n::Int) = Iterators.prod(big, n:-2:1)
+
+"Gaussian integral for n ∈ Integers (Ref: Wolfram MathWorld + simplifications)"
+gauss_int(a, n) = double_factorial(n - 1) / (2 * a)^((n + 1)/2) * (iseven(n) ? sqrt(π / 2) : 1)
+
+"Gives ∫dp p u' u where u' and u' may have different l (Ref: worked out in Mathematica)"
+function integral(np, lp, n, l, μω)
+    s = [(-1)^(m + mp) * gauss_int(1, 2 * m + 2 * mp + l + lp + 3) * N_lnk(l, n, m) * N_lnk(lp, np, mp) for (m, mp) in Iterators.product(0:n, 0:np)]
+    return 2 * sqrt(μω) * better_sum(s)
+end
+
+"Gives <n' l'|| p ||n l> for the HO basis"
+reduced_matrix_element(np, lp, n, l, μω) = (-1)^lp * sqrt(2*lp + 1) * sqrt(2*l + 1) * wigner3j(Float64, lp, 1, l, 0, 0, 0) * integral(np, lp, n, l, μω)
+
+"Matrix element <n1p l1p n2p l2p| p1⋅p2 |n1 l1 n2 l2> (Ref: de-Shalit & Talmi, Eq 15.5)"
+function racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)
+    val = wigner6j(Float64, l1p, l2p, Λ, l2, l1, 1)
+    if val != 0; val *= reduced_matrix_element(n1p, l1p, n1, l1, μ1ω1); end
+    if val != 0; val *= reduced_matrix_element(n2p, l2p, n2, l2, μ2ω2); end
+    return (-1)^(l1 + l2p + Λ) * val
+end

test/racah_test_data.csv

@@ -0,0 +1,384 @@
+2,0,2,3,7,3,5,1,3,2.2274833794413214,0.45049378982090493,0.
+8,1,3,3,2,2,0,0,2,2.202162451867089,2.3467827533568713,0.
+4,1,7,1,0,2,0,3,2,1.1507774013957701,1.7553231930689552,0.
+3,3,7,1,7,0,3,3,3,1.7207018027964498,0.5156168080476098,0.
+10,3,2,3,2,2,8,1,2,1.8402526798628767,2.8981293881514993,0.
+4,3,8,1,0,3,0,3,3,2.471579266664868,0.9515272650259292,0.
+9,2,8,3,10,3,1,2,1,2.1675261170354094,0.38546777353112605,0.
+2,2,7,3,4,3,7,3,3,2.76860332907484,0.7916486882071214,0.
+2,3,3,3,5,1,1,0,1,1.6146747022554235,2.8377955838430013,0.
+4,3,5,0,7,3,2,3,3,0.9198951336314427,0.9950080709108984,0.
+7,3,4,2,2,2,8,0,2,0.2364795732399818,1.839193911487496,0.
+7,3,3,3,5,1,6,2,3,0.29514953289458434,0.3780481641738338,0.
+5,2,4,1,9,3,1,0,3,2.8086957499249543,2.1260743936378903,0.
+4,1,6,1,1,0,10,0,0,2.6628402589565523,2.425950546585856,0.
+8,3,1,0,6,3,0,3,3,2.1328726647005265,0.9431361454676961,0.
+4,3,5,0,5,1,0,3,3,2.576495013458395,0.37293785437117544,0.
+10,3,4,0,1,2,2,3,3,2.0030278179077214,1.1187270744696427,0.
+9,1,8,2,1,1,5,3,2,2.789813153068671,2.6763038014457665,0.
+5,0,6,1,5,3,9,3,1,1.9319197502299517,0.6216330166798851,0.
+8,0,6,2,7,3,8,2,2,1.8084358691372344,2.0181750785519634,0.
+3,1,5,1,9,0,7,1,1,1.599191638971793,1.552252944957802,0.
+1,0,7,2,6,3,4,2,2,2.4217191310288078,1.2381047980553532,0.
+9,3,2,0,1,2,4,2,3,1.261927284387386,0.39708064622835604,0.
+0,1,5,1,5,0,0,2,2,0.407989632316895,2.4021451804173477,0.
+3,1,6,3,9,2,2,1,2,1.2053920704239278,0.25375827636999704,0.
+9,2,7,0,5,2,5,3,2,2.539812963473447,0.7382167735661622,0.
+10,2,7,3,4,3,0,0,3,0.8427632747449358,0.7622120453149561,0.
+5,2,0,2,7,3,5,0,3,0.460843656858851,0.583650962079461,0.
+8,1,8,0,6,2,0,2,1,2.056601089239093,0.22262414697020771,0.
+4,2,4,2,6,1,0,3,3,1.4808637079563942,1.358717819227813,0.
+10,3,4,3,4,0,5,2,2,1.5284998408295385,2.925365977036159,0.
+10,1,1,1,5,1,10,1,0,1.1846342029623314,1.186468360104492,0.
+9,0,5,2,1,2,2,3,2,1.733114516642603,2.8283467464705456,0.
+8,3,8,1,10,3,4,3,3,0.8918790824078111,0.9129770894549143,0.
+10,0,9,3,5,3,8,3,3,0.6417078137732939,1.965954991753998,0.
+10,2,8,2,8,2,8,3,4,0.6600917142427227,0.987969186819146,0.
+8,2,4,1,3,2,9,1,3,0.8935573241720247,2.1162172596314806,0.
+4,2,10,3,7,2,7,2,2,1.5052737158592535,2.7171501665361504,0.
+0,1,3,1,5,2,5,0,2,1.577080591236359,0.34164973402053134,0.
+0,2,10,3,6,3,6,0,3,2.7925439641102736,1.831313118698632,0.
+1,0,9,3,9,0,0,3,3,2.3331462089933934,1.6945470525813482,0.
+8,2,7,3,7,2,3,1,1,1.4988059394628408,1.3530552540483156,0.
+9,2,8,1,8,3,6,1,2,0.6337531145346973,1.5268564417495032,0.
+10,1,0,0,7,2,7,2,1,1.3790583258395,1.5669277734464737,0.
+3,0,0,2,6,0,4,2,2,2.4730048216795737,2.70759381187844,0.
+8,1,3,1,5,3,8,2,2,0.8399302875607697,0.32436518041750384,0.
+4,1,1,2,4,3,10,1,3,0.8437127682376069,1.7940889253407324,0.
+9,0,3,1,4,3,2,2,1,0.8755496194878809,0.9675651711989977,0.
+9,0,2,2,3,1,8,2,2,1.4950146011492311,2.8152388609955157,0.
+7,0,8,3,2,3,1,3,3,2.2587603015484863,2.0666085413750848,0.
+10,3,9,1,4,3,0,2,4,1.5489892079797567,1.7151221413658782,0.
+3,0,1,2,2,3,4,1,2,0.25514014111567285,2.06281752109991,0.
+5,1,1,0,9,1,3,2,1,1.3329695370285548,0.5298097905264418,0.
+8,1,1,0,5,0,9,1,1,2.95254173854617,2.1687131876307753,0.
+6,1,8,0,5,2,1,2,1,1.5801597077900533,2.257296177275017,0.
+10,0,9,1,9,3,4,3,1,2.5352867964035397,1.7981404960338239,0.
+10,2,1,2,0,0,9,3,3,2.06463437201601,1.4295850747781111,0.
+1,0,1,3,6,3,6,3,3,0.8872777463742612,1.9002353924129451,0.
+6,3,10,0,2,3,2,1,3,1.2410196118280048,1.8568476669071874,0.
+8,2,1,2,2,3,6,3,2,0.5302021471052729,1.2216709678909177,0.
+1,2,8,1,0,3,3,0,3,2.8876975946711045,2.717322643102454,0.
+4,2,4,1,0,1,7,0,1,2.0026155982539864,0.7640140889966269,0.
+1,0,1,1,4,2,3,2,1,1.016582775534538,2.3841546013418045,0.
+2,0,3,1,10,1,1,1,1,1.6404532954231827,2.1972778830632356,0.
+7,3,1,3,2,2,4,3,4,2.648817905371118,2.798529790102262,0.
+9,0,4,1,6,0,2,1,1,2.3138662599243407,1.9450229178773055,0.
+0,1,8,1,1,1,1,0,1,0.6524221878214105,2.9034613327001972,0.
+6,1,6,3,0,1,6,2,3,1.831532710689828,2.890598375121809,0.
+4,1,0,0,0,2,7,3,1,0.626504505885948,2.6340691505629463,0.
+10,3,2,3,10,3,1,1,2,1.6591641837977482,2.452206213376977,0.
+3,1,6,2,2,2,7,1,3,1.1999874026225852,2.795586000355012,3.3573329331683617
+8,3,1,2,1,3,5,2,5,0.6830778416708307,1.1113967340568016,0.
+7,2,3,1,5,3,8,2,3,1.714181314975816,1.3518067184764888,0.
+0,2,4,2,0,1,4,1,2,1.060932248034836,0.8036155418186053,-1.1371377045290691
+4,3,9,3,10,3,10,3,5,0.42490400067847744,1.9191837942593741,0.
+0,0,6,1,2,2,7,2,1,1.4606282497051195,2.390102981429126,0.
+8,0,0,1,9,0,10,1,1,2.9779138627844803,1.482942802025069,0.
+1,2,2,3,3,3,8,3,4,1.540918746101882,0.9621051195430765,0.
+2,1,1,3,2,0,1,3,3,2.0584800719627525,1.7796724319652748,0.
+10,2,5,2,4,0,0,2,2,1.2759898818923578,0.6080747387544254,0.
+1,0,9,3,2,2,10,2,3,2.835152414163625,1.2470542223023808,0.
+1,2,2,0,6,0,7,2,2,1.8642054302207116,0.48490065380939873,0.
+5,1,10,2,5,0,10,3,3,0.6428303158043476,0.5501263433952537,2.105493662515726
+1,1,1,1,4,2,0,1,1,2.91586289712962,0.2030446107349677,0.
+4,2,0,2,1,3,1,0,3,0.31239590050870536,2.3530015462264577,0.
+10,0,0,3,0,1,4,3,3,1.4066735483371904,0.3492189542356794,0.
+3,1,9,3,6,2,4,2,4,2.039761320004641,2.0376049586912632,0.
+4,1,3,3,2,2,8,2,2,2.485347427066465,2.2425974670112305,0.
+8,0,10,3,7,3,3,0,3,1.145305638226675,0.7430070058758362,0.
+1,0,4,3,1,3,3,2,3,2.236729779510484,2.518871850291589,0.
+9,3,10,2,0,0,9,3,3,2.2231614378407896,2.934256387523143,0.
+1,3,7,1,6,2,7,1,3,1.739591417445606,2.852892536289553,0.
+9,0,0,2,7,1,0,2,2,2.628613403491599,1.948716181798507,0.
+2,0,5,3,3,2,4,3,3,0.7244831235767495,0.307815884940033,0.
+9,2,2,3,4,0,9,2,2,1.4712965684171477,0.7152116770614816,0.
+10,0,7,1,7,2,2,2,1,2.0944053463871253,0.7877074060643912,0.
+10,2,8,2,4,1,4,3,3,1.2400259243830605,1.533481685252715,0.
+7,2,2,0,6,3,2,1,2,0.7570369287424406,1.5727691945398687,2.4154000428935993
+9,0,9,2,4,2,5,0,2,2.9903195228568853,1.7342728120955062,0.
+0,2,2,1,2,2,2,2,3,2.871118456330721,1.8840932443901108,0.
+1,3,4,1,9,2,8,1,3,2.2820457592186356,2.8820550176739523,0.
+4,3,9,1,5,3,7,1,2,0.5157719636123437,2.2186070481213065,0.
+3,3,9,3,7,2,1,1,1,1.806277095079202,2.303389106635813,0.
+4,2,0,0,8,0,10,2,2,0.4100439145977339,1.0520725633546872,0.
+10,1,6,0,4,2,2,3,1,1.126129220884609,1.5550003535454513,0.
+3,2,10,1,6,1,10,1,1,0.2658540542443206,0.608227371386918,0.
+6,3,5,1,4,1,0,3,3,2.705056186004474,0.42463741187699755,0.
+10,0,5,2,5,1,1,2,2,0.8756297484714799,0.7584762754047034,0.
+7,1,7,3,5,1,0,3,4,2.914828464578812,2.261191790401713,0.
+10,0,3,1,10,2,4,3,1,2.50303873714987,1.5051451493521428,0.
+1,1,10,0,1,3,0,3,1,1.3072742851724795,1.9423953478998177,0.
+1,1,3,2,3,2,0,3,2,2.3584710757323544,1.1305990899916512,0.
+10,3,2,3,3,1,6,1,0,1.4654067038271226,1.985658804098854,0.
+7,2,0,2,9,2,0,1,1,1.3873373238347515,0.7703862068685314,0.
+8,3,7,0,8,2,5,1,3,0.8253295318767577,2.899984116328765,0.
+1,1,0,0,0,1,10,2,1,2.3641832402927276,2.1083504299969444,0.
+8,2,7,3,10,1,0,2,2,1.3652741505836983,2.0949097737438462,0.
+8,3,10,0,4,2,5,1,3,0.24020243315960244,1.9022194572611664,0.
+9,1,9,2,0,0,4,1,1,2.5376357730342507,1.8526159863451914,0.
+7,3,9,3,3,3,7,0,3,1.6807661985105629,2.687649747358564,0.
+7,2,2,0,9,2,8,1,2,2.4565884018951802,0.8982169296897196,0.
+9,2,0,3,10,3,5,1,2,0.6353524665165655,1.0371641317201128,0.
+3,0,7,3,2,2,1,1,3,2.7958304824728186,2.4767642915078376,0.
+9,0,7,1,10,1,4,0,1,0.25616825512317254,1.308344122398832,0.
+1,2,1,3,10,1,1,1,2,2.224763269329724,2.231499665082513,0.
+8,0,4,3,0,2,5,3,3,0.770816371156068,2.3322630361929546,0.
+5,1,1,2,2,2,3,1,2,0.8022984805467428,1.4426168425194534,0.
+7,0,5,0,0,2,10,2,0,1.9511722131189702,0.4950560622323814,0.
+9,3,7,1,7,0,10,2,2,2.775167770563364,0.45822632372936667,0.
+7,2,7,3,9,0,6,1,1,2.8779651825740737,2.4436285639972706,0.
+1,1,6,3,0,0,3,3,3,0.5442622594206528,2.9853682012412737,0.
+10,2,8,1,4,2,7,3,2,1.329367888136959,0.9413429526552286,0.
+7,3,0,1,4,1,9,1,2,1.3457614299418887,1.5531678736115984,0.
+4,2,8,2,10,1,3,3,3,0.8248650313839705,0.6616884801498313,0.
+0,0,8,3,3,1,10,2,3,2.0673228373578825,1.4463899838479248,0.
+0,3,0,2,9,3,0,3,5,0.906810420580336,2.811065359991593,0.
+2,1,2,0,4,3,8,2,1,1.4719829696395097,1.5045170591445123,0.
+8,3,4,2,7,0,6,2,2,2.934501293001735,2.475691074854602,0.
+4,1,7,3,9,1,2,2,3,1.5229908742702083,2.7910336359576577,0.
+3,0,4,3,9,2,8,2,3,0.7395750631044415,2.7437950070306316,0.
+0,3,10,0,8,3,5,0,3,2.6218495656463876,0.917324005589486,0.
+0,1,10,3,3,3,7,3,4,1.2601693603785153,2.8311848091672465,0.
+5,3,3,1,0,2,0,1,2,1.5322630072331016,2.4612811983034675,0.
+2,0,10,3,4,3,8,2,3,2.581504744416619,2.3352397148922988,0.
+8,2,10,1,8,2,10,1,1,1.1752269484482167,1.7646051998689414,0.
+0,2,4,1,10,2,7,0,2,1.1544888515337508,2.637958865288499,0.
+4,0,6,3,6,1,7,3,3,1.4304562327889254,1.498204529248202,0.
+8,2,6,1,3,3,6,0,3,2.3823736912605336,2.0477368833060705,0.
+4,2,2,3,10,1,3,1,1,2.5205545321391254,0.7220948301952252,0.
+10,0,7,1,8,3,9,2,1,1.2260538418022682,2.451257871532359,0.
+0,0,0,3,10,1,10,3,3,2.006553045177734,1.173388489802477,0.
+6,3,3,3,1,3,1,2,5,2.9612859268853233,1.4353422269551883,0.
+2,1,2,0,2,3,1,2,1,1.3130842056605396,1.1311416250752515,0.
+3,3,9,3,5,1,3,3,3,0.4024708536057622,2.0642879809097643,0.
+8,2,7,2,3,2,7,2,0,1.2953268563208105,0.5693758931801445,0.
+0,3,9,1,8,1,8,1,2,2.152069679767771,1.41039953394532,0.
+10,3,10,1,2,1,2,3,3,1.922953978768665,1.3543149100065444,0.
+4,2,5,1,10,3,8,0,3,0.6098612038176467,0.5440300294800884,0.
+8,2,7,1,3,1,8,2,2,0.608144934570964,2.548777086726803,0.
+9,2,3,0,6,0,10,2,2,1.7250046113327357,1.415947697499866,0.
+4,1,3,2,7,0,3,3,3,1.5376308517695847,2.8988547754211718,0.
+5,3,4,0,6,3,4,2,3,0.7636088061341804,0.7149031039865208,0.
+4,3,9,3,2,0,6,2,2,1.8376914709492826,1.479579605670744,0.
+2,3,10,2,5,2,10,0,2,0.8658425814507988,1.1218947490443987,0.
+6,0,9,1,2,3,10,2,1,2.7269511334702985,1.2873996903274714,0.
+5,3,1,1,4,1,2,2,3,2.0472755894627443,1.4404427131274424,0.
+7,3,5,1,5,3,1,0,3,2.1329904508839563,1.7676149035568383,0.
+0,1,6,1,3,2,1,3,1,1.1352856084425453,2.0153824999161527,0.
+9,1,4,2,9,2,1,0,2,0.5735267931557,0.883812410190215,0.
+10,2,8,0,2,2,0,0,2,2.286451077255511,2.3140961306516,0.
+3,3,3,3,2,2,0,2,2,0.7388292443690174,1.6936545691041545,0.
+10,2,2,1,10,0,10,2,2,1.4443453550869583,1.112093007841184,0.
+7,3,6,1,10,3,2,2,4,1.8267825436432625,1.0428273092046476,0.
+0,0,1,3,9,2,9,1,3,2.4360668218156167,1.0158379194857559,0.
+9,0,6,3,8,1,5,2,3,2.4459546670789223,2.5273734233477834,0.
+8,2,0,2,9,3,2,2,1,0.9849892765217048,0.31610510276469483,0.
+2,0,6,0,3,1,9,1,0,1.646375524767966,0.39139222602227974,0.
+2,3,0,2,6,0,10,1,1,0.3040006591059585,1.784871103068122,0.
+7,2,2,2,4,3,6,0,3,1.5667275423919191,1.7334318105948183,0.
+5,3,2,1,1,2,4,0,2,0.4724440845832545,2.0191544118568805,0.
+2,1,3,2,6,1,8,3,2,1.6569706450379993,1.22551530354268,0.
+4,0,8,2,10,3,6,2,2,1.5517249658736345,1.763047858243108,0.
+6,2,4,2,0,3,3,0,3,2.056076795383781,1.693356584470254,0.
+1,3,6,1,4,1,3,3,3,1.9168275535161614,2.2496624307522994,0.
+6,1,1,1,8,1,3,0,1,0.7505340273835843,1.1586206830861343,0.
+1,2,9,2,8,2,1,0,2,2.5465024812722747,0.40185310698150234,0.
+6,2,1,3,8,3,5,3,2,1.0806918086518902,2.785171921529357,0.
+2,3,9,1,0,1,7,3,2,2.649104146816554,0.22330923947941184,0.
+3,3,2,2,8,2,7,1,1,1.4153027357986057,0.7592362643035888,0.
+1,2,4,1,2,0,8,1,1,0.3716613874614527,2.741327801112848,0.
+8,2,4,0,7,3,8,2,2,2.772853450624477,0.7349370680695975,0.
+4,2,10,1,9,1,10,2,2,2.3205630100773416,0.7942049271730163,0.
+1,3,6,3,6,1,4,3,2,1.3281998811443967,1.0363866114967322,0.
+3,2,7,3,3,1,3,3,4,0.8164145918403629,1.064457528831932,0.
+5,1,6,0,6,1,5,0,1,2.7770950000390737,1.3710296853982933,0.
+4,3,3,3,3,1,2,2,3,2.714056705000763,2.4602516368227603,0.
+4,1,5,2,3,2,7,3,1,0.35146521178282075,0.3571498607132635,0.
+9,1,8,2,9,1,6,1,1,1.2274240317357688,1.6715623227186507,0.
+1,1,3,3,8,3,7,0,3,1.9317099576956442,1.552216114675943,0.
+9,0,3,1,8,0,2,1,1,2.9043448611385188,2.7338758740571345,0.
+6,0,8,2,2,2,4,1,2,1.675816848063465,0.8679109639840443,0.
+9,1,9,1,0,1,1,0,1,0.596382672070606,1.7119233696062492,0.
+7,2,3,1,5,1,0,3,3,1.374644910756329,2.8427817664580024,0.
+6,3,10,3,7,1,10,1,2,2.6167816000584327,0.8461412388948104,0.
+5,2,7,1,6,3,2,3,2,1.372612635683653,2.0961606258142194,0.
+8,1,7,2,8,3,1,2,3,0.9231592095990986,0.4730174634532238,0.
+3,1,1,1,10,1,10,3,2,0.3461817829038991,1.5964106156876117,0.
+8,0,7,3,0,2,4,1,3,2.096022436900264,1.9588093081620253,0.
+4,3,9,3,5,1,9,1,0,2.449961912097761,0.4920648568697943,0.
+4,0,8,1,6,2,1,3,1,2.7092656850999868,1.2113036872224647,0.
+2,2,4,3,2,2,1,0,2,0.7822138935161549,0.8243720955516407,0.
+5,3,9,3,8,0,7,3,3,1.245926765783301,2.995825714452799,0.
+9,1,9,0,5,3,1,3,1,0.9447956603373777,1.9688135363735153,0.
+5,2,2,2,3,3,7,1,2,2.2346399253675413,0.7417391778003699,0.
+0,1,8,3,0,2,2,3,2,1.0232322455743281,1.8925470263446593,0.
+5,3,2,3,8,2,9,0,2,2.9817835096507697,2.111249932049714,0.
+9,2,8,2,1,1,4,3,4,0.4083625562135289,0.40305961637541055,0.
+4,3,9,2,8,2,10,1,3,2.7964161333721353,2.78601941028056,0.
+1,1,10,0,10,0,10,1,1,2.8918303445485822,2.7282191833543576,0.
+3,1,1,3,0,1,9,2,3,0.5883986968362152,1.1131593470756491,0.
+3,1,5,3,6,3,10,0,3,1.29122181890919,2.1687719665184835,0.
+5,1,1,0,10,1,9,1,1,2.7757481174210463,2.3835512920725517,0.
+3,3,1,3,10,0,6,1,1,1.127749467977892,1.3592596107053452,0.
+0,1,10,2,10,2,4,0,2,1.403190963591391,1.7426120757596841,0.
+8,2,0,3,2,3,4,2,3,2.21376917027845,1.2373565371204496,0.
+6,3,2,2,0,2,6,2,4,2.0909192232076848,2.9595518129553477,0.
+2,0,4,3,0,3,2,0,3,2.9389053651691066,0.2571903449038695,0.
+3,0,8,2,2,2,8,3,2,1.1292930350844212,2.662958886705143,0.
+4,2,9,3,3,1,8,2,2,1.2069732981501304,0.6603448818737543,0.
+8,1,1,2,9,2,2,3,2,1.176331991745739,2.9348246139214273,0.
+8,2,9,1,3,2,6,0,2,1.2579465199621849,2.8708297284639936,0.
+2,1,0,0,4,2,0,2,1,0.5068840995027317,1.4926314796061328,0.
+8,0,5,2,0,1,6,1,2,2.083459207268281,2.044350098956075,0.
+8,1,3,1,8,3,3,1,2,0.27182349339562206,0.25278077977555524,0.
+9,0,9,3,8,2,7,1,3,2.0842851912303697,1.4617238592678117,0.
+6,0,7,3,4,2,5,3,3,0.5438422711246416,1.1960642527951837,0.
+3,0,0,3,1,3,5,2,3,2.142935683266975,0.3991912324294198,0.
+5,2,1,1,6,3,6,2,2,1.7422616583207757,1.1999207722763399,0.
+1,2,1,0,8,2,9,0,2,0.32666758600821844,2.6968272363747605,0.
+3,2,0,0,10,1,6,1,2,1.6221941314416704,1.1434410020543084,0.
+0,3,1,2,10,1,3,1,2,0.9726067645442806,0.5905895434956969,0.
+0,1,0,0,1,0,7,1,1,2.900837043628563,2.949017267115262,0.
+8,3,8,1,2,0,0,2,2,1.5234948494086824,0.3368477864264441,0.
+7,2,4,2,8,1,6,1,1,0.650243384315508,0.5631520583276046,0.
+3,3,9,3,2,0,9,0,0,2.3807738065036146,0.579866727058429,0.
+4,3,7,0,10,2,6,3,3,1.3173017138545768,0.5962764653299519,0.
+8,3,0,3,0,3,2,1,4,0.6053424778386112,0.8254011764061571,0.
+9,1,3,1,1,0,2,2,2,0.9352316673114225,2.648521223427407,0.
+4,3,6,2,1,3,8,3,3,1.1811263754353778,0.27949826874059935,0.
+6,2,9,1,0,3,8,3,2,1.7891686343587905,2.5475273097794897,0.
+3,3,6,3,1,1,4,2,1,0.8652651806579512,2.144133642324886,0.
+8,2,5,1,10,3,3,1,3,1.3153722390010687,2.370685736109543,0.
+9,1,5,2,7,3,2,1,2,1.4002015422875607,1.5704174280943155,0.
+4,3,8,3,9,1,4,2,1,1.1067444816569871,2.7633308770019767,0.
+8,1,5,1,7,0,7,1,1,1.3033564335609653,1.8205597872081913,0.
+6,2,10,2,0,2,9,0,2,1.417347861506408,0.8228017043867935,0.
+3,3,5,2,8,1,5,1,1,2.3178061623921504,2.446585257893128,0.
+0,2,5,2,5,1,0,0,1,0.417745684372842,0.962045539464802,0.
+8,1,8,1,8,1,5,0,1,2.241624712766236,2.117102285343411,0.
+5,2,1,1,2,0,2,3,3,1.452790633495452,2.3926515680383504,0.
+7,1,0,3,9,1,9,1,2,1.320602559814661,0.5022153836743799,0.
+10,2,8,3,1,1,4,3,4,0.8992253656072662,2.3674150876262807,0.
+2,3,9,0,10,1,9,2,3,1.0158467671907352,1.376506498457696,0.
+1,3,1,1,4,1,4,1,2,1.8583872135299577,2.4330946543470144,0.
+8,0,9,1,3,0,7,1,1,2.9549925009110574,0.7408650898965572,0.
+5,1,3,1,2,2,0,0,2,0.7685503537061233,1.8557635823419316,0.
+10,1,7,0,3,1,7,1,1,2.91588044339084,0.4955239475168276,0.
+6,2,10,0,10,3,1,1,2,1.2165093404056426,0.48447018353308113,0.
+8,3,0,2,9,0,9,3,3,1.2458573416728456,1.8984850604019927,0.
+5,1,6,1,5,2,2,3,2,1.6678698181298168,0.27271548312207194,0.
+1,0,5,1,6,3,7,3,1,2.4353366447686176,2.251622116348118,0.
+4,1,6,2,10,3,8,1,3,0.9785359746307756,2.2569876645542513,0.
+7,1,8,0,3,0,8,1,1,2.6254607678979696,1.4870987629995316,0.
+7,2,4,0,4,1,6,2,2,2.4860925929709445,2.078702116156186,0.
+4,3,10,0,2,3,4,1,3,2.7833993206816183,2.4091684481860574,0.
+10,2,3,2,10,0,3,2,2,2.600869455586518,2.7023313585291886,0.
+10,0,3,2,1,1,5,1,2,0.6308140012941346,2.192072324433969,0.
+0,2,1,3,8,3,9,2,5,0.47211349059914465,1.2835269519982067,0.
+1,3,5,3,9,0,8,1,1,0.39051019613873894,0.6191713045924181,0.
+8,2,1,2,5,2,10,1,1,2.8083327076207043,2.7889141797722985,0.
+0,1,6,2,4,3,1,3,3,2.576880910721891,2.8940089415951142,0.
+10,2,9,2,3,1,7,1,2,0.4344721828159619,2.6268647157878835,0.
+4,1,0,3,8,0,10,3,3,2.3791812481389023,0.41047900744269716,0.
+9,0,0,3,9,2,2,1,3,0.3591915419041567,2.665840985111095,0.
+9,3,6,2,5,1,5,2,3,2.1661396893720006,1.1610985541712369,0.
+7,3,2,3,2,1,5,2,2,1.2178868457103484,1.906937613481415,0.
+9,1,4,2,1,2,5,2,1,1.76351959987508,0.310172698567285,0.
+10,1,4,2,5,1,5,1,2,2.9057790676388597,2.246935216027948,0.
+9,3,6,1,7,0,5,2,2,1.8228283848861624,0.5136268936144317,0.
+4,1,2,1,10,2,1,2,2,2.4146521536083805,0.4173254666656194,0.
+9,1,2,1,3,2,8,0,2,2.091075002425989,1.6190364383142475,0.
+5,2,7,1,2,2,7,1,2,1.9881551763498084,0.9065836043584428,0.
+5,3,6,2,4,1,7,2,1,0.2437619648775846,0.4633502649952983,0.
+6,1,8,1,10,2,5,1,1,1.979325271197399,1.1385489265368367,0.
+3,2,6,2,4,0,8,3,3,2.761222760478871,2.281187219578926,0.
+7,3,2,2,3,0,4,1,1,1.387675975897055,2.978542623739007,0.
+0,2,8,0,0,1,2,2,2,0.7224763006498796,1.564980597384034,0.
+3,3,7,1,10,2,10,2,3,2.9525286612356716,0.8408929108171121,0.
+0,1,7,2,2,1,9,0,1,0.9126033819378172,2.8275866134399896,0.
+9,2,1,0,8,0,0,2,2,1.0314013900322792,0.6732978424388585,0.
+8,2,9,1,7,1,9,0,1,1.3865069980063334,2.540865404640523,0.
+6,2,0,2,2,1,9,2,3,0.4235665273207312,1.6846937203861359,0.
+5,2,5,2,5,0,3,1,1,2.7659439645673256,1.0115836172242267,0.
+0,1,9,3,6,3,2,1,3,1.2065541331332832,1.0405234817364937,0.
+8,1,6,1,4,0,1,1,1,0.7227473981225692,1.155456043230863,0.
+2,3,8,0,1,0,3,3,3,1.9299731710816523,0.6103164598521023,0.
+4,0,2,3,0,2,1,2,3,2.8621863382346,2.9310316681902524,0.
+8,3,5,1,8,2,6,3,4,0.45896263960829176,1.9911804294619255,0.
+1,0,3,2,6,2,4,0,2,2.335589089980126,1.0652452461608548,0.
+3,1,1,3,1,2,10,3,3,1.0154942249838212,1.8603629108124498,0.
+2,1,6,0,3,2,2,2,1,0.6444259802836747,1.4787144338375233,0.
+7,1,10,2,2,1,3,2,3,0.7147755411709733,0.741536139217259,0.
+0,2,1,3,2,1,8,1,1,1.8857744159033487,2.4763969783476716,0.
+2,3,8,0,10,3,4,0,3,1.9960210966818677,2.9863811890951855,0.
+2,3,2,3,2,1,5,2,2,1.4578550110452344,2.0655616115569337,0.
+8,3,7,2,0,0,0,2,2,2.1690558898141488,1.6547046613975827,0.
+7,2,3,1,6,3,3,0,3,2.553795269242924,2.98657012675977,-5.8584894789029045
+5,2,2,2,8,0,7,0,0,2.111885434964618,1.9937513659932842,0.
+8,3,8,0,8,2,0,3,3,2.8705609121506805,1.3331637714091094,0.
+3,2,6,3,1,2,0,1,1,2.34251648574082,0.5435767229545014,0.
+4,0,6,2,6,0,9,2,2,1.20338522636726,2.788327210207373,0.
+7,1,1,0,10,1,0,0,1,2.3923670262306302,2.4029578275837586,0.
+0,1,8,3,6,2,0,3,2,1.203543583831212,2.036143250730075,0.
+5,3,3,1,10,3,9,1,2,2.7314904241739955,0.21670230017016534,0.
+7,3,0,3,4,0,9,3,3,0.7221556673410454,2.504574860198268,0.
+0,1,9,0,7,2,2,1,1,0.6129040077484231,0.9377073179769542,0.
+8,3,8,3,8,0,6,1,1,0.2916837472356941,1.6706587000955677,0.
+5,0,8,1,2,1,9,1,1,1.1071388188914164,2.907305251123775,0.
+4,2,3,3,2,1,9,0,1,2.6049414020919137,1.1423144552963707,0.
+0,2,3,2,6,0,3,1,1,0.6589475462254843,0.3458167300233934,0.
+5,1,7,2,1,1,9,1,2,0.37152940436379245,2.257037185691721,0.
+5,2,5,2,9,3,3,2,4,1.1411180980151068,1.360756047176233,0.
+7,2,5,3,4,3,6,3,5,2.885010245830152,0.6362691751312788,0.
+5,2,9,3,10,1,0,1,2,0.37940369466801593,1.2334735865845734,0.
+0,0,7,2,1,2,7,1,2,1.5479759698453543,1.9801733480094859,0.
+2,2,4,0,10,0,10,2,2,1.3293709201965624,1.2689927819499371,0.
+9,0,3,1,0,2,8,2,1,2.164068098414817,1.2453603926890962,0.
+0,2,9,2,3,1,4,0,1,0.6696321137903438,2.6084879932864853,0.
+8,1,3,3,4,1,6,1,2,2.042466846990008,1.8892178256161216,0.
+2,3,2,1,10,3,10,2,3,1.4323483837153708,1.0105488574564596,0.
+5,0,10,3,7,3,10,2,3,2.726039671844071,2.8549137956975468,0.
+2,3,1,2,3,3,10,1,3,2.351151357868736,0.36675214929021926,0.
+2,0,7,0,9,1,3,1,0,0.7394282604433373,0.2998036281738994,0.
+4,3,0,3,6,2,9,2,4,1.4189068677052719,2.9473449090470467,0.
+9,2,1,3,9,2,7,0,2,2.2125653578966116,0.6449925605947713,0.
+3,3,9,3,7,0,2,0,0,2.6527765476432785,0.7017535912522956,0.
+3,2,1,3,10,3,10,2,3,1.3545764993952725,1.4463176767351857,0.
+6,0,6,3,7,3,10,1,3,0.5231215491635299,1.5951002259931988,0.
+10,1,0,3,0,0,3,2,2,0.932146690038874,0.6396108542280041,0.
+1,3,4,3,3,1,5,1,2,2.6568428894223652,1.0041895351349086,0.
+2,0,6,3,7,2,5,2,3,2.1388490384614993,1.084038962927921,0.
+1,0,4,3,3,3,6,1,3,1.1954487430499983,1.3480169360596266,0.
+4,2,10,3,7,0,8,3,3,0.9989698852671967,0.31492500891608577,0.
+0,1,5,1,4,1,7,1,2,0.2915533255281222,0.8488169939496006,0.
+5,3,8,2,6,3,10,0,3,1.3273506939181288,1.9576844972216119,0.
+4,1,1,1,8,1,3,1,0,2.5512768516580167,1.1776908814863933,0.
+0,1,7,2,4,3,5,2,3,2.55835426097922,0.7999378464474947,0.
+7,1,3,2,4,3,4,3,3,2.9469964384408645,1.0130439824706348,0.
+10,1,10,1,4,0,9,2,2,2.29282320065841,1.5895923699679324,0.
+3,3,8,2,9,2,8,3,5,1.4288368272508416,0.32914872565251896,0.
+5,2,10,0,7,1,4,2,2,2.0963238426388164,2.5184909838350205,0.
+7,3,4,0,1,2,10,3,3,0.7462285283562706,1.9678766341114837,0.
+5,1,6,0,6,1,1,0,1,2.742103144450633,1.0100914398280834,0.
+6,2,8,2,1,0,8,0,0,2.9686558639239102,2.571359181752765,0.
+6,2,1,0,8,3,0,2,2,2.8725089493368277,2.890444927773988,0.
+3,2,6,3,9,1,2,2,1,0.9140224473432359,1.7614216539656447,0.
+9,1,3,1,2,3,8,2,2,1.5224977087087135,1.7130595192975742,0.
+2,1,5,1,1,1,2,1,0,0.620215788594543,2.125090314126055,0.
+3,3,2,2,9,3,10,3,5,2.1571418399260214,0.2352162581649102,0.
+4,1,3,3,8,2,1,0,2,2.5710875717760064,1.8947511401060053,0.
+8,2,5,2,6,1,8,1,0,1.5775238404669114,2.9176906027237717,0.
+10,2,1,2,2,3,3,3,2,2.9869659703459304,0.6494176627258623,0.
+9,1,7,1,7,0,5,0,0,0.35843660067087235,1.1434109912448518,0.
+8,0,0,3,0,1,7,3,3,2.5846819654664843,0.30192303129821685,0.
+9,0,10,1,3,0,4,1,1,2.2400726780346583,2.76156913953512,0.
+0,2,6,1,3,1,3,0,1,1.2200817282296477,0.8137895784475973,0.
+4,3,0,1,10,1,3,2,3,1.1171419500565256,2.8469536184972846,0.
+6,1,3,2,2,3,3,0,3,0.8843975271072293,2.0058645383535723,0.
+0,2,4,3,3,1,1,2,1,2.612997115865329,0.34104461891724425,0.
+3,0,1,3,6,3,10,2,3,2.716076032829566,1.791883854821867,0.
+0,3,1,3,1,2,7,2,2,2.989190167535284,0.6804444348264811,0.
+7,0,7,0,6,1,8,1,0,1.4624896265248113,1.7583815597493881,0.
+6,3,6,3,1,1,9,2,3,2.4493465867634194,1.3172740005089132,0.
+3,2,1,2,5,3,8,1,2,2.9637052407119135,1.386678142676876,0.

test/racahs_reduction_formula.jl

@@ -0,0 +1,43 @@
+using QuadGK
+include("../math.jl")
+
+# Testing ∫dp p u' u where u' and u' may have different l
+
+function explicit_integral(np, lp, n, l, μω)
+    integrand(p) = p * (-1)^(n + np) * 1/sqrt(μω) * ho_basis(l, n, 1/sqrt(μω) * p) * ho_basis(lp, np, 1/sqrt(μω) * p)
+    (integral, _) = quadgk(integrand, 0, Inf; atol=1e-10)
+    return integral
+end
+
+n_list = 0 : 10
+l_list = 0 : 5
+μω_list = [1, 2.66]
+
+for np in n_list
+    for n in n_list
+        for lp in l_list
+            for l in l_list
+                for μω in μω_list
+                    println("Checking n=$n, np=$np, l=$l, lp=$lp, μω=$μω")
+                    expl = explicit_integral(np, lp, n, l, μω)
+                    anal = integral(np, lp, n, l, μω)
+                    @assert isapprox(expl, anal; atol=1e-8) "Explicit=$(expl) and analytical=$(anal)"
+                end
+            end
+        end
+    end
+end
+
+# Testing the full Racah's reduction formula
+
+test_data = DelimitedFiles.readdlm("test/racah_test_data.csv", ',') # format = n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2, racah_matrix_element
+test_data_ints = Int.(test_data[:, 1:(end - 3)])
+test_data_floats = test_data[:, (end - 2):end]
+
+for i in 1 : size(test_data)[1]
+    n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ = test_data_ints[i, :]
+    μ1ω1, μ2ω2, racah_matrix_element = test_data_floats[i, :]
+    println("Checking n1'=$n1p, l1'=$l1p, n2'=$n2p, l2'=$l2p, n1=$n1, l1=$l1, n2=$n2, l2=$l2, Λ=$Λ, μ1ω1=$μ1ω1, μ2ω2=$μ2ω2")
+    val = racahs_reduction_formula(n1p, l1p, n2p, l2p, n1, l1, n2, l2, Λ, μ1ω1, μ2ω2)
+    @assert isapprox(val, racah_matrix_element; atol=1e-8) "Formula=$(val) and test reference=$racah_matrix_element"
+end