Optimize polynomial evaluation of coefficients in Gauss-Laguerre (#112)

* use evalpoly for gl_bessel * add gl_bulk * fix bad comment * remove @evalpoly and update coefs * Remove explicit 1.8 tests * optimize airy region coefs Co-authored-by: Sheehan Olver <solver@mac.com>
JuliaApproximation · Oct 31, 2022 · 03b6fb5 · 03b6fb5
1 parent 15cf476
commit 03b6fb5
Show file tree

Hide file tree

Showing 4 changed files with 126 additions and 186 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -17,7 +17,6 @@ jobs:
         version:
           - '1.6'
           - '1'
-          - '^1.8.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

diff --git a/src/besselroots.jl b/src/besselroots.jl
@@ -156,12 +156,12 @@ function besselZeroRoots(m)
     @inbounds for jj = 21:min(m, 47)
         ak = π * (jj - .25)
         ak82 = (.125 / ak)^2
-        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5], p[3])
+        jk[jj] = ak + .125 / ak * evalpoly(ak82, (1.0, p[7], p[5], p[3]))
     end
     @inbounds for jj = 48:min(m, 344)
         ak = π * (jj - .25)
         ak82 = (.125 / ak)^2
-        jk[jj] = ak + .125 / ak * @evalpoly(ak82, 1.0, p[7], p[5])
+        jk[jj] = ak + .125 / ak * evalpoly(ak82, (1.0, p[7], p[5]))
     end
     @inbounds for jj = 345:min(m,13191)
         ak = π * (jj - .25)
@@ -189,20 +189,17 @@ function besselJ1(m)
     @inbounds for jj = 11:min(m, 15)
         ak = π * (jj - .25)
         ak2 = (1 / ak)^2
-        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3],
-                                                  c[2], c[1]), ak2^2, c[7])
+        Jk2[jj] = 1 / (π * ak) * muladd(evalpoly(ak2, (c[5], c[4], c[3], c[2], c[1])), ak2^2, c[7])
     end
     @inbounds for jj = 16:min(m, 21)
         ak = π * (jj - .25)
         ak2 = (1 / ak)^2
-        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3], c[2]),
-                                        ak2^2, c[7])
+        Jk2[jj] = 1 / (π * ak) * muladd(evalpoly(ak2, (c[5], c[4], c[3], c[2])), ak2^2, c[7])
     end
     @inbounds for jj = 22:min(m,55)
         ak = π * (jj - .25)
         ak2 = (1 / ak)^2
-        Jk2[jj] = 1 / (π * ak) * muladd(@evalpoly(ak2, c[5], c[4], c[3]),
-                                        ak2^2, c[7])
+        Jk2[jj] = 1 / (π * ak) * muladd(evalpoly(ak2, (c[5], c[4], c[3])), ak2^2, c[7])
     end
     @inbounds for jj = 56:min(m,279)
         ak = π * (jj - .25)

diff --git a/src/gausslaguerre.jl b/src/gausslaguerre.jl
@@ -267,107 +267,116 @@ end
 ## The bulk
 
 # Note: there is always one division by an integer, placed such that it preserves the type of `d`
-gl_bulk_x3(t, d, α) = -(12*α^2 + 5*(1-t)^(-2) - 4*(1-t)^(-1) - 4) * d / 12
-gl_bulk_x5(t, d, α) = d^3*(1-t)/t/720*(1600*(1-t)^(-6) - 3815*(1-t)^(-5)
-    + 480*α^4 +2814*(1-t)^(-4) - 576*(1-t)^(-3) - 960*α^2 - 48*(15*α^4 - 30*α^2 + 7)*(1-t)^(-1)
-    -16*(1-t)^(-2) + 224)
-gl_bulk_x7(t, d, α) = -d^5/181440*(1-t)^2/t^2*(10797500*(1-t)^(-10) - 43122800*(1-t)^(-9)
-    + 66424575*(1-t)^(-8) -48469876*(1-t)^(-7) + 193536*α^6 + 16131880*(1-t)^(-6)
-    + 80*(315*α^4 - 630*α^2 -221)*(1-t)^(-4) - 1727136*(1-t)^(-5)
-    - 967680*α^4 - 320*(63*α^4 - 126*α^2 +43)*(1-t)^(-3)
-    + 384*(945*α^6 - 4620*α^4 + 6405*α^2 - 1346)*(1-t)^(-2)
-    + 1354752*α^2   - 23040*(21*α^6 - 105*α^4 + 147*α^2 - 31)*(1-t)^(-1) -285696)
-gl_bulk_x9(t, d, α) = d^7/10886400*(1-t)^3/t^3*(43222750000*(1-t)^(-14) - 241928673000*(1-t)^(-13)
-    + 566519158800*(1-t)^(-12) -714465642135*(1-t)^(-11) + 518401904799*(1-t)^(-10)
-    + 672*(12000*α^4 - 24000*α^2 +64957561)*(1-t)^(-8) - 212307298152*(1-t)^(-9)
-    + 24883200*α^8 - 192*(103425*α^4 -206850*α^2 + 15948182)*(1-t)^(-7)
-    + 3360*(4521*α^4 - 9042*α^2 - 7823)*(1-t)^(-6) -232243200*α^6
-    - 1792*(3375*α^6 - 13905*α^4  + 17685*α^2 - 1598)*(1-t)^(-5)
-    + 16128*(450*α^6 - 2155*α^4 + 2960*α^2 - 641)*(1-t)^(-4)
-    + 812851200*α^4 -768*(70875*α^8 - 631260*α^6 + 2163630*α^4
-    - 2716980*α^2 +555239)*(1-t)^(-3)  + 768*(143325*α^8 - 1324260*α^6
-    + 4613070*α^4 -5826660*α^2 + 1193053)*(1-t)^(-2) - 1028505600*α^2
-    - 5806080*(15*α^8 -140*α^6 + 490*α^4 - 620*α^2 + 127)*(1-t)^(-1) + 210677760)
-
-gl_bulk_w3(t, d, α) = d^2/6*(2*t + 3)/(t-1)^3
-gl_bulk_w5(t, d, α) = (1-t)^2/720/t^2*d^4*(8000*(1-t)^(-8) - 24860*(1-t)^(-7) + 27517*(1-t)^(-6)
-    - 12408*(1-t)^(-5) + 1712*(1-t)^(-4) +16*(15*α^4 - 30*α^2 + 7)*(1-t)^(-2) + 32*(1-t)^(-3))
-gl_bulk_w7(t, d, α) = -(1-t)^3/90720/t^3*d^6*(43190000*(1-t)^(-12) -204917300*(1-t)^(-11)
-    + 393326325*(1-t)^(-10) - 386872990*(1-t)^(-9) + 201908326*(1-t)^(-8)
-    + 80*(315*α^4 - 630*α^2 + 53752)*(1-t)^(-6)  - 50986344*(1-t)^(-7)
-    - 320*(189*α^4 -378*α^2 - 89)*(1-t)^(-5) + 480*(63*α^4 - 126*α^2
-    + 43)*(1-t)^(-4)  -384*(315*α^6 - 1470*α^4 + 1995*α^2 - 416)*(1-t)^(-3)
-    + 2304*(21*α^6 -105*α^4 + 147*α^2 - 31)*(1-t)^(-2) )
-
-# And for α = 0
-gl_bulk_x3(t, d) = -d/12*(5*(1-t)^(-2) - 4*(1-t)^(-1) - 4)
-gl_bulk_x5(t, d) = d^3*(1-t)/t/720*(1600*(1-t)^(-6) - 3815*(1-t)^(-5) + 2814*(1-t)^(-4)
-    - 576*(1-t)^(-3) - 48*7*(1-t)^(-1) -16*(1-t)^(-2) + 224)
-gl_bulk_x7(t, d) = -d^5/181440*(1-t)^2/t^2*(10797500*(1-t)^(-10)
-    - 43122800*(1-t)^(-9) + 66424575*(1-t)^(-8) -48469876*(1-t)^(-7)
-    + 16131880*(1-t)^(-6) - 80*221*(1-t)^(-4) - 1727136*(1-t)^(-5) - 320*43*(1-t)^(-3)
-    - 384*1346*(1-t)^(-2) + 23040*31*(1-t)^(-1) -285696)
-gl_bulk_x9(t, d) = d^7/10886400*(1-t)^3/t^3*(43222750000*(1-t)^(-14)
-    - 241928673000*(1-t)^(-13) + 566519158800*(1-t)^(-12) -714465642135*(1-t)^(-11)
-    + 518401904799*(1-t)^(-10) + 672*64957561*(1-t)^(-8)   - 212307298152*(1-t)^(-9)
-    - 192*15948182*(1-t)^(-7)  - 3360*7823*(1-t)^(-6) + 1792*1598*(1-t)^(-5)
-    + 16128*(- 641)*(1-t)^(-4)  -768*555239*(1-t)^(-3)  + 768*1193053*(1-t)^(-2)
-    - 5806080*127*(1-t)^(-1) + 210677760)
-
-gl_bulk_w3(t, d) = d^2/6*(2*t + 3)/(t-1)^3
-gl_bulk_w5(t, d) = (1-t)^2/720/t^2*d^4*(8000*(1-t)^(-8) - 24860*(1-t)^(-7)
-    + 27517*(1-t)^(-6) - 12408*(1-t)^(-5) + 1712*(1-t)^(-4) +16*7*(1-t)^(-2)
-    + 32*(1-t)^(-3))
-gl_bulk_w7(t, d) = -(1-t)^3/90720/t^3*d^6*(43190000*(1-t)^(-12)
-    - 204917300*(1-t)^(-11) + 393326325*(1-t)^(-10) - 386872990*(1-t)^(-9)
-    + 201908326*(1-t)^(-8) +80*53752*(1-t)^(-6)
-    - 50986344*(1-t)^(-7) + 320*89*(1-t)^(-5)
-    + 480*43*(1-t)^(-4) + 384*416*(1-t)^(-3) - 2304*31*(1-t)^(-2) )
+function gl_bulk(t, d, α)
+    α² = α * α
+    tinv = inv(1-t)
+    _t = (1-t) / t
+    _t² = _t * _t
+    _t³ = _t² * _t
+    d² = d * d
+    d⁴ = d² * d²
+    d⁶ = d⁴ * d²
+
+    c1 = evalpoly(α², (-4, 12))
+    x3 = -evalpoly(tinv, (c1, -4, 5)) * d / 12
+    c1, c2 = evalpoly(α², (224, -960, 480)), evalpoly(α², (7, -30, 15))
+    x5 = d² * d * _t / 720 * evalpoly(tinv, (c1, -48*c2, -16, -576, 2814, -3815, 1600))
+    c1, c2, c3 = evalpoly(α², (-285696, 1354752, -967680, 193536)), evalpoly(α², (-31, 147, -105, 21)), evalpoly(α², (-1346, 6405, -4620, 945))
+    c4, c5 = evalpoly(α², (43, -126, 63)), evalpoly(α², (-221, -630, 315))
+    x7 = -d⁴ * d / 181440 * _t² * evalpoly(tinv, (c1, -23040*c2, 384*c3, -320*c4, 80*c5, -1727136, 16131880, -48469876, 66424575, -43122800, 10797500))
+    c1, c2, c3 = evalpoly(α², (210677760, -1028505600, 812851200, -232243200, 24883200)), evalpoly(α², (127, -620, 490, -140, 15)), evalpoly(α², (1193053, -5826660, 4613070, -1324260, 143325))
+    c4, c5, c6 = evalpoly(α², (555239, -2716980, 2163630, -631260, 70875)), evalpoly(α², (-641, 2960, -2155, 450)), evalpoly(α², (-1598, 17685, -13905, 3375))
+    c7, c8, c9 = evalpoly(α², (-7823, -9042, 4521)), evalpoly(α², (15948182, -206850, 103425)), evalpoly(α², (64957561, -24000, 12000))
+    x9 = d⁶ * d / 10886400 * _t³ * evalpoly(tinv, (c1, -5806080*c2, 768*c3, -768*c4, 16128*c5, -1792*c6, 3360*c7, -192*c8, 672*c9, -212307298152, 518401904799, -714465642135, 566519158800, -241928673000, 43222750000))
+
+    w3 = d² / 6 * (2*t + 3) / (t-1)^3
+    c1 = evalpoly(α², (7, -30, 15))
+    w5 = d⁴ * _t² / 720  * evalpoly(tinv, (0, 0, 16*c1, 32, 1712, -12408, 27517, -24860, 8000))
+    c1, c2, c3 = evalpoly(α², (-31, 147, -105, 21)), evalpoly(α², (-416, 1995, -1470, 315)), evalpoly(α², (43, -126, 63))
+    c4, c5 = evalpoly(α², (-89, -378, 189)), evalpoly(α², (53752, -630, 315))
+    w7 = -d⁶ * _t³ / 90720 * evalpoly(tinv, (0, 0, 2304*c1, -384*c2, 480*c3, -320*c4, 80*c5, -50986344, 201908326, -386872990, 393326325, -204917300, 43190000))
+    return (x3, x5, x7, x9), (w3, w5, w7)
+end
 
+function gl_bulk(t, d)
+    tinv = inv(1-t)
+    _t = (1-t) / t
+    _t² = _t * _t
+    _t³ = _t² * _t
+    d² = d * d
+    d⁴ = d² * d²
+    d⁶ = d⁴ * d²
+
+    x3 = -d / 12 * evalpoly(tinv, (-4, -4, 5))
+    x5 = d² * d * _t / 720 * evalpoly(tinv, (224, -336, -16, -576, 2814, -3815, 1600))
+    x7 = -d⁴ * d / 181440 * _t² * evalpoly(tinv, (-285696, 714240, -516864, -13760, -17680, -1727136, 16131880, -48469876, 66424575, -43122800, 10797500))
+    x9 = d⁶ * d / 10886400 * _t³ * evalpoly(tinv, (210677760, -737372160, 916264704, -426423552, -10338048, 2863616, -26285280, -3062050944, 43651480992, -212307298152, 518401904799, -714465642135, 566519158800, -241928673000, 43222750000))
+
+    w3 = d² / 6 * (2*t + 3) / (t-1)^3
+    w5 = d⁴ / 720 * _t² * evalpoly(tinv, (0, 0, 112, 32, 1712, -12408, 27517, -24860, 8000))
+    w7 = -d⁶ * _t³ / 90720 * evalpoly(tinv, (0, 0, -71424, 159744, 20640, 28480, 4300160, -50986344, 201908326, -386872990, 393326325, -204917300, 43190000))
+    return (x3, x5, x7, x9), (w3, w5, w7)
+end
 
 ## The hard edge (Bessel region)
-
-gl_bessel_x3(jak, d, α) = (jak^2 + 2*α^2 - 2)*d^2 / 3
-gl_bessel_x5(jak, d, α) = (11*jak^4 +3*jak^2*(11*α^2-19) +46*α^4 -140*α^2 +94)*d^4 / 45
-gl_bessel_x7(jak, d, α) = (657*jak^6 +36*jak^4*(73*α^2-181) +2*jak^2*(2459*α^4 -10750*α^2 +14051)
-    + 4*(1493*α^6 -9303*α^4 +19887*α^2 - 12077) )*d^6 / 2835
-gl_bessel_x9(jak, d, α) = (10644*jak^8 + 60*(887*α^2 - 2879)*jak^6 + (125671*α^4 -729422*α^2 + 1456807)*jak^4
-    + 3*(63299*α^6 - 507801*α^4 + 1678761*α^2 - 2201939)*jak^2 + 2*(107959*α^8
-    - 1146220*α^6 + 5095482*α^4 -10087180*α^2 + 6029959) )*d^8 / 42525
-
-gl_bessel_w3(jak, d, α) = (α^2 + jak^2 -1)*2*d^2 / 3
-gl_bessel_w5(jak, d, α) = (46*α^4 + 33*jak^4 +6*jak^2*(11*α^2 -19) -140*α^2 +94)*d^4 / 45
-gl_bessel_w7(jak, d, α) = (1493*α^6 + 657*jak^6 + 27*(73*α^2 - 181)*jak^4 - 9303*α^4
-    + (2459*α^4 -10750*α^2 + 14051)*jak^2 + 19887*α^2 - 12077)*4*d^6 / 2835
-gl_bessel_w9(jak, d, α) = (215918*α^8 + 53220*jak^8 + 240*(887*α^2 - 2879)*jak^6 -2292440*α^6 +
-    3*(125671*α^4 - 729422*α^2 + 1456807)*jak^4 + 10190964*α^4  +
-    6*(63299*α^6 - 507801*α^4 + 1678761*α^2 -2201939)*jak^2 -
-    20174360*α^2 + 12059918)*d^8 / 42525
+function gl_bessel(jak, d, α)
+    jak² = jak * jak
+    d² = d * d
+    d⁴ = d² * d²
+    d⁸ = d⁴ * d⁴
+    α² = α * α
+
+    c1 = evalpoly(α², (-2, 2))
+    x3 = d² * evalpoly(jak², (c1, 1)) / 3
+    w3 = d² * evalpoly(jak², (c1, 2)) / 3
+    c1, c2 = evalpoly(α², (94, -140, 46)), evalpoly(α², (-19, 11)) 
+    x5 = d⁴ * evalpoly(jak², (c1, 3*c2, 11)) / 45
+    w5 = d⁴ * evalpoly(jak², (c1, 6*c2, 33)) / 45
+    c1, c2, c3 = evalpoly(α², (-12077, 19887, -9303, 1493)), evalpoly(α², (14051, -10750, 2459)), evalpoly(α², (-181, 73))
+    x7 = d⁴ * d² * evalpoly(jak², (4*c1, 2*c2, 36*c3, 657)) / 2835
+    w7 = 4 * d⁴ * d² * evalpoly(jak², (c1, c2, 27*c3, 657)) / 2835
+    c1, c2, = evalpoly(α², (6029959, -10087180, 5095482, -1146220, 107959)), evalpoly(α², (-2201939, 1678761, -507801, 63299))
+    c3, c4 = evalpoly(α², (1456807, -729422, 125671)), evalpoly(α², (-2879, 887))
+    x9 = d⁸ * evalpoly(jak², (2*c1, 3*c2, c3, 60*c4, 10644)) / 42525
+    w9 = d⁸ * evalpoly(jak², (2*c1, 6*c2, 3*c3, 240*c4, 53220)) / 42525
+    return (x3, x5, x7, x9), (w3, w5, w7, w9)
+end
 
 # And for α = 0:
-gl_bessel_x3(jak, d) = (jak^2 - 2)*d^2 / 3
-gl_bessel_x5(jak, d) = (11*jak^4 - 57*jak^2 + 94)d^4 / 45
-gl_bessel_x7(jak, d) = (657*jak^6 - 6516*jak^4 + 28102*jak^2 - 48308)*d^6 / 2835
-gl_bessel_x9(jak, d) = (10644*jak^8 - 172740*jak^6 + 1456807*jak^4 -  6605817*jak^2 + 12059918)*d^8 / 42525
-gl_bessel_x11(jak, d) = (410649*jak^10 -  9908262*jak^8 + 138902061*jak^6 - 1248722004*jak^4 + 6028914206*jak^2 - 11427291076)*d^10 / 1403325
-
-gl_bessel_w3(jak, d) = (jak^2 - 1)*2*d^2 / 3
-gl_bessel_w5(jak, d) = (33*jak^4 -114*jak^2 + 94)*d^4 / 45
-gl_bessel_w7(jak, d) = (657*jak^6 - 4887*jak^4 + 14051*jak^2 - 12077)*4*d^6 / 2835
-gl_bessel_w9(jak, d) = (53220*jak^8 - 690960*jak^6 + 4370421*jak^4 - 13211634*jak^2 + 12059918)*d^8 / 42525
-gl_bessel_w11(jak, d) = (1231947*jak^10 - 24770655*jak^8 + 277804122*jak^6 - 1873083006*jak^4 + 6028914206*jak^2 - 5713645538)*2*d^10 / 1403325
-
+function gl_bessel(jak, d)
+    jak² = jak * jak
+    d² = d * d
+    d⁴ = d² * d²
+    d⁸ = d⁴ * d⁴
+
+    x3 = d² * evalpoly(jak², (-2, 1)) / 3
+    x5 = d⁴ * evalpoly(jak², (94, -57, 11)) / 45
+    x7 = d⁴ * d² * evalpoly(jak², (-48308, 28102, -6516, 657)) / 2835
+    x9 = d⁸ * evalpoly(jak², (12059918, -6605817, 1456807, -172740, 10644)) / 42525
+    x11 = d⁸ * d² * evalpoly(jak², (-11427291076, 6028914206, -1248722004, 138902061, -9908262, 410649)) / 1403325
+
+    w3 = 2 * d² * evalpoly(jak², (-1, 1)) / 3
+    w5 = d⁴ * evalpoly(jak², (94, -114, 33)) / 45
+    w7 = 4 * d⁴ * d² * evalpoly(jak², (-12077, 14051, -4887, 657)) / 2835
+    w9 = d⁸ * evalpoly(jak², (12059918, -13211634, 4370421, -690960, 53220)) / 42525
+    w11 = 2 * d⁸ * d² * evalpoly(jak², (-5713645538, 6028914206, -1873083006, 277804122, -24770655, 1231947)) / 1403325
+    return (x3, x5, x7, x9, x11), (w3, w5, w7, w9, w11)
+end
 
 ## The soft edge (Airy region)
-
-gl_airy_x1(ak, d, α) = 1/d + ak*(d/4)^(-1/3)
-gl_airy_x3(ak, d, α) = ak^2*(d*16)^(1/3)/5 + (11/35-α^2-12/175*ak^3)*d + (16/1575*ak+92/7875*ak^4)*2^(2/3)*d^(5/3)
-gl_airy_x5(ak, d, α) = -(15152/3031875*ak^5+1088/121275*ak^2)*2^(1/3)*d^(7/3)
-
-gl_airy_x1(ak, d) = 1/d + ak*(d/4)^(-1/3)
-gl_airy_x3(ak, d) = ak^2*(d*16)^(1/3)/5 + (11/35-12/175*ak^3)*d + (16/1575*ak+92/7875*ak^4)*2^(2/3)*d^(5/3)
-gl_airy_x5(ak, d) = -(15152/3031875*ak^5+1088/121275*ak^2)*2^(1/3)*d^(7/3)
-
+function gl_airy(ak, d, α)
+    T = eltype(α)
+    cbrt_d = cbrt(d)
+    cbrt_d2 = cbrt_d * cbrt_d
+    ak² = ak * ak
+    ak³ = ak² * ak
+    ak⁴ = ak² * ak²
+
+    x1 = 1 / d + ak  * cbrt(T(4)) / cbrt_d
+    x3 = ak² * cbrt(T(16)) / 5 * cbrt_d + (T(11)/35 - α^2 - T(12)/175*ak³)*d + (T(16)/1575 * ak + T(92)/7875*ak⁴) * cbrt(T(4)) * d * cbrt_d2
+    x5 = -(T(15152)/3031875 * ak⁴ * ak + T(1088)/121275 * ak²) * cbrt(T(2)) * d * d * cbrt_d
+    return x1, x3, x5
+end
 
 # Sum the first Q elements of vals, and return the absolute value of the next
 # value in the list (or of the last value in the list)
@@ -420,16 +429,7 @@ function gausslaguerre_asy_bulk(n, α, k, d, T)
     end
 
     t = gl_bulk_solve_t(n, k, d)
-    x3 = gl_bulk_x3(t, d, α)
-    x5 = gl_bulk_x5(t, d, α)
-    x7 = gl_bulk_x7(t, d, α)
-    x9 = gl_bulk_x9(t, d, α)
-    w3 = gl_bulk_w3(t, d, α)
-    w5 = gl_bulk_w5(t, d, α)
-    w7 = gl_bulk_w7(t, d, α)
-
-    xs = (x3, x5, x7, x9)
-    ws = (w3, w5, w7)
+    xs, ws = gl_bulk(t, d, α)
 
     xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
     wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)
@@ -446,16 +446,8 @@ end
 
 function gausslaguerre_asy0_bulk(n, k, d, T)
     t = gl_bulk_solve_t(n, k, d)
-    x3 = gl_bulk_x3(t, d)
-    x5 = gl_bulk_x5(t, d)
-    x7 = gl_bulk_x7(t, d)
-    x9 = gl_bulk_x9(t, d)
-    w3 = gl_bulk_w3(t, d)
-    w5 = gl_bulk_w5(t, d)
-    w7 = gl_bulk_w7(t, d)
+    xs, ws = gl_bulk(t, d)
 
-    xs = (x3, x5, x7, x9)
-    ws = (w3, w5, w7)
 
     xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
     wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)
@@ -474,17 +466,8 @@ function gausslaguerre_asy_bessel(n, α, jak, d, T)
     if α == 0
         return gausslaguerre_asy0_bessel(n, jak, d, T)
     end
-    x3 = gl_bessel_x3(jak, d, α)
-    x5 = gl_bessel_x5(jak, d, α)
-    x7 = gl_bessel_x7(jak, d, α)
-    x9 = gl_bessel_x9(jak, d, α)
-    w3 = gl_bessel_w3(jak, d, α)
-    w5 = gl_bessel_w5(jak, d, α)
-    w7 = gl_bessel_w7(jak, d, α)
-    w9 = gl_bessel_w9(jak, d, α)
-
-    xs = (x3, x5, x7, x9)
-    ws = (w3, w5, w7, w9)
+
+    xs, ws = gl_bessel(jak, d, α)
 
     xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
     wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)
@@ -503,20 +486,8 @@ function gausslaguerre_asy_bessel(n, α, jak, d, T)
 end
 
 function gausslaguerre_asy0_bessel(n, jak, d, T)
-    x3 = gl_bessel_x3(jak, d)
-    x5 = gl_bessel_x5(jak, d)
-    x7 = gl_bessel_x7(jak, d)
-    x9 = gl_bessel_x9(jak, d)
-    x11 = gl_bessel_x11(jak, d)
-    w3 = gl_bessel_w3(jak, d)
-    w5 = gl_bessel_w5(jak, d)
-    w7 = gl_bessel_w7(jak, d)
-    w9 = gl_bessel_w9(jak, d)
-    w11 = gl_bessel_w11(jak, d)
-
-    xs = (x3, x5, x7, x9, x11)
-    ws = (w3, w5, w7, w9, w11)
-
+    xs, ws = gl_bessel(jak, d)
+
     xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
     wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)
 
@@ -537,22 +508,14 @@ function compute_airy_root(n, k)
         ak = AIRY_ROOTS[index]
     else
         t = 3 * π/2 * (index-0.25)
-        ak = -t^(2/3)*(1 + 5/48/t^2 - 5/36/t^4 + 77125/82944/t^6 -10856875/6967296/t^8)
+        ak = -t^(2/3) * evalpoly(inv(t*t), (6967296, 725760, -967680, 6478500, -10856875)) / 6967296
     end
     ak
 end
 
 function gausslaguerre_asy_airy(n, α, k, d, T)
-    if α == 0
-        return gausslaguerre_asy0_airy(n, k, d, T)
-    end
-
     ak = compute_airy_root(n, k)
-    x1 = gl_airy_x1(ak, d, α)
-    x3 = gl_airy_x3(ak, d, α)
-    x5 = gl_airy_x5(ak, d, α)
-
-    xs = (x1, x3, x5)
+    xs = gl_airy(ak, d, α)
 
     xk, xdelta = (T > 0) ? sum_explicit(xs, (T+1)>>1) : sum_adaptive(xs)
 
@@ -562,24 +525,6 @@ function gausslaguerre_asy_airy(n, α, k, d, T)
     return xk, wk, max(xdelta,wdelta)
 end
 
-function gausslaguerre_asy0_airy(n, k, d, T)
-    ak = compute_airy_root(n, k)
-
-    x1 = gl_airy_x1(ak, d)
-    x3 = gl_airy_x3(ak, d)
-    x5 = gl_airy_x5(ak, d)
-
-    xs = (x1, x3, x5)
-
-    xk, xdelta = (T > 0) ? sum_explicit(xs, (T+1)>>1) : sum_adaptive(xs)
-
-    wk = 4^(1/3) * xk^(1/3) * exp(-xk) / (airyaiprime(ak))^2
-    wdelta = abs(wk)
-
-    return xk, wk, max(xdelta,wdelta)
-end
-
-
 """
 Calculate Gauss-Laguerre nodes and weights from the eigenvalue decomposition of
 the Jacobi matrix.