# Algorith integration rules
#
# References:
#
#  - J. Berntsen, T.O. Espelid and A. Genz. "An adaptive algorithm for the
#    approximate calculation of multiple integrals". ACM Trans. Math. Soft.,
#    17, 437-451, 1991.
#
#  - J. Berntsen and T.O. Espelid. "On the construction of higher degree
#    three-dimensional embedded integrations rules". SIAM J. Numer. Anal., 25,
#    222-234, 1988.
#
#  - A. C. Genz and A. A. Malik. "An imbedded family of fully symmetric
#    numerical integration rules". SIAM J. Numer. Anal., 20, 580-588, 1983.


function points_per_generator(n)
    p₁ = 2n
    p₂ = p₁ * (n - 1)
    p₃ = 2p₂
    p₅ = p₃ * (n - 2)
    p₄ = p₅ ÷ 3
    p₆ = 1 << n
    return (1, p₁, p₂, p₃, p₄, p₅, p₆)
end


### 11th degree rule
#
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 5, 2, 0, 3, 2, 0); (n = 3)

# function evaluations per region
# sum(map(*, (K₀, K₁, K₃, K₄, K₅, K₆), points_per_generator(3)))
const L₁₁ = 127

const α₁ = 0.19
const α₂ = 0.50
const α₃ = 0.75
const α₄ = 0.80
const α₅ = 0.995

const β₁ = 0.9987344998351400
const β₂ = 0.7793703685672423

const ɛ₁ = 0.9999698993088767
const ɛ₂ = 0.7902637224771788
const ɛ₃ = 0.4403396687650737

const ζ₁ = 0.4378478459006862
const η₁ = 0.9549373822794593

const ζ₂ = 0.9661093133630747
const η₂ = 0.4577105877763134

const α = (α₁, α₂, α₃, α₄, α₅)
const β = (β₁, β₂)
const ɛ = (ɛ₁, ɛ₂, ɛ₃)
const ζ = (ζ₁, ζ₂)
const η = (η₁, η₂)

# weight for (0, 0, ..., 0)
const p₀ = (  0.007923078151097301  ,  1.7150062482246855    ,
              1.9360149789495251    ,  0.517082819560578     ,
              2.0544045038185197    )
# weights for α
const a₁ = (  0.06797177392788224   , -0.3755893815889209    ,
             -0.3673449403754268    ,  0.01445269144914044   ,
              0.013777599884901202  )
const a₂ = (  0.0010869865388070375 ,  0.1488632145140549    ,
              0.029297786578981855  , -0.3601489663995932    ,
             -0.576806291790441     )
const a₃ = (  0.18386336622127977   , -0.2497046640620823    ,
             -0.1151883520260315    ,  0.3628307003418485    ,
              0.03726835047700328   )
const a₄ = (  0.0336211977782927    ,  0.1792501419135204    ,
              0.050866582208722165  ,  0.007148802650872727  ,
              0.0068148789397772195 )
const a₅ = (  0.010137511233340512  ,  0.0034461267589738897 ,
              0.04453911087786472   , -0.09222852896022965   ,
              0.0572316973385185    )
# weights for β
const b₁ = (  0.0016876486839852163 , -0.005140483185555828  ,
             -0.02287828257125902   ,  0.01719339732471725   ,
             -0.044930187438112855  )
const b₂ = (  0.13464685645128072   ,  0.006536017839876424  ,
              0.029089262163458317  , -0.102141653746035     ,
              0.027292365738663484  )
# weights for ɛ
const d₁ = (  0.0017501458846003852 , -0.00065134549392297   ,
             -0.002898884350669207  , -0.007504397861080488  ,
              0.0003547473950556992 )
const d₂ = (  0.07752336383837458   , -0.006304672433547204  ,
             -0.028059634133074954  ,  0.01648362537726711   ,
              0.01571366799739551   )
const d₃ = (  0.2461864902770248    ,  0.012669593997882627  ,
              0.05638741361145884   ,  0.05234610158469334   ,
              0.04990099219278569   )
# weights for ζ and η
const e₁ = (  0.06797944868483038   , -0.005454241018647931  ,
             -0.02427469611942451   ,  0.014454323316130661  ,
              0.013779155526667704  )
const e₂ = (  0.014199628233007138  ,  0.004826995274768427  ,
              0.021483070341828822  ,  0.003019236275367777  ,
              0.0028782064230998723 )

const wa = (a₁, a₂, a₃, a₄, a₅)
const wb = (b₁, b₂)
const wd = (d₁, d₂, d₃)
const we = (e₁, e₂)

# unique(-(w[i+1] / w[i]) for w in (p₀, wa..., wb..., wd..., we...), i = 2:4)
const μ₁ = ( -1.1288675950618956    ,  -0.9780493229637744   ,
             -0.19681011641875891   ,  -0.46129836003941455  ,
             -0.28377429253731273   , -12.924397154539543    ,
             -4.45060935819115      )
const μ₂ = ( -0.2670861667822143    ,   0.0393436518667435   ,
             12.29270222952484      ,   3.1498905397991295   ,
             -0.14054025925191632   ,   2.0707312548995196   ,
              0.7515160839177919    ,   3.5113181342338917   ,
             -2.5887192979423617    ,   0.5874497614292568   ,
             -0.9283295372507698    ,   0.5954481673022622   )
const μ₃ = ( -3.9730666463921054    ,  -0.9532895608672675   ,
             -1.6015769739859858    ,  -0.10271553769262119  ,
              0.6205422333386417    ,   2.613223354846873    ,
              0.26720113428477693   ,   0.047271933288012324 )

# norms for each μ in μᵢ
#
# let ws = ((p₀,), wa, wb, wd, we)
#     ns = points_per_generator(3)[[1, 2, 3, 5, 6]]
#     μs = (μ₁, μ₂, μ₃)
#     s = Float64[]
#     for i in 1:length(μs), μ in μs[i]
#         sum = 0.0
#         for j in 1:length(n), w in ws[j]
#             sum += n[j] * abs(μ[i] * w[i+1] + w[i+2])
#         end
#         push!(s, 1/sum)
#     end
#     (s...)
# end
const s₁ = (  0.19485823831443944   ,   0.21677927747604941  ,
              0.155629450215002     ,   0.19365063964604642  ,
              0.16991906305920448   ,   0.010482710226773457 ,
              0.035771397525237954  )
const s₂ = (  0.11602998413860789   ,   0.12595256004163577  ,
              0.010803324347691502  ,   0.043690435761900195 ,
              0.12171190392783988   ,   0.0592433918191275   ,
              0.09766409013104478   ,   0.03936956866824965  ,
              0.037753982678781134  ,   0.10512211588403095  ,
              0.07529459306375538   ,   0.10482954809408154  )
const s₃ = (  0.038152895079104944  ,   0.1253631725603185   ,
              0.10318930777035908   ,   0.13160544230723548  ,
              0.08802194500708062   ,   0.03776436580400292  ,
              0.11024221575948996   ,   0.1221778013467861   )

# heuristic constants
const c = (4, 4, 0.5, 3, 0.5, 0.25)

### 9th degree rule
#
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 4, 2, 1, 0, 0, 0); (n = 2)
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 4, 1, 1, 1, 0, 1); (n ≥ 3)

# δ² = big"0.4707"
#const δ  = parse(BigFloat, "0.686075797561756291400285229032935720681" *
#                             "3963449254746025939952564970604511034176")

# λ₁² = 4 / (15 - 5 / δ²)
#const λ₁ = parse(BigFloat, "0.955907315804539012385720839434366012946" *
#                             "392550072657165493502911217854564875694")

# s = p = λ₁²
# d = (1 - λ₁² / δ²) / 27
# λ₂² = (1 // 7 - s / 5 - d) / (1 // 5 - s / 3 - d / δ²)
#const λ₂ = parse(BigFloat, "0.406057174738239735599606880474655826651" *
#                             "5592938007285994340533769219899878305832")

# s += λ₂²
# p *= λ₂²
# d *= (1 - λ₂² / δ²) / 3
# λ₃² = (1 // 9 - s / 7 + p / 5 - d) / (1 // 7 - s / 5 + p / 3 - d / δ²)
#const λ₃ = parse(BigFloat, "0.895254709252356257641540961306199707857" *
#                             "0302757591422795174722770397905378692193")
