# precise values were computed with Maple 2018.2.1
@testset "Smallest parameter is equal to 1" begin
  a = Float64[1, 10, 100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    @test SimpleBirthDeathProcess.log_meixner_ortho_poly(i, one(Float64), log(x)) ≈ -log(z[u])
  end
end

@testset "Smallest parameter is equal to 10" begin
  a = Float64[10, 100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F = [-.52628324337533965e-3 -.91738185249562661e-3 -.99103028876012253e-3 -.99905081268556093e-3;
       -2.0317801262571448 -3.6572455938562417 -4.0102491632437109 -4.0501571465136512;
       -3.3504104733399496 -6.1641594451500663 -6.8430557799176588 -6.9224879687084634;
       -6.9013533188221590 -13.505607286760857 -15.749584728918338 -16.058540058411344;
       -8.9931244812245734 -18.252192303278967 -22.280051500153905 -22.945569932182599;
       -13.364668417083967 -28.631101920969909 -40.723125251220806 -45.230948827869905]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    @test SimpleBirthDeathProcess.log_meixner_ortho_poly(i, Float64(10), log(x)) ≈ F[u, v]
  end
end

@testset "Smallest parameter is equal to 100" begin
  a = Float64[100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F = [-.50248118964360310e-2 -.90986888595390629e-2 -.99014706210901467e-2;
       -19.349403266385754 -36.223694926634441 -40.058734895006031;
       -31.843688770301594 -60.988773615097764 -68.343834334910785;
       -65.106469020297578 -133.11309240174212 -157.16644812648946;
       -84.255629726916470 -179.31620948692071 -222.11471731657940;
       -120.73661754514218 -277.84032552047403 -403.55806307301357]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    @test SimpleBirthDeathProcess.log_meixner_ortho_poly(i, Float64(100), log(x)) ≈ F[u, v]
  end
end

@testset "Smallest parameter is equal to 1_000" begin
  a = Float64[1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F = [-.50021886485016083e-1 -.90912435206967968e-1;
       -192.57640014113831 -361.89231624299969;
       -316.86775046147991 -609.24410302537978;
       -647.41688830508753 -1329.2305629085790;
       -837.33067117580264 -1790.0472246195170;
       -1196.8270424842711 -2770.4365897006691]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    @test SimpleBirthDeathProcess.log_meixner_ortho_poly(i, Float64(1_000), log(x)) ≈ F[u, v]
  end
end

@testset "k in [1, 0, -1]: Smallest parameter is equal to 2" begin
  a = Float64[2, 10, 100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    r1, r2, r3 = SimpleBirthDeathProcess.log_meixner_ortho_poly_joint(i, Float64(2), log(x))
    @test r1 ≈ log(i + 2 * z[u] - 1) - 2 * log(z[u]) - log1p(i)
    @test r2 ≈ log(i + z[u] - 1) - log(z[u]) - log(i)
    @test r3 === zero(Float64)
  end
end

@testset "k in [1, 0, -1]: Smallest parameter is equal to 10" begin
  a = Float64[10, 100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F1 = [-.52628324337533965e-3 -.91738185249562661e-3 -.99103028876012253e-3 -.99905081268556093e-3;
        -2.0317801262571448 -3.6572455938562417 -4.0102491632437109 -4.0501571465136512;
        -3.3504104733399496 -6.1641594451500663 -6.8430557799176588 -6.9224879687084634;
        -6.9013533188221590 -13.505607286760857 -15.749584728918338 -16.058540058411344;
        -8.9931244812245734 -18.252192303278967 -22.280051500153905 -22.945569932182599;
        -13.364668417083967 -28.631101920969909 -40.723125251220806 -45.230948827869905]

  F2 = [-.44997154621306883e-3 -.82495557250215841e-3 -.89191929546486659e-3 -.89914565055050351e-3;
        -1.7268423647221549 -3.2876137827946026 -3.6091750283602051 -3.6451409270498669;
        -2.8344139388261806 -5.5392105465265657 -6.1586334908790476 -6.2302379617871075;
        -5.7391285605322099 -12.112750882906229 -14.173760567088424 -14.452676449894550;
        -7.3620201029852153 -16.326937625328312 -20.048881722140631 -20.650973873851784;
        -10.125594848132590 -25.159927670658581 -36.558405822355922 -40.704713967442946]

  F3 = [-.37644616098559351e-3 -.73267076952353524e-3 -.79281005797932914e-3 -.79924050637141279e-3;
        -1.4345092663145545 -2.9187743959822365 -3.2081117219182366 -3.2401248196531340;
        -2.3417773974599485 -4.9159716647216521 -5.4742368237939885 -5.5379882234321887;
        -4.6484173208971883 -10.728016496952154 -12.598124000963975 -12.846814969489424;
        -5.8601107881637621 -14.420485859284776 -17.818384875997870 -18.356386451070455;
        -7.6035926506422622 -21.822978690338372 -32.411557448038426 -36.179147716752838]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    r1, r2, r3 = SimpleBirthDeathProcess.log_meixner_ortho_poly_joint(i, Float64(10), log(x))
    @test r1 ≈ F1[u, v]
    @test r2 ≈ F2[u, v]
    @test r3 ≈ F3[u, v]
  end
end

@testset "k in [1, 0, -1]: Smallest parameter is equal to 100" begin
  a = Float64[100, 1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F1 = [-.50248118964360310e-2 -.90986888595390629e-2 -.99014706210901467e-2;
        -19.349403266385754 -36.223694926634441 -40.058734895006031;
        -31.843688770301594 -60.988773615097764 -68.343834334910785;
        -65.106469020297578 -133.11309240174212 -157.16644812648946;
        -84.255629726916470 -179.31620948692071 -222.11471731657940;
        -120.73661754514218 -277.84032552047403 -403.55806307301357]

  F2 = [-.49496903336308133e-2 -.90068897338010612e-2 -.98024463036859674e-2;
        -19.049868084336226 -35.856894960464313 -39.658088245693884;
        -31.337725781208335 -60.369013437295331 -67.660255610319402;
        -63.975125304549095 -131.73475269960223 -155.59375275504627;
        -82.680492341128495 -177.41449504455335 -219.88985501919318;
        -117.81301911039979 -274.45170026360770 -399.42270183276153]

  F3 = [-.48748212962385064e-2 -.89151055509698546e-2 -.97034221785626762e-2;
        -18.751468651711550 -35.490177672488245 -39.257442779810132;
        -30.833857028198276 -59.749429848752658 -66.976679679969060;
        -62.850044157717020 -130.35722132504727 -154.02107758909523;
        -81.116336230254627 -175.51460195041096 -217.66506396756262;
        -114.93595962960294 -271.07491295298205 -395.28907222661447]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    r1, r2, r3 = SimpleBirthDeathProcess.log_meixner_ortho_poly_joint(i, Float64(100), log(x))
    @test r1 ≈ F1[u, v]
    @test r2 ≈ F2[u, v]
    @test r3 ≈ F3[u, v]
  end
end

@testset "k in [1, 0, -1]: Smallest parameter is equal to 1_000" begin
  a = Float64[1_000, 10_000]
  z = Float64[1.0001, 1.5, 2.0, 5.0, 10.0, 100.0]

  F1 = [-.50021886485016083e-1 -.90912435206967968e-1;
        -192.57640014113831 -361.89231624299969;
        -316.86775046147991 -609.24410302537978;
        -647.41688830508753 -1329.2305629085790;
        -837.33067117580264 -1790.0472246195170;
        -1196.8270424842711 -2770.4365897006691]

  F2 = [-.49946878041294134e-1 -.90820697729011826e-1;
        -192.27737387819252 -361.52579376188207;
        -316.36272606764623 -608.62485012554081;
        -646.28835292904761 -1327.8536325701771;
        -835.76046016206889 -1788.1477816590733;
        -1193.9243639196877 -2767.0554727792707]

  F3 = [-.49871894621978197e-1 -.90728961752889464e-1;
        -191.97846010652304 -361.15927958063496;
        -315.85790900945579 -608.00561493593598;
        -645.16043639731818 -1326.4767830232228;
        -834.19133195599179 -1786.2485202922418;
        -1191.0262000516550 -2763.6755269221370]

  @testset for (v, i) = enumerate(a), (u, x) = enumerate(z)
    r1, r2, r3 = SimpleBirthDeathProcess.log_meixner_ortho_poly_joint(i, Float64(1_000), log(x))
    @test r1 ≈ F1[u, v]
    @test r2 ≈ F2[u, v]
    @test r3 ≈ F3[u, v]
  end
end
