
using FastGaussQuadrature
using Polynomials

using DecFP: Dec128


function getCoefficientsGLRK1(T=Float64)
    a = reshape(Array([1//2]), (1,1))
    b = Vector([1//1])
    c = Vector([1//2])
    o = 2

    CoefficientsRK(T, :GLRK1, o, a, b, c)
end

function getCoefficientsGLRK2(T=Float64)
    a = @dec128 [
          [1/4       1/4-√3/6]
          [1/4+√3/6  1/4     ]
        ]
    b = @dec128 [1/2,      1/2     ]
    c = @dec128 [1/2-√3/6, 1/2+√3/6]
    o = 4

    CoefficientsRK(T, :GLRK2, o, a, b, c)
end

function getCoefficientsGLRK3(T=Float64)
    a = @dec128 [
          [5/36         2/9-√15/15  5/36-√15/30]
          [5/36+√15/24  2/9         5/36-√15/24]
          [5/36+√15/30  2/9+√15/15  5/36       ]
        ]
    b = @dec128 [5/18,        4/9,        5/18       ]
    c = @dec128 [1/2-√15/10,  1/2,        1/2+√15/10 ]
    o = 6

    CoefficientsRK(T, :GLRK3, o, a, b, c)
end

function getCoefficientsGLRK4(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
         [(-√30/144 + 1/8)
          (-13*√42/336 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/28 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/14 - √(60*√30 + 450)/144 - √42/168 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√(60*√30 + 450)/144 - √3/12 + √10/16 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (√30/144 + 1/8)
          (-√35/14 + √42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-13*√42/336 - √35/28 - √(60*√30 + 450)/144 - √3/12 + √10/16 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√10/16 - √(60*√30 + 450)/144 + √3/12 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (-√42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8 + √35/14)/√(-2*√30 + 15)
          (√30/144 + 1/8)
          (-13*√42/336 - √35/28 - √10/16 - √(60*√30 + 450)/144 + √3/12 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√(60*√30 + 450)/144 + √42/168 + √35/14 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (-13*√42/336 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/28 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√30/144 + 1/8)]
        ]), (4,4))'

    b = @dec128 [
          49/(216+12*√30),
          49/(216-12*√30),
          49/(216-12*√30),
          49/(216+12*√30)
        ]

    c = @dec128 [
          1/2-√(3/7+2*√30/35)/2,
          1/2-√(3/7-2*√30/35)/2,
          1/2+√(3/7-2*√30/35)/2,
          1/2+√(3/7+2*√30/35)/2
        ]

    o = 8

    CoefficientsRK(T, :GLRK4, o, a, b, c)
end

function getCoefficientsGLRK5(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
         [-13*√70/3600 + 161/1800
          -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 161/1800
          -92*√(14*√70 + 245)/8505 + 4*√(20*√70 + 350)/1215 + 32/225
          -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 161/1800
          -2*√(14*√70 + 245)/315 - 13*√70/3600 + √(20*√70 + 350)/360 + 161/1800];
         [-113*√(-14*√70 + 245)/34020 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 11*√(20*√70 + 350)/6480 + 59*√(-20*√70 + 350)/19440 + 161/1800
          13*√70/3600 + 161/1800
          -92*√(-14*√70 + 245)/8505 - 4*√(-20*√70 + 350)/1215 + 32/225
          -2*√(-14*√70 + 245)/315 - √(-20*√70 + 350)/360 + 13*√70/3600 + 161/1800
          -11*√(20*√70 + 350)/6480 - 113*√(-14*√70 + 245)/34020 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 59*√(-20*√70 + 350)/19440 + 161/1800];
         [-23*√(20*√70 + 350)/5760 - 13*√70/3600 + 161/1800 + 11*√(14*√70 + 245)/1440
          13*√70/3600 + 23*√(-20*√70 + 350)/5760 + 11*√(-14*√70 + 245)/1440 + 161/1800
          32/225
          -11*√(-14*√70 + 245)/1440 - 23*√(-20*√70 + 350)/5760 + 13*√70/3600 + 161/1800
          -11*√(14*√70 + 245)/1440 - 13*√70/3600 + 161/1800 + 23*√(20*√70 + 350)/5760];
         [-59*√(-20*√70 + 350)/19440 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 11*√(20*√70 + 350)/6480 + 161/1800
          13*√70/3600 + √(-20*√70 + 350)/360 + 2*√(-14*√70 + 245)/315 + 161/1800
          4*√(-20*√70 + 350)/1215 + 92*√(-14*√70 + 245)/8505 + 32/225
          13*√70/3600 + 161/1800
          -59*√(-20*√70 + 350)/19440 - 11*√(20*√70 + 350)/6480 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 161/1800];
         [-√(20*√70 + 350)/360 - 13*√70/3600 + 161/1800 + 2*√(14*√70 + 245)/315
          -11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
          -4*√(20*√70 + 350)/1215 + 32/225 + 92*√(14*√70 + 245)/8505
          -4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
          -13*√70/3600 + 161/1800]
        ]), (5,5))'

    b = @dec128 [
          1/(4600/729+1300*√70/5103),
          1/(4600/729-1300*√70/5103),
          64/225,
          1/(4600/729-1300*√70/5103),
          1/(4600/729+1300*√70/5103)
        ]

    c = @dec128 [
          1/2-√(5/9+2*√70/63)/2,
          1/2-√(5/9-2*√70/63)/2,
          1/2,
          1/2+√(5/9-2*√70/63)/2,
          1/2+√(5/9+2*√70/63)/2
        ]

    o = 10

    CoefficientsRK(T, :GLRK5, o, a, b, c)
end

function getCoefficientsGLRK6(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
            ["+0.0428311230947925862600740355431832233817056253710109955996588599497364402",
             "-0.0147637259971974124753725910605206514420001878003408840631515394063317448",
             "+0.00932505070647775119143888450800314858828816322281640330116548464376180630",
             "-0.00566885804948351190092125641621650656214351346713436514621339322133564885",
             "+0.00285443331509933513092928583011602153367099528247654439080898011903375405",
             "-0.000812780171264762112299135651562540066903951464973606518543200348171357290"];
            ["+0.0926734914303788631865122917633203161433485534502333679517368248345204102",
             "+0.0901903932620346518924583784594290279153804731866863705724348100592850095",
             "-0.0203001022932395859524940805242724601067259796418104502230488160090068057",
             "+0.0103631562402464237307199458065599778725040265335003155034035217850387283",
             "-0.00488719292803767146341420376578964407137592740187386736001668777722125914",
             "+0.00135556105548506177551787075080010874364457167567922814822699418129216854"];
            ["+0.0822479226128438738077716511411289215554396236655129640997987178532804788",
             "+0.196032162333245006055759781563801382788760315019545751706035725357898777",
             "+0.116978483643172761847467585997387748702913901442302633827906329990978550",
             "-0.0204825277456560976298590118654006438219872235078237396303198827234514190",
             "+0.00798999189966233579720442148030827079362758194981571971201211474213850758",
             "-0.00207562578486633419359528915758164772805951363936400480613058750798236207"];
            ["+0.0877378719744515067137433602439480944914707643813859960054483074074552416",
             "+0.172390794624406967987712335438549785037133364423557021432857505376431512",
             "+0.254439495032001621324794183860176141227815026392429007286132542705408521",
             "+0.116978483643172761847467585997387748702913901442302633827906329990978551",
             "-0.0156513758091757022708430246449433269579993686461730105611661052393287579",
             "+0.00341432357674129871237641994523752520797162707650902709951900204619240137"];
            ["+0.0843066851341001107446302003355663380197666790663427630510907257181807085",
             "+0.185267979452106975248330960684647699902136873775246608504886307895791278",
             "+0.223593811046099099964215226188215519533323776351104952152409138196918370",
             "+0.254257069579585109647429252519047957512553782526415717878861475990963909",
             "+0.0901903932620346518924583784594290279153804731866863705724348100592850092",
             "-0.00701124524079369066636422067695386937993730270821137675241910493504753033"];
            ["+0.0864750263608499346324472067379289868303152022069955977178609202476442335",
             "+0.177526353208969968653987471088742034297089951090896196754060639999536273",
             "+0.239625825335829035595856428410992003967971316351739632802026053203292742",
             "+0.224631916579867772503496287486772348817539639661788864354647175338195299",
             "+0.195144512521266716260289347979378707272761134173713625208021159524901768",
             "+0.0428311230947925862600740355431832233817056253710109955996588599497364397"]
        ]), (6,6))'

    b = @dec128 [
            "0.0856622461895851725201480710863664467634112507420219911993177198994728803",
            "0.180380786524069303784916756918858055830760946373372741144869620118570019",
            "0.233956967286345523694935171994775497405827802884605267655812659981957101",
            "0.233956967286345523694935171994775497405827802884605267655812659981957101",
            "0.180380786524069303784916756918858055830760946373372741144869620118570019",
            "0.0856622461895851725201480710863664467634112507420219911993177198994728804"
        ]

    c = @dec128 [
            "0.0337652428984239860938492227530026954326171311438550875637251917366932495891",
            "0.169395306766867743169300202490047326496775717802414964592736647073908251704",
            "0.380690406958401545684749139159644032290694684929989324909302417712862532843",
            "0.619309593041598454315250860840355967709305315070010675090697582287137467157",
            "0.830604693233132256830699797509952673503224282197585035407263352926091748261",
            "0.966234757101576013906150777246997304567382868856144912436274808263306750376"
        ]

    o = 12

    CoefficientsRK(T, :GLRK6, o, a, b, c)
end


function getCoefficientsGLRK(s::Int; T=Float64, high_precision=true)

    if high_precision == true
        if s == 1
            return getCoefficientsGLRK1(T)
        elseif s == 2
            return getCoefficientsGLRK2(T)
        elseif s == 3
            return getCoefficientsGLRK3(T)
        elseif s == 4
            return getCoefficientsGLRK4(T)
        elseif s == 5
            return getCoefficientsGLRK5(T)
        elseif s == 6
            return getCoefficientsGLRK6(T)
        end
    end


    # order
    o = 2s

    # obtain Gauss-Legendre nodes and weights
    q = GaussLegendreQuadrature(s, T)

    # create Lagrange polynomial
    vdm = vandermonde_matrix_inverse(nodes(q))

    # compute monomial basis functions and corresponding integrals
    poly_ints = []
    for i in 1:s
        y = zeros(s)
        y[i] = 1
        mon = *(vdm, y)
        push!(poly_ints, polyint(Poly(mon)))
    end

    # compute Runge-Kutta coefficients
    a = zeros(s,s)
    for i in 1:s
        for j in 1:s
            a[i,j] = poly_ints[j](nodes(q)[i])
        end
    end

    CoefficientsRK(Symbol("GLRK", s), o, a, weights(q), nodes(q))
end
