using OrthogonalPolynomials
using Test, Base.Math

# Laguerre
# Abramowitz and Stegun pg. 799
const TriangleCoeffsLaguerre = [
    1 1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600;
    1 1 1 -4 -18 -96 -600 -4320 -35280 -322560 -3265920 -36288000 -439084800 -6748019200;
    2 2 -4 2 18 144 1200 10800 105840 1128960 13063680 163296000 2195424000 31614106600;
    6 6 -18 9 6 -96 -1200 -14400 -176400 -2257920 -30481920 -435456000 -6586272000 -105380352000;
    24 24 -96 72 -16 24 600 10800 176400 2822400 45722880 762048000 13172544000 237105792000;
    120 120 -600 600 -200 25 120 -4320 -106840 -2257920 -45722880 -914457600 -18441561600 -379369267200;
    720 720 -4320 5400 -2400 450 -36 720 35280 1128960 30481920 762048000 18441561600 442597478400;
    5040 5040 -35280 52920 -29400 7350 -882 49 5040 -322560 -13063680 -435456000 -13172544000 -379369267200;
    40320 40320 -322560 564480 -376320 117600 -18816 1568 -64 40320 3269520 163296000 6586272000 237105792000;
    362880 362880 -3265920 6531840 -5080320 1905120 -381024 42336 -2592 81 362880 -36288000 -2195424000 -105380352000;
    3623800 3628800 -36288000 81648000 -72576000 31752000 -7620480 1058400 -86400 4050 -100 3628800 439084800 31614106600;
    39916800 39916800 -439084800 1097712000 -1097712000 648806000 -103679680 25613280 -2613600 163350 -6050 121 39916800 -5748019200;
    479001600 479001600 -5748019200 15807052800 -17563392000 9879408000 -3161410560 614718720 -75271680 5880600 -290400 8712 -144 479001600]
L0(x) = 1.0
    L1(x) = -x + 1.0
    L2(x) = .5 * (x^2 - 4x + 2)
    L3(x) = (1/6)*(-x^3 +9x^2 -18x + 6)
    L4(x) = (1/24)*(x^4 - 16x^3 + 72x^2 - 96x + 24)
    L5(x) = (1/120)*(-x^5 + 25x^4 -200x^3 + 600x^2 - 600x + 120)
    L6(x) = (1/720)*(x^6 - 36x^5 + 450x^4 - 2400x^3 + 5400x^2 - 4320x + 720)
    L7(x) = (1/5040)*(-x^7 + 49x^6 - 882x^5 + 7350x^4 - 29_400x^3 + 52_920x^2 - 35_280x + 5040)
    L8(x) = (1/40320)*(x^8 - 64x^7 + 1568x^6 - 18_816x^5 + 117_600x^4 - 376_320x^3 + 564_480x^2 - 322_560x + 40320)
    L9(x) = (1/362880)*(-x^9 + 81x^8 - 2592x^7 + 42_336x^6 - 381_024x^5 + 1_905_120x^4 - 5_080_320x^3 + 6_531_840x^2 - 3_265_920x + 362_880)
    L10(x) = (1/3628800)*(x^10 - 100x^9 + 4_050x^8 - 86_400x^7 + 1_058_400x^6 - 7_620_480x^5 + 31_752_000x^4 - 72_576_000x^3  + 81_648_000x^2 - 36_288_000x + 3_628_800)
    L11(x) = (1/39916800)*(-x^11 + 121x^10 - 6_050x^9 + 163_350x^8 - 2_613_600x^7 + 25_613_280x^6 - 153_679_680x^5 + 548_856_000x^4 - 1_097_712_000x^3 + 1_097_712_000x^2 - 439_084_800x + 39_916_800)
    L12(x) = (1/479001600)*(x^12 - 144x^11 + 8_712x^10 - 290_400x^9 + 5_880_600x^8 - 75_271_680x^7 + 614_718_720x^6 - 3_161_410_560x^5 + 9_879_408_000x^4 - 17_563_392_000x^3 + 15_807_052_800x^2 - 5_748_019_200x + 479_001_600)

    Ls = [L0, L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12]

    # Credit to @yingboma
    function laguerre(n, x)

        n == 0 ? (return 1) : nothing

        l0, l1 = 1, 1-x
        for k in 1:n-1
            l0, l1 = l1, ( (2k+1-x)*l1 - k*l0 )/(k+1)
        end
        l1
    end
#n\z is the first column - Abramowitz and Stegun p. 800
const testvalslaguerre = [0. .5 1 3 5 10];
const TestTableLaguerre = [
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000;
    1 0.5000000000 0.0000000000 -2.0000000000 -4.0000000000 -9.0000000000;
    2 0.1250000000 -0.5000000000 -0.5000000000 3.5000000000 31.0000000000;
    3 -0.1458333333 -0.6666666667 1.0000000000 2.6666666667 -45.6666666667;
    4 -0.3307291667 -0.6250000000 1.3750000000 -1.2916666667 11.0000000000;
    5 -0.4455729167 -0.4666666667 0.8500000000 -3.1666666667 34.3333333333;
    6 -0.5041449653 -0.2569444444 -0.0125000000 -2.0902777778 -3.4444444444;
    7 -0.5183392237 -0.0404761905 -0.7464285714 0.3253968254 -30.9047619048;
    8 -0.4983629984 0.1539930556 -1.1087053571 2.2357390873 -16.3015873016;
    9 -0.4529195204 0.3097442681 -1.0611607143 2.6917438272 14.7918871252;
    10 -0.3893744141 0.4189459325 -0.7000223214 1.7562761795 27.9841269841;
    11 -0.3139072988 0.4801341791 -0.1807995130 0.1075436909 14.5369568703;
    12 -0.2316496389 0.4962122235 0.3403546063 -1.4486042948 -9.9037464593];

# Hermite - Obtained via the HTable function in utils.jl
const TriangleCoeffsHermite =[
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    -2 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 -12 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    12 0 -48 0 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 120 0 -160 0 32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; # H5
    -120 0 720 0 -480 0 64 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 -1680 0 3360 0 -1344 0 128 0 0 0 0 0 0 0 0 0 0 0 0 0;
    1680 0 -13440 0 13440 0 -3584 0 256 0 0 0 0 0 0 0 0 0 0 0 0;
    0 30240 0 -80640 0 48384 0 -9216 0 512 0 0 0 0 0 0 0 0 0 0 0;
    -30240 0 302400 0 -403200 0 161280 0 -23040 0 1024 0 0 0 0 0 0 0 0 0 0; #H10
    0 -665280 0 2217600 0 -1774080 0 506880 0 -56320 0 2048 0 0 0 0 0 0 0 0 0;
    665280 0 -7983360 0 13305600 0 -7096320 0 1520640 0 -135168 0 4096 0 0 0 0 0 0 0 0; # H12
    0 17297280 0 -69189120 0 69189120 0 -263577600 0 4392960 0 -319488 0 8192 0 0 0 0 0 0 0;
    -17297280 0 242161920 0 -484323840 0 322882560 0 -92252160 0 12300288 0 -745472 0 16384 0 0 0 0 0 0;
    0 -5189184000 0 2421619200 0 -2905943040 0 1383782400 0 -307507200 0 33546240 0 -1720320 0 32768 0 0 0 0 0;
    518918400 0 -8302694400 0 19372953600 0 -15498362880 0 5535129600 0 -984023040 0 89456640 0 -3932160 0 65536 0 0 0 0;
    0 17643225600 0 -94097203200 0 131736084480 0 -75277762560 0 20910489600 0 -3041525760 0 233963520 0 -8912896 0 131072 0 0 0;
    -17643225600 0 317578060800 0 -846874828800 0 790416506880 0 -338749931520 0 75277762560 0 -9124577280 0 601620480 0 -20054016 0 262144 0 0;
    0 -670442572800 0 4022655436800 0 -6436248698880 0 4290832465920 0 -1430277488640 0 260050452480 0 -26671841280 0 1524105216 0 -44826624 0 524288 0;
    670442572800 0 -13408851456000 0 40226554368000 0 -42908324659200 0 21454162329600 0 -5721109954560 0 866834841600 0 -76205260800 0 3810263040 0 -99614720 0 1048576]
H0(x) = 1.0
    H1(x) = 2x
    H2(x) = -2 + 4x^2
    H3(x) = -12x + 8x^3
    H4(x) = 12 - 48x^2 + 16x^4
    H5(x) = 120x -160x^3 + 32x^5
    H6(x) = -120 + 720x^2 - 480x^4 + 64x^6
    H7(x) = -1680x + 3360x^3 - 1344x^5 + 128x^7
    H8(x) =1680 - 13440x^2 + 13440x^4 - 3584x^6 + 256x^8
    H9(x) =30240x - 80640x^3 + 48384x^5 - 9216x^7 + 512x^9
    H10(x) =-30240 + 302400x^2 - 403200x^4 + 161280x^6 - 23040x^8 + 1024x^10
    H11(x) =-665280x + 2217600x^3 - 1774080x^5 + 506880x^7 - 56320x^9 + 2048x^11
    H12(x) =665280 - 7983360x^2 + 13305600x^4 - 7096320x^6 + 1520640x^8 - 135168x^10 + 4096x^12

    Hs = [H0, H1, H2, H3, H4, H5, H6, H7, H8, H9, H10, H11, H12]

    function hermite(x,n)
        if n == 0
            return 1
        elseif n == 1
            return 2x
        else
            return 2x*hermite(x,n-1) - 2*(n-1)*hermite(x,n-2)
        end
    end

# A and S p. 802
const testvalshermite = [0 0.5 1 3 5 10]
const TestTableHermite = [
    0 1.00000 1.00000 1.0000000 1.0000000000 1.0000000000;
    1 1.00000 2.00000 6.0000000 1.0000000000e1 2.0000000000e1
    2 -1.00000 2.00000 3.4000000e1 9.8000000000e1 3.9800000000e2;
    3 -5.00000 -4.00000 1.8000000e2 9.4000000000e2 7.8800000000e3;
    4 1.00000 -2.00000e1 8.7600000e2 8.8120000000e3 1.5521200000e5;
    5 4.10000e1 -8.00000 3.8160000e3 8.0600000000e4 3.0412000000e6;
    6 3.10000e1 1.84000e2 1.4136000e4 7.1788000000e5 5.9271880000e7;
    7 -4.61000e2 4.64000e2 3.9024000e4 6.2116000000e6 1.1489432000e9;
    8 -8.95000e2 -1.64800e3 3.6240000e4 5.2065680000e7 2.2149057680e10;
    9 6.48100e3 -1.07200e4 -4.0694400e5 4.2127120000e8 4.2459806240e11;
    10 2.25910e4 8.22400e3 -3.0939840e6 3.2755297600e9 8.0932782098e12;
    11 -1.07029e5 2.30848e5 -1.0425024e7 2.4329873600e10 1.5337360295e14;
    12 -6.04031e5 2.80768e5 5.5175040e6 1.7123708128e11 2.8894199383e15];

# Legendre - Table made with script from utils.jl
const TriangleCoeffsLegendre = [
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    6 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    20 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    70 60 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    252 280 60 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; # P5
    924 1260 420 20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    3432 5544 2520 280 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    12870 24024 13860 2520 70 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    48620 102960 72072 18480 1260 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    184756 437580 360360 120120 13860 252 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; #P10
    705432 1847560 1750320 720720 120120 5544 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    2704156 7759752 8314020 4084080 900900 72072 924 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    10400600 32449872 38798760 22170720 6126120 720720 24024 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    40116600 135207800 178474296 11639628038798760 6126120 360360 3432 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    155117520 561632400 811246800 594914320 232792560 46558512 4084080 102960 0 0 0 0 0 0 0 0 0 0 0 0 0; # P15
    601080390 2326762800 3650610600 2974571600 1338557220 325909584 38798760 1750320 12870 0 0 0 0 0 0 0 0 0 0 0 0;
    23336062209617286240 16287339600 14602442400 7436429000 2141691552 325909584 22170720 437580 0 0 0 0 0 0 0 0 0 0 0 0 0 ;
    9075135300 39671305740 72129646800 70578471600 40156716600 13385572200 2498640144 232792560 8314020 48620 0 0 0 0 0 0 0 0 0 0 0;
    35345263800 163352435400 317370445920 336605018400 211735414800 80313433200 17847429600 2141691552116396280 1847560 0 0 0 0 0 0 0 0 0 0 0 0 ;
    137846528820 671560012200 1388495700900 1586852229600 1093966309800 465817912560 120470149800 17847429600 1338557220 38798760 184756 0 0 0 0 0 0 0 0 0 0]
P0(x) = 1
    P1(x) =  x
    P2(x) = .5*(3x^2 -1)
    P3(x) = .5*(5x^3 -3x)
    P4(x) = (1/8)*(35x^4 - 30x^2 + 3)
    P5(x) = (1/8)*(63x^5 - 70x^3 + 15x)
    P6(x) = (1/16)*(231x^6 - 315x^4 + 105x^2 - 5)
    P7(x) = (1/16)*(429x^7 - 693x^5 + 315x^3 - 35x)
    P8(x) = (1/128)*(6435x^8 - 12012x^6 + 6930x^4 - 1260x^2 + 35)
    P9(x) = (1/128)*(12155x^9 - 25740x^7 + 18018x^5 - 4620x^3 + 315x)
    P10(x) = (1/256)*(46189x^10 - 109395x^8 + 90090x^6- 30030x^4 + 3465x^2 - 63)
    P11(x) = (1/256)*(2^11)*(705432x^11 - 1847560x^9 + 1750320x^7 - 720720x^5 + 120120x^3 - 5544x)
    P12(x) = (1/512)*(2^12)*(2704156x^12 - 7759752x^10 + 8314020x^8 - 4084080x^6 + 900900x^4 - 72072x^2 + 924)

    Ps = [P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12]
    testPs = [P2, P3, P9, P10]

    function legendre(x, n)
        if n == 0
            return 1
        elseif n == 1
            return x
        else
            return ((2n-1)*x*legendre(x,n-1) - (n - 1)*legendre(x, n-2))/n
        end
    end

# Credit: Abramowitz and Stegun, Table 8.1, p. 342
# X arccos x P2(x) P3(x) P9(x) P10(x)
const TestTableLegendre =[
    0.00 90.0000000 -0.50000  0.0000000 0.00000000 -0.2460937;
    0.01 89.4270326 -0.49985 -0.0149975 0.02457330 -0.2447414;
    0.02 88.8540080 -0.49940 -0.0299800 0.04893045 -0.2406984;
    0.03 88.2808687 -0.49865 -0.0449325 0.07285701 -0.2340069;
    0.04 87.7075572 -0.49760 -0.0598400 0.09614188 -0.2247364;
    0.05 87.1340160 -0.49625 -0.0746875 0.11857899 -0.2129835;
    0.06 86.5601872 -0.49460 -0.0894600 0.13996890 -0.1988711;
    0.07 85.9860128 -0.49265 -0.1041425 0.16012040 -0.1825468;
    0.08 85.4114343 -0.49040 -0.1187200 0.17885206 -0.1641820;
    0.09 84.8363929 -0.48785 -0.1331775 0.19599366 -0.1439702;
    0.10 84.2608295 -0.48500 -0.1475000 0.21138764 -0.1221250;
    0.11 83.6846844 -0.48185 -0.1616725 0.22489042 -0.0988786;
    0.12 83.1078974 -0.47840 -0.1756800 0.23637363 -0.0744793;
    0.13 82.5304077 -0.47465 -0.1895075 0.24572526 -0.0491890;
    0.14 81.9521537 -0.47060 -0.2031400 0.25285070 -0.0232812;
    0.15 81.3730734 -0.46625 -0.2165625 0.25767367 0.0029618;
    0.16 80.7931038 -0.46160 -0.2297600 0.26013706 0.0292520;
    0.17 80.2121810 -0.45665 -0.2427175 0.26020358 0.0552981;
    0.18 79.6302402 -0.45140 -0.2554200 0.25785632 0.0808085;
    0.19 79.0472158 -0.44585 -0.2678525 0.25309918 0.1054942;
    0.20 78.4630410 -0.44000 -0.2800000 0.24595712 0.1290720;
    0.21 77.8776477 -0.43385 -0.2918475 0.23647631 0.1512674;
    0.22 77.2909670 -0.42740 -0.3033800 0.22472407 0.1718175;
    0.23 76.7029282 -0.42065 -0.3145825 0.21078870 0.1904736;
    0.24 76.1134596 -0.41360 -0.3254400 0.19477914 0.2070049;
    0.25 75.5224878 -0.40625 -0.3359375 0.17682442 0.2212002;
    0.26 74.9299379 -0.39860 -0.3460600 0.15707305 0.2328714;
    0.27 74.3357331 -0.39065 -0.3557925 0.13569215 0.2418552;
    0.28 73.7397953 -0.38240 -0.3651200 0.11286642 0.2480162;
    0.29 73.1420440 -0.37385 -0.3740275 0.08879707 0.2512481;
    0.30 72.5423969 -0.36500 -0.3825000 0.06370038 0.2514763;
    0.31 71.9407695 -0.35585 -0.3905225 0.03780634 0.2486591;
    0.32 71.3370751 -0.34640 -0.3980800 0.01135691 0.2427889;
    0.33 70.7312245 -0.33665 -0.4051575 -0.01539566 0.2338937;
    0.34 70.1231259 -0.32660 -0.4117400 -0.04219085 0.2220373;
    0.35 69.5126849 -0.31625 -0.4178125 -0.06876185 0.2073200;
    0.36 68.8998039 -0.30560 -0.4233600 -0.09483780 0.1898783;
    0.37 68.2843827 -0.29465 -0.4283675 -0.12014608 0.1698848;
    0.38 67.6663173 -0.28340 -0.4328200 -0.14441472 0.1475472;
    0.39 67.0455006 -0.27185 -0.4367025 -0.16737489 0.1231073;
    0.40 66.4218215 -0.26000 -0.4400000 -0.18876356 0.0968391;
    0.41 65.7951652 -0.24785 -0.4426975 -0.20832609 0.0690471;
    0.42 65.1654125 -0.23540 -0.4447800 -0.22581900 0.0400639;
    0.43 64.5324399 -0.22265 -0.4462325 -0.24101269 0.0102469;
    0.44 63.8961188 -0.20960 -0.4470400 -0.25369426 -0.0200245;
    0.45 63.2563161 -0.19625 -0.4471875 -0.26367022 -0.0503530;
    0.46 62.6128925 -0.18260 -0.4466600 -0.27076932 -0.0803272;
    0.47 61.9657035 -0.16865 -0.4454425 -0.27484521 -0.1095264;
    0.48 61.3145980 -0.15440 -0.4435200 -0.27577908 -0.1375251;
    0.49 60.6594184 -0.13985 -0.4408775 -0.27348225 -0.1638987;
    0.50 60.0000000 -0.12500 -0.4375000 -0.26789856 -0.1882286;
    0.51 59.3361703 -0.10985 -0.4333725 -0.25900667 -0.2101083;
    0.52 58.6677485 -0.09440 -0.4284800 -0.24682215 -0.2291492;
    0.53 57.9945451 -0.07865 -0.4228075 -0.23139939 -0.2449873;
    0.54 57.3163611 -0.06260 -0.4163400 -0.21283321 -0.2572892;
    0.55 56.6329870 -0.04625 -0.4090625 -0.19126025 -0.2657585;
    0.56 55.9442022 -0.02960 -0.4009600 -0.16686000 -0.2701428;
    0.57 55.2497742 -0.01265 -0.3920175 -0.13985552 -0.2702397;
    0.58 54.5494574  0.00460 -0.3822200 -0.11051366 -0.2659030;
    0.59 53.8429918 0.02215 -0.3715525 -0.07914497 -0.2570492;
    0.60 53.1301024 0.04000 -0.3600000 -0.04610304 -0.2436627;
    0.61 52.4104970 0.05815 -0.3475475 -0.01178332 -0.2258016;
    0.62 51.6838655 0.07660 -0.3341800 0.02337862 -0.2036019;
    0.63 50.9498775 0.09535 -0.3198825 0.05890951 -0.1772816;
    0.64 50.2081805 0.11440 -0.3046400 0.09430141 -0.1471441;
    0.65 49.4583981 0.13375 -0.2884375 0.12901554 -0.1135805;
    0.66 48.7001272 0.15340 -0.2712600 0.16248693 -0.0770701;
    0.67 47.9329352 0.17335 -0.2530925 0.19412981 -0.0381808;
    0.68 47.1563569 0.19360 -0.2339200 0.22334410 0.0024330;
    0.69 46.3698911 0.21415 -0.2137275 0.24952270 0.0440337;
    0.70 45.5729960 0.23500 -0.1925000 0.27205993 0.0858058;
    0.71 44.7650847 0.25615 -0.1702225 0.29036111 0.1268631;
    0.72 43.9455196 0.27760 -0.1468800 0.30385323 0.1662589;
    0.73 43.1136059 0.29935 -0.1224575 0.31199698 0.2029976;
    0.74 42.2685844 0.32140 -0.0969400 0.31430004 0.2360508;
    0.75 41.4096221 0.34375 -0.0703125 0.31033185 0.2643745;
    0.76 40.5358021 0.36640 -0.0425600 0.29973981 0.2869319;
    0.77 39.6461111 0.38935 -0.0136675 0.28226712 0.3027179;
    0.78 38.7394246 0.41260 0.0163800 0.25777224 0.3107893;
    0.79 37.8144885 0.43615 0.0475975 0.22625012 0.3102979;
    0.80 36.8698976 0.46000 0.0800000 0.18785528 0.3005298;
    0.81 35.9040686 0.48415 0.1136025 0.14292678 0.2809487;
    0.82 34.9152062 0.50860 0.1484200 0.09201529 0.2512452;
    0.83 33.9012620 0.53335 0.1844675 0.03591226 0.2113919;
    0.84 32.8598804 0.55840 0.2217600 -0.02431874 0.1617050;
    0.85 31.7883306 0.58375 0.2603125 -0.08730820 0.1029123;
    0.86 30.6834171 0.60940 0.3001400 -0.15134456 0.0362291;
    0.87 29.5413605 0.63535 0.3412575 -0.21433544 -0.0365586;
    0.88 28.3576366 0.66160 0.3836800 -0.27376627 -0.1130029;
    0.89 27.1267531 0.68815 0.4274225 -0.32665610 -0.1898929;
    0.90 25.8419328 0.71500 0.4725000 -0.36951049 -0.2631456;
    0.91 24.4946485 0.74215 0.5189275 -0.39827146 -0.3276858;
    0.92 23.0739181 0.76960 0.5667200 -0.40826421 -0.3773158;
    0.93 21.5651850 0.79735 0.6158925 -0.39414060 -0.4045743;
    0.94 19.9484436 0.82540 0.6664600 -0.34981919 -0.4005829;
    0.95 18.1948723 0.85375 0.7184375 -0.26842182 -0.3548803;
    0.96 16.2602047 0.88240 0.7718400 -0.14220642 -0.2552433; #this last 3 is a 4 in the book :(
    0.97 14.0698677 0.91135 0.8266825 0.03750397 -0.0874940;
    0.98 11.4783409 0.94060 0.8829800 0.28039609 0.1647081;
    0.99  8.1096144 0.97015 0.9407475 0.59724553 0.5200890;
    1.00  0.0000000 1.00000 1.0000000 1.00000000 1.0000000;]

# Chebyshev T First kind - Credit to Michael Somos, https://oeis.org/A049310
const TriangleCoeffsChebyshevT = float([
    1    0    0     0     0      0      0      0      0      0      0       0     0; # T0
    0    1    0     0     0      0      0      0      0      0      0       0     0;
   -1    0    2     0     0      0      0      0      0      0      0       0     0;
    0   -3    0     4     0      0      0      0      0      0      0       0     0;
    1    0   -8     0     8      0      0      0      0      0      0       0     0;
    0    5    0   -20     0     16      0      0      0      0      0       0     0; # T5
   -1    0   18     0   -48      0     32      0      0      0      0       0     0;
    0   -7    0    56     0   -112      0     64      0      0      0       0     0;
    1    0  -32     0   160      0   -256      0    128      0      0       0     0;
    0    9    0  -120     0    432      0   -576      0    256      0       0     0;
   -1    0   50     0  -400      0   1120      0  -1280      0    512       0     0; # T10
    0  -11    0   220     0  -1232      0   2816      0  -2816      0    1024     0;
    1    0  -72     0   840      0  -3584      0   6912      0  -6144       0  2048;])
T0(x) = 1
    T1(x) = x
    T2(x) = -1 + 2x^2
    T3(x) = -3x + 4x^3
    T4(x) = 1 - 8x^2 + 8x^4
    T5(x) = 5x - 20x^3 + 16x^5
    T6(x) = -1 + 18x^2 - 48x^4 + 32x^6
    T7(x) = -7x + 56x^3 - 112x^5 + 64x^7
    T8(x) = 1 - 32x^2 + 160x^4 - 256x^6 + 128x^8
    T9(x) = 9x - 120x^3 + 432x^5 - 576x^7 + 256x^9
    T10(x) = -1 + 50x^2 - 400x^4 + 1120x^6 - 1280x^8 + 512x^10
    T11(x) = -11x + 220x^3 - 1232x^5 + 2816x^7 - 2816x^9 + 1024x^11
    T12(x) = 1 - 72x^2 + 840x^4 - 3584x^6 + 6912x^8 - 6144x^10 + 2048x^12

    Ts = [T0 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10 T11 T12]

    function chebyshevfirst(x, n)
        if n == 0
            1
        elseif n == 1
            x
        else
            2x*chebyshevfirst(x,n-1) - chebyshevfirst(x, n-2)
        end
    end


# Credit: Abramowitz and Stegun
const testvalschebyshev = [0.2 0.4 0.6 0.8 1.0]
const TestTableChebyshevT = float([
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1;
    1 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1;
    2 -0.920000000 -0.680000000 -0.280000000 0.2800000000 1;
    3 -0.568000000 -0.944000000 -0.936000000 -0.352000000 1;
    4 0.6928000000 -0.075200000 -0.843200000 -0.843200000 1;
    5 0.8451200000 0.8838400000 -0.075840000 -0.997120000 1;
    6 -0.354752000 0.7822720000 0.7521920000 -0.752192000 1;
    7 -0.987020800 -0.258022400 0.9784704000 -0.206387200 1;
    8 -0.040056320 -0.988689920 0.4219724800 0.4219724800 1;
    9 0.9709982720 -0.532929536 -0.4721034240 0.881543168 1;
    10 0.4284556288 0.5623462912 -0.9884965888 0.9884965888 1;
    11 -0.7996160205 0.9828065690 -0.7140924826 0.7000513741 1;
    12 -0.7483020370 0.2238989640 0.1315856097 0.1315856097 1])

# Chebyshev U Second kind
const TriangleCoeffsChebyshevU = [
    1    0    0     0     0      0      0      0      0      0       0     0     0; # T0
    0    2    0     0     0      0      0      0      0      0       0     0     0;
    -1   0    4     0     0      0      0      0      0      0       0     0     0;
    0   -4    0     8     0      0      0      0      0      0       0     0     0;
    1    0  -12     0    16      0      0      0      0      0       0     0     0;
    0    6    0   -32     0     32      0      0      0      0       0     0     0; # T5
    -1   0   24     0   -80      0     64      0      0      0       0     0     0;
    0   -8    0    80     0   -192      0    128      0      0       0     0     0;
    1    0  -40     0   240      0   -448      0    256      0       0     0     0;
    0   10    0  -160     0    672      0  -1024      0    512       0     0     0;
    -1   0   60     0  -560      0   1792      0  -2304      0    1024     0     0; # T10
    0  -12    0   280     0  -1792      0   4608      0  -5120       0  2048     0;
    1    0  -84     0  1120      0  -5376      0  11520      0  -11264     0  4096;]
U0(x) = 1
    U1(x) = 2x
    U2(x) = 4x^2 - 1
    U3(x) = 8x^3 - 4x
    U4(x) = 16x^4 - 12x^2 + 1
    U5(x) = 32x^5 - 32x^3 + 6x
    U6(x) = 64x^6  - 80x^4 + 24x^2 - 1
    U7(x) = 128x^7 - 192x^5 + 80x^3 - 8x
    U8(x) = 256x^8 - 448x^6 + 240x^4 - 40x^2 + 1
    U9(x) = 512x^9 - 1024x^7 + 672x^5 - 160x^3 + 10x
    U10(x) = 1024x^10 - 2304x^8 + 1792x^6 - 560x^4 + 60x^2 - 1
    U11(x) = 2048x^11 - 5120x^9 + 4608x^7 - 1792x^5 + 280x^3 - 12x
    U12(x) = 4096x^12 - 11264x^10 + 11520x^8 - 5376x^6 + 1120x^4 - 84x^2 + 1

    Us = [ U0 U1 U2 U3 U4 U5 U6 U7 U8 U9 U10 U11 U12]

    function chebyshevsecond(x, n)
        if n == 0
            1
        elseif n == 1
            2x
        else
            2x*chebyshevsecond(x,n-1) - chebyshevsecond(x, n-2)
        end
    end

# Credit: Abramowitz and Stegun
const TestTableChebyshevU = [
    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1;
    1 0.4000000000 0.8000000000 1.2000000000 1.6000000000 2;
    2 -0.840000000 -0.360000000 0.4400000000 1.5600000000 3;
    3 -0.736000000 -1.088000000 -0.672000000 0.8960000000 4;
    4 0.5456000000 -0.510400000 -1.246400000 -0.126400000 5;
    5 0.9542400000 0.6796800000 -0.823680000 -1.098240000 6;
    6 -0.163904000 1.0541440000 0.2579840000 -1.630784000 7;
    7 -1.019801600 0.1636352000 1.1332608000 -1.511014400 8;
    8 -0.244016640 -0.923235840 1.1019289600 -0.786839040 9;
    9 0.9221949440 -0.902223872 0.1890539520 0.2520719360 10;
    10 0.6128946176 0.2014567424 -0.8750642176 1.1901541376 11;
    11 -0.6770370970 1.0633892659 -1.2391310131 1.6521746842 12;
    12 -0.8837094564 0.6492546703 -0.6118929981 1.4533253571 13;]

# Gegenbauer
function poch(x,n)
    if n == 0
        return 1
    end
    res = x
    for i in 1:n-1
        res *= x+i
    end
    return res
end
TriangleCoeffsGegenbauer = []
Cα0(x, α) = 1
    Cα1(x, α) = 2α*x
    Cα2(x, α) = -α + 2*poch(α, 2)*x^2
    Cα3(x, α) = (1/3) * (-6poch(α, 2)x + 4poch(α, 3)x^3)
    Cα4(x, α) = (1/6) * (3poch(α, 2) - 12poch(α, 3)x^2 + 4poch(α, 4)x^4)
    Cα5(x, α) = (1/15) * (15poch(α, 3)x - 20poch(α, 4)x^3 + 4poch(α, 5)x^5)
    Cα6(x, α) = (1/90) * (-15poch(α, 3) + 90poch(α, 4)x^2 - 60poch(α, 5)x^4 + 8poch(α, 6)x^6)

    Cαs = [Cα0 Cα1 Cα2 Cα3 Cα4 Cα5 Cα6]

    # Abramowitz and Stegun, example 3, Table 22.2, p. 789
    # For C(.5,n+1)(2.5) == [5(n+.25)*C(.5,n)(2.5) - (n-.5)*C(.5, n-1)(2.5)]*(1/(n+1))

# Obtained from ApproxFun with Ultraspherical(λ == 2)
# using ApproxFun
# cs = [Fun(Ultraspherical(2), [zeros(i);1]) for i in 0:12]
# [cs[i](j) for i in 1:7, j in -1:.125:1]
const TestTableGegenbauer = [
      1.0    1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0   1.0       1.0       1.0       1.0;
     -4.0   -3.5      -3.0      -2.5      -2.0  -1.5      -1.0      -0.5       0.0   0.5       1.0       1.5       2.0   2.5       3.0       3.5       4.0;
     10.0    7.1875    4.75      2.6875    1.0  -0.3125   -1.25     -1.8125   -2.0  -1.8125   -1.25     -0.3125    1.0   2.6875    4.75      7.1875   10.0;
    -20.0  -10.9375   -4.5      -0.3125    2.0   2.8125    2.5       1.4375    0.0  -1.4375   -2.5      -2.8125   -2.0   0.3125    4.5      10.9375   20.0;
     35.0   13.1445    1.3125   -3.54297  -4.0  -2.16797   0.3125    2.26953   3.0   2.26953   0.3125   -2.16797  -4.0  -3.54297   1.3125   13.1445   35.0;
    -56.0  -12.291     3.9375    5.75195   2.0  -1.98633  -3.6875   -2.69336   0.0   2.69336   3.6875    1.98633  -2.0  -5.75195  -3.9375   12.291    56.0;
     84.0    7.56812  -8.64063  -3.66431   3.0   4.62866   1.73438  -2.24048  -4.0  -2.24048   1.73438   4.62866   3.0  -3.66431  -8.64063   7.56812  84.0;]

# Jacobi
# Table 22.1 of Abramowitz and Stegun, p. 793
Pαβ0(α , β, x) = 1
    Pαβ1(α , β, x) = .5 * (2*(α + 1) + (α + β + 2)*(x - 1))
    Pαβ2(α , β, x) = .125 * (4*poch(α + 1, 2) + 4*(α + β + 3)*poch(α + 2, 1)*(x - 1) + poch(α + β + 3, 2)*(x - 1)^2)
    Pαβ3(α , β, x) = (1/48) * (8*poch(α + 1, 3) + 12*(α + β + 4)*poch(α + 2, 2)*(x - 1) + 6*poch(α + β + 4, 2)*poch(α + 3, 1)*(x - 1)^2 + poch(α + β + 4, 3)*(x - 1)^3)
    Pαβ4(α , β, x) = (1/384) * (16*poch(α + 1, 4) + 32*(α + β + 5)*poch(α + 2, 3)*(x - 1) + 24*poch(α + β + 5, 2)*poch(α + 3, 2)*(x - 1)^2 + 8*poch(α + β + 5, 3)*poch(α + 4, 1)*(x - 1)^3 + poch(α + β + 5,4)*(x - 1)^4)
    Pαβ5(α , β, x) = (1/3840) * (32*poch(α + 1, 5) + 80*(α + β + 6)*poch(α + 2, 4)*(x - 1) + 80*poch(α + β + 6, 2)*poch(α + 3, 3)*(x - 1)^2 + 40*poch(α + β + 6, 3)*poch(α + 4, 2)*(x - 1)^3 + 10*poch(α + β + 6, 4)*poch(α + 5, 1)*(x - 1)^4 + poch(α + β + 6, 5)*(x-1)^5)
    Pαβ6(α , β, x) = (1/46080) * (64*poch(α + 1, 6) + 192*(α + β + 7)*poch(α + 2, 5)*(x - 1) + 240*poch(α + β + 7, 2)*poch(α + 3, 4)*(x - 1)^2 + 160*poch(α + β + 7, 3)*poch(α + 4, 3)*(x - 1)^3 + 60*poch(α + β + 7, 4)*poch(α + 5, 2)*(x - 1)^4 + 12*poch(α + β + 7, 5)*(α + 6)*(x - 1)^5 + poch(α + β + 7, 6)*(x-1)^6)

    Pαβs = [Pαβ0 Pαβ1 Pαβ2 Pαβ3 Pαβ4 Pαβ5 Pαβ6]

# Obtained from ApproxFun.jl for Jacobi(2 == β, 3 == α)
# Note that we do α THEN β whereas they do it the other way around.
# using ApproxFun
# ps = [Fun(Jacobi(2,3), [zeros(i); 1]) for i in 0:12]
# [ps[i](j) for i in 1:7, j in -1:.125:1]
const TestTableJacobi = [
      1.0   1.0       1.0        1.0         1.0        1.0        1.0         1.0        1.0       1.0        1.0        1.0        1.0        1.0        1.0        1.0       1.0;
     -3.0  -2.5625   -2.125     -1.6875     -1.25      -0.8125    -0.375       0.0625     0.5       0.9375     1.375      1.8125     2.25       2.6875     3.125      3.5625    4.0;
      6.0   4.14063   2.5625     1.26563     0.25      -0.484375  -0.9375     -1.10938   -1.0      -0.609375   0.0625     1.01563    2.25       3.76563    5.5625     7.64063  10.0;
    -10.0  -5.21362  -1.94336    0.0524902   1.01563    1.18774    0.810547    0.125732  -0.625    -1.19995   -1.35742   -0.855713   0.546875   3.09204    7.02148   12.5769   20.0;
     15.0   5.4113    0.487061  -1.31355    -1.26953   -0.398026   0.545654    1.06804    0.9375    0.18425   -0.899658  -1.76033   -1.58203    0.712814   6.46362   17.2717   35.0;
    -21.0  -4.59868   1.14267    1.60243     0.280273  -0.859185  -1.04269    -0.344714   0.65625   1.22132    0.794891  -0.651664  -2.27637   -2.02295    3.72281   20.4768   56.0;
     28.0   2.90217  -2.18814   -0.735861    0.943359   1.01195   -0.0296936  -0.952181  -0.875     0.18675    1.30673    1.16287   -0.916016  -3.3725    -0.311676  21.1388   84.0]


@testset "OP Table Cases" begin
    @testset "Laguerre" begin
        @testset "Hand Implemenation" begin
            @test all(f(0) ≈ 1.0 for f in Ls)
            @test all(Ls[i](testvalslaguerre[j]) ≈ TestTableLaguerre[i,j] for i in 1:13, j in 2:6)
            @test all([isapprox.(Ls[i](testvalslaguerre[j]),TestTableLaguerre[i,j]) for i in 1:13, j in 2:6])
        end
        @testset "Naive Implementation" begin
            @test all(laguerre(i,0.0) ≈ 1.0 for i in 0:12)
            @test all([laguerre(i-1,testvalslaguerre[j]) ≈ TestTableLaguerre[i,j] for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            α = 0
            @test_broken all(eval.([a(0,i) for i in 1:12])  .≈ 1)
            @test_broken all(eval.([a(testvalslaguerre[j],i) for i in 0:12, j in 2:6]) .≈ TestTableLaguerre[:,2:6])
        end
    end

    @testset "Hermite" begin
        @testset "Hand Implemented" begin
            @test all(Hs[i](testvalshermite[j]) ≈ TestTableHermite[i,j] for i in 1:13, j in 2:6)
        end
        @testset "Naive Implementation" begin
            @test all(hermite(i,0.0) ≈ 1.0 for i in 1:12)
            @test all([hermite(testvalshermite[j], i-1) ≈ TestTableHermite[i,j] for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Legendre" begin
        @testset "Hand Implemented" begin
            @test all([isapprox.(testPs[j-2](TestTableLegendre[i,1]),TestTableLegendre[i,j]; atol= 6.9e-8) for i in 1:101, j in 3:6]) # Do not the godawful atol keyword argument to shimmy all tests into passing.
        end
        @testset "Naive Recursion" begin
            @test all([legendre(i,0) for i in 1:10] .≈ 1)
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 2), TestTableLegendre[i,3]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 3), TestTableLegendre[i,4]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 9), TestTableLegendre[i,5]; atol=5.e-6) for i in 1:101])
            @test all([isapprox.(legendre(TestTableLegendre[i,1], 10), TestTableLegendre[i,6]; atol=5.e-6) for i in 1:101])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 ==2
        end
    end

    @testset "Chebyshev T First Kind" begin
        @testset "Hand Implemented" begin
            @test all(Ts[i](testvalschebyshev[j]) ≈ TestTableChebyshevT[i,j+1] for i in 1:13, j in 1:5)
        end
        @testset "Naive Recursion" begin
            @test all([isapprox.(chebyshevfirst(testvalschebyshev[j-1], i-1), TestTableChebyshevT[i,j]) for i in 1:13, j in 2:6])
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Chebyshev U Second Kind" begin
        @testset "Hand Implemented" begin
            @test all(Us[i](testvalschebyshev[j-1]) ≈ TestTableChebyshevU[i,j] for i in 1:13, j in 2:6)
        end
        @testset "Naive Recursion" begin
            @test_broken 1 == 2
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Gegenbauer/Ultraspherical" begin
        @testset "Hand Implementation" begin
            testcs() = [Cαs[i](j, 2) for i in 1:7, j in -1:.125:1]
            @test isapprox.(testcs(), TestTableGegenbauer; atol = 4e-5) |> all
        end
        # Abramowitz and Stegun, example 3, Table 22.2, p. 789
        @testset "Naive Recursion" begin
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end

    @testset "Jacobi" begin
        @testset "Hand Implementation" begin
            test() = [Pαβs[i](3,2,j) for i in 1:7, j in -1:.125:1]
            @test isapprox.(test(), TestTableJacobi; atol = 5e-5) |> all
        end
        @testset "Naive Recursion" begin
            @test_broken 1==2
        end
        @testset "Stubborn Way" begin
            @test_broken 1 == 2
        end
    end


        j = Jacobi(.5,1.5,8)
            ge = Gegenbauer(3,4)
            go = Gegenbauer(4,5)
            ue = ChebyshevSecondKind(4)
            uo = ChebyshevSecondKind(5)
            te = ChebyshevFirstKind(4)
            to = ChebyshevFirstKind(5)
            le = Legendre(4)
            lo = Legendre(5)
            he = Hermite(4)
            ho = Hermite(5)
            lae = Laguerre(2,4)
            lao = Laguerre(3,5)

    @testset "Constructors and Parity" begin
        @test j == Jacobi{.5,1.5,8}()
        @test ge == GegenbauerEven{3,2}()
        @test go == GegenbauerOdd{4,2}()
        @test ue == ChebyshevSecondKindEven{2}()
        @test uo == ChebyshevSecondKindOdd{2}()
        @test te == ChebyshevFirstKindEven{2}()
        @test to == ChebyshevFirstKindOdd{2}()
        @test le == Legendre{4}()
        @test lo == Legendre{5}()
        @test he == HermiteEven{2}()
        @test ho == HermiteOdd{2}()
        @test lae == Laguerre{2,4}()
        @test lao == Laguerre{3,5}()
    end

    # Abramowitz pg. 790, Example 2.
    @testset "Jacobi{.5, 1.5, 8}(2) for a, b, c, d, f" begin
        @test_broken all([a(j,2,i) for i in 8:-1:0] .≈ [1, 1.132353, 1.366667, 1.841026, 3.008392, 6.849651, 26.44156, 223.1091, 6545.533])
        @test all([b(j,i) for i in 8:-1:0] .== [18.0, 34.0, 48.0, 60.0, 70.0, 78.0, 84.0, 88.0, 90.0])
        @test all([c(j,i) for i in 8:-1:0] .== [136.0, 105.0, 78.0, 55.0, 36.0, 21.0, 10.0, 3.0, 0.0])
        @test d(j,2) ≈ 3.338470458984375
        @test f(j,2) == -1
        @test a(j,2) ≈ 21852.07 atol=.01
    end

    # Abramowitz pg. 790, Example 3
    @testset "Gegenbauer{.5, 0:6}(2.5)" begin
        g0 = Gegenbauer(.5, 0)
        g1 = Gegenbauer(.5, 1)
        g2 = Gegenbauer(.5, 2)
        g3 = Gegenbauer(.5, 3)
        g4 = Gegenbauer(.5, 4)
        g5 = Gegenbauer(.5, 5)
        g6 = Gegenbauer(.5, 6)
        @test isapprox(a(g0,2.5), 1; atol=5e-6)
        @test_broken isapprox(a(g1,2.5), 1.25; atol=5e-6)
        @test_broken isapprox(a(g2,2.5) , 3.65625; atol=5e-6)
        @test_broken isapprox(a(g3,2.5) , 13.08594; atol=5e-6)
        @test_broken isapprox(a(g4,2.5) , 50.87648; atol=5e-6)
        @test_broken isapprox(a(g5,2.5) , 207.0649; atol=5e-6)
        @test_broken isapprox(a(g6,2.5), 867.7516; atol=5e-6)

    end

end
