# order 2
h = [1//2, 1, 1, 1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0]
d1 = [-3//2, 2, -1//2, 0, 0, 0, 0]
d2 = [1//1, -2, 1, 0, 0, 0, 0]
d3 = [-1//1, 3, -3, 1, 0, 0, 0]
H = Diagonal(h)

q1 = [0 1//2 0 0 0 0 0]
q2 = [-1//2 0 1//2 0 0 0 0]
q3 = [0 -1//2 0 1//2 0 0 0]
q4 = [0 0 -1//2 0 1//2 0 0]
q5 = [0 0 0 -1//2 0 1//2 0]
q6 = [0 0 0 0 -1//2 0 1//1]
q7 = [0 0 0 0 0 -1//2 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7)

m1 = [1//1 -1 0 0 0 0 0]
m2 = [m1[2] 2 -1 0 0 0 0]
m3 = [0 -1 2 -1 0 0 0]
m4 = [0 0 -1 2 -1 0 0]
m5 = [0 0 0 -1 2 -1 0]
m6 = [0 0 0 0 -1 2 -1]
m7 = [0 0 0 0 0 -1 2 ]
M = vcat(m1, m2, m3, m4, m5, m6, m7)

r1 = [0 -13//16 7//8 -1//16 0 0 0]
r2 = [-r1[2] 0 -23//16 5//8 0 0 0]
r3 = [-r1[3] -r2[3] 0 -17//16 1//2 0 0]
r4 = [0 -1//2 2//2 0 -2//2 1//2 0]
r5 = [0 0 -1//2 2//2 0 -2//2 1//2]
r6 = [0 0 0 -1//2 2//2 0 -2//2]
r7 = [0 0 0 0 -1//2 2//2 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7)

n1 = [13//10 -12//5 9//10 1//5 0 0 0]
n2 = [n1[2] 26//5 -16//5 2//5 0 0 0]
n3 = [n1[3] n2[3] 47//10 -17//5 1 0 0]
n4 = [n1[4] n2[4] n3[4] -29//5 -4 1 0]
n5 = [0 0 1//1 -4 6 -4 1]
n6 = [0 0 0 1//1 -4 6 -4]
n7 = [0 0 0 0 1//1 -4 6]
N = vcat(n1, n2, n3, n4, n5, n6, n7)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)

x=1//1:length(h)
ks = 0:5
for k in ks
    u=x.^k
    println("k = ", k)
    println(D1*u - k*x.^(k-1))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D2*u - k*(k-1)*x.^(k-2))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D3*u - k*(k-1)*(k-2)*x.^(k-3))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(D4*u - k*(k-1)*(k-2)*(k-3)*x.^(k-4))
end


# order 4
h = [35809//100800,
    13297//11200,
    5701//5600,
    45109//50400,
    35191//33600,
    33503//33600,
    1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
d1 = [-11//6, 18//6, -9//6, 2//6, 0, 0, 0, 0, 0, 0]
d2 = [2//1, -5, 4, -1, 0, 0, 0, 0, 0, 0]
d3 = [-1//1, 3, -3, 1, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = [0 526249//907200 -10819//777600 -50767//907200 -631//28800 91//7776 0 0 0 0]
q2 = [-q1[2] 0 1421209//2721600 16657//201600 -8467//453600 -33059//5443200 0 0 0 0]
q3 = [-q1[3] -q2[3] 0 631187//1360800 400139//5443200 -8789//302400 0 0 0 0]
q4 = [-q1[4] -q2[4] -q3[4] 0 496403//907200 -308533//5443200 0 0 0 0]
q5 = [-q1[5] -q2[5] -q3[5] -q4[5] 0 1805647//2721600 -1//12 0 0 0]
q6 = [-q1[6] -q2[6] -q3[6] -q4[6] -q5[6] 0 8//12 -1//12 0 0]
q7 = [0 0 0 0 1//12 -8//12 0 8//12 -1//12 0]
q8 = [0 0 0 0 0 1//12 -8//12 0 8//12 -1//12]
q9 = [0 0 0 0 0 0 1//12 -8//12 0 8//12]
q10 = [0 0 0 0 0 0 0 1//12 -8//12 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10)

m1 = [2386127//2177280 -515449//453600 -10781//777600 61567//1360800 6817//403200 -1069//136080 0 0 0 0]
m2 = [m1[2] 4756039//2177280 -1270009//1360800 -3751//28800 3067//680400 119459//10886400 0 0 0 0]
m3 = [m1[3] m2[3] 111623//60480 -555587//680400 -551339//5443200 8789//453600 0 0 0 0]
m4 = [m1[4] m2[4] m3[4] 1025327//544320 -464003//453600 222133//5443200 0 0 0 0]
m5 = [m1[5] m2[5] m3[5] m4[5] 5074159//2177280 -1784047//1360800 1//12 0 0 0]
m6 = [m1[6] m2[6] m3[6] m4[6] m5[6] 1812749//725760 -16//12 1//12 0 0]
m7 = [0 0 0 0 1//12 -16//12 30//12 -16//12 1//12 0]
m8 = [0 0 0 0 0 1//12 -16//12 30//12 -16//12 1//12]
m9 = [0 0 0 0 0 0 1//12 -16//12 30//12 -16//12]
m10 = [0 0 0 0 0 0 0 1//12 -16//12 30//12]
M = vcat(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)

r1 = [0 -88471//67200 58139//33600 -1167//2800 -89//11200 7//640 0 0 0 0]
r2 = [-r1[2] 0 -43723//16800 46783//33600 -191//3200 -1567//33600 0 0 0 0]
r3 = [-r1[3] -r2[3] 0 -4049//2400 29083//33600 -71//1400 0 0 0 0]
r4 = [-r1[4] -r2[4] -r3[4] 0 -8591//5600 10613//11200 -1//8 0 0 0]
r5 = [-r1[5] -r2[5] -r3[5] -r4[5] 0 -108271//67200 8//8 -1//8 0 0]
r6 = [-r1[6] -r2[6] -r3[6] -r4[6] -r5[6] 0 -13//8 8//8 -1//8 0]
r7 = [0 0 0 1//8 -8//8 13//8 0 -13//8 8//8 -1//8]
r8 = [0 0 0 0 1//8 -8//8 13//8 0 -13//8 8//8]
r9 = [0 0 0 0 0 1//8 -8//8 13//8 0 -13//8]
r10 = [0 0 0 0 0 0 1//8 -8//8 13//8 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10)

n1 = [4596181//1814400 -10307743//1814400 160961//43200 -535019//907200 109057//1814400 -29273//604800 0 0 0 0]
n2 = [n1[2] 8368543//604800 -9558943//907200 2177057//907200 -11351//86400 204257//1814400 0 0 0 0]
n3 = [n1[3] n2[3] 4938581//453600 -786473//151200 1141057//907200 -120619//907200 0 0 0 0]
n4 = [n1[4] n2[4] n3[4] 3146581//453600 -4614143//907200 24587//14400 -1//6 0 0 0]
n5 = [n1[5] n2[5] n3[5] n4[5] 185709//22400 -11293343//1814400 12//6 -1//6 0 0]
n6 = [n1[6] n2[6] n3[6] n4[6] n5[6] 16787381//1814400 -39//6 12//6 -1//6 0]
n7 = [0 0 0 -1//6 12//6 -39//6 56//6 -39//6 12//6 -1//6]
n8 = [0 0 0 0 -1//6 12//6 -39//6 56//6 -39//6 12//6]
n9 = [0 0 0 0 0 -1//6 12//6 -39//6 56//6 -39//6]
n10 = [0 0 0 0 0 0 -1//6 12//6 -39//6 56//6]
N = vcat(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)

x=1//1:length(h)
ks = 0:7
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D1*u - k*x.^(k-1)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D2*u - k*(k-1)*x.^(k-2)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D3*u - k*(k-1)*(k-2)*x.^(k-3)))
end
for k in ks
    u=x.^k
    println("k = ", k)
    println(Float64.(D4*u - k*(k-1)*(k-2)*(k-3)*x.^(k-4)))
end


# order 6
h = [318365//1016064,
    145979//103680,
    139177//241920,
    964969//725760,
    593477//725760,
    52009//48384,
    141893//145152,
    1019713//1016064,
    1, 1, 1, 1]
e = [1//1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
d1 = [-25//12, 48//12, -36//12, 16//12, -3//12, 0, 0, 0, 0, 0, 0, 0]
d2 = [35//12, -104//12, 114//12, -56//12, 11//12, 0, 0, 0, 0, 0, 0, 0]
d3 = [-5//2, 18//2, -24//2, 14//2, -3//2, 0, 0, 0, 0, 0, 0, 0]
H = Diagonal(h)

q1 = reshape(
    [0, 1547358409//2421619200, -422423//11211200, -1002751721//8717829120,
    -15605263//484323840, 1023865//24216192, 291943739//21794572800, -24659//2534400,
    0, 0, 0, 0], 1, length(h))
q2 = reshape(
    [-q1[2], 0, 23031829//62899200, 10784027//34594560,
    2859215//31135104, -45982103//345945600, -26681//1182720, 538846039//21794572800,
    0, 0, 0, 0], 1, length(h))
q3 = reshape(
    [-q1[3], -q2[3], 0, 28368209//69189120,
    -9693137//69189120, 1289363//17740800, -39181//5491200, -168647//24216192,
    0, 0, 0, 0], 1, length(h))
q4 = reshape(
    [-q1[4], -q2[4], -q3[4], 0,
    5833151//10644480, 4353179//69189120, 2462459//155675520, -215471//10762752,
    0, 0, 0, 0], 1, length(h))
q5 = reshape(
    [-q1[5], -q2[5], -q3[5], -q4[5], 0, 7521509//13837824, -1013231//11531520, 103152839//8717829120,
    0, 0, 0, 0], 1, length(h))
q6 = reshape(
    [-q1[6], -q2[6], -q3[6], -q4[6], -q5[6], 0, 67795697//98841600, -17263733//151351200,
    1//60, 0, 0, 0], 1, length(h))
q7 = reshape(
    [-q1[7], -q2[7], -q3[7], -q4[7], -q5[7], -q6[7], 0, 1769933569//2421619200,
    -9//60, 1//60, 0, 0], 1, length(h))
q8 = reshape(
    [-q1[8], -q2[8], -q3[8], -q4[8], -q5[8], -q6[8], -q7[8], 0,
    45//60, -9//60, 1//60, 0], 1, length(h))
q9 = [0 0 0 0 0 -1//60 9//60 -45//60 0 45//60 -9//60 1//60]
q10 = [0 0 0 0 0 0 -1//60 9//60 -45//60 0 45//60 -9//60]
q11 = [0 0 0 0 0 0 0 -1//60 9//60 -45//60 0 45//60]
q12 = [0 0 0 0 0 0 0 0 -1//60 9//60 -45//60 0]
Q = vcat(q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12)

m1 = reshape(
    [4347276223//3736212480, -1534657609//1210809600, 68879//3057600, 1092927401//13076743680,
    18145423//968647680, -1143817//60540480, -355447739//65383718400, 56081//16473600,
    0, 0, 0, 0], 1, length(h))
m2 = reshape(
    [m1[2], 42416226217//18681062400, -228654119//345945600, -12245627//34594560,
    -2995295//46702656, 52836503//691891200, 119351//12812800, -634102039//65383718400,
    0, 0, 0, 0], 1, length(h))
m3 = reshape(
    [m1[3], m2[3], 5399287//4193280, -24739409//34594560,
    7878737//69189120, -1917829//31449600, 39727//3660800, 10259//4656960,
    0, 0, 0, 0], 1, length(h))
m4 = reshape(
    [m1[4], m2[4], m3[4], 7780367599//3736212480,
    -70085363//69189120, -500209//6289920, -311543//17962560, 278191//21525504,
    0, 0, 0, 0], 1, length(h))
m5 = reshape(
    [m1[5], m2[5], m3[5], m4[5],
    7116321131//3736212480, -545081//532224, 811631//11531520, -84101639//13076743680,
    0, 0, 0, 0], 1, length(h))
m6 = reshape(
    [m1[6], m2[6], m3[6], m4[6],
    m5[6], 324760747//138378240, -65995697//49420800, 1469203//13759200,
    -2//180, 0, 0, 0], 1, length(h))
m7 = reshape(
    [m1[7], m2[7], m3[7], m4[7],  m5[7], m6[7], 48284442317//18681062400, -1762877569//1210809600,
    27//180, -2//180, 0, 0], 1, length(h))
m8 = reshape(
    [m1[8], m2[8], m3[8], m4[8],  m5[8], m6[8], m7[8], 10117212851//3736212480,
    -270//180, 27//180, -2//180, 0], 1, length(h))
m9 = [0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180 27//180 -2//180]
m10 = [0 0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180 27//180]
m11 = [0 0 0 0 0 0 0 -2//180 27//180 -270//180 490//180 -270//180]
m12 = [0 0 0 0 0 0 0 0 -2//180 27//180 -270//180 490//180]
M = vcat(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12)

r1 = reshape(
    [0, -10882810591//5811886080, 398713069//132088320, -1746657571//1162377216,
    56050639//145297152, -11473393//1162377216, -38062741//1452971520, 30473//4392960,
    0, 0, 0, 0], 1, length(h))
r2 = reshape(
    [-r1[2], 0, -3720544343//830269440, 767707019//207567360,
    -1047978301//830269440, 1240729//14826240, 6807397//55351296, -50022767//1452971520,
    0, 0, 0, 0], 1, length(h))
r3 = reshape(
    [-r1[3], -r2[3], 0, -2870078009//830269440,
    74962049//29652480, -12944857//30750720, -17846623//103783680, 68707591//1162377216,
    0, 0, 0, 0], 1, length(h))
r4 = reshape(
    [-r1[4], -r2[4], -r3[4], 0,
    -727867087//276756480, 327603877//207567360, -175223717//830269440, 1353613//726485760,
    0, 0, 0, 0], 1, length(h))
r5 = reshape(
    [-r1[5], -r2[5], -r3[5], -r4[5],
    0, -1804641793//830269440, 311038417//207567360, -1932566239//5811886080,
    7//240, 0, 0, 0], 1, length(h))
r6 = reshape(
    [-r1[6], -r2[6], -r3[6], -r4[6], -r5[6], 0, -1760949511//830269440, 2105883973//1452971520,
    -72//240, 7//240, 0, 0], 1, length(h))
r7 = reshape(
    [-r1[7], -r2[7], -r3[7], -r4[7], -r5[7], -r6[7], 0, -1081094773//528353280,
    338//240, -72//240, 7//240, 0], 1, length(h))
r8 = reshape(
    [-r1[8], -r2[8], -r3[8], -r4[8], -r5[8], -r6[8], -r7[8], 0,
    -488//240, 338//240, -72//240, 7//240], 1, length(h))
r9 = [0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240 338//240 -72//240]
r10 = [0 0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240 338//240]
r11 = [0 0 0 0 0 0 -7//240 72//240 -338//240 488//240 0 -488//240]
r12 = [0 0 0 0 0 0 0 -7//240 72//240 -338//240 488//240 0]
R = vcat(r1, r2, r3, r4, r5, r6, r7, r8, r9, r10, r11, r12)

n1 = reshape(
    [40833734273//10761070320, -162181998421//16397821440, 4696168417//521748864, -245714671483//68870850048,
    2185939219//2869618752, -15248255797//114784750080, 345156907//12298366080, 6388381//1093188096,
    0, 0, 0, 0], 1, length(h))
n2 = reshape(
    [n1[2], 147281127041//5380535160, -3072614435609//114784750080, 320122985851//28696187520,
    -768046031383//344354250240, 7861605187//14348093760, -803762437//4251287040, 167394281//86088562560,
    0, 0, 0, 0], 1, length(h))
n3 = reshape(
    [n1[3], n2[3], 139712483333//4782697920, -1634124842747//114784750080,
    90855193447//28696187520, -26412188989//38261583360, 668741173//1793511720, -132673781//2342545920,
    0, 0, 0, 0], 1, length(h))
n4 = reshape(
    [n1[4], n2[4], n3[4], 437353997177//43044281280,
    -172873969321//38261583360, 34759553483//28696187520, -98928859751//344354250240, 295000207//3587023440,
    0, 0, 0, 0], 1, length(h))
n5 = reshape(
    [n1[5], n2[5], n3[5], n4[5],
    126711914423//21522140640, -520477408939//114784750080, 49581230003//28696187520, -99640101991//344354250240,
    7//240, 0, 0, 0], 1, length(h))
n6 = reshape(
    [n1[6], n2[6], n3[6], n4[6],
    n5[6], 19422074929//2391348960, -772894368601//114784750080, 10579712849//4099455360,
    -96//240, 7//240, 0, 0], 1, length(h))
n7 = reshape(
    [n1[7], n2[7], n3[7], n4[7],
    n5[7], n6[7], 456715296239//43044281280, -915425403107//114784750080,
    676//240, -96//240, 7//240, 0], 1, length(h))
n8 = reshape(
    [n1[8], n2[8], n3[8], n4[8], n5[8], n6[8], n7[8], 488029542379//43044281280,
    -1952//240, 676//240, -96//240, 7//240], 1, length(h))
n9 = [0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240 676//240 -96//240]
n10 = [0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240 676//240]
n11 = [0 0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240 -1952//240]
n12 = [0 0 0 0 0 0 0 7//240 -96//240 676//240 -1952//240 2730//240]
N = vcat(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12)

D1 = H \ (Q - 1//2*e*e'); display(D1)
D2 = H \ (-M - e*d1'); display(D2)
D3 = H \ (R + 1//2*d1*d1' - e*d2'); display(D3)
D4 = H \ (N - e*d3' + d1*d2'); display(D4)
