optimize U-polynomials

Author: heltonmc

src/U_polynomials.jl

@@ -1,81 +1,102 @@
-function Uk_poly_Kn(p, v, p2, ::Type{Float32})
-    u0 = one(p)
-    u1 = 1 / 24 * evalpoly(p2, (3, -5))
-    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
-    return evalpoly(-p/v, (u0, u1, u2))
-end
-function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: Float64
-    u0 = one(T)
-    u1 = 1 / 24 * evalpoly(p2, (3, -5))
-    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
-    u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425))
-    u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4))
-end
-function Uk_poly_In(p, v, p2, ::Type{T}) where T <: Float64
-    u0 = one(T)
-    u1 = -1 / 24 * evalpoly(p2, (3, -5))
-    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
-    u3 = -1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425))
-    u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4))
-end
-function Uk_poly_In(p, v, p2, ::Type{Float32})
-    u0 = one(p)
-    u1 = -1 / 24 * evalpoly(p2, (3, -5))
-    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
-    return evalpoly(-p/v, (u0, u1, u2))
- end
- 
-
-function Uk_poly_Jn(p, v, p2, x, ::Type{T}) where T <: Float64
+function Uk_poly_Jn(p, v, p2, x::T) where T <: Float64
     if v > 5.0 + 1.00033*x + (1427.61*x)^(1/3)
-        return Uk_poly_Jn10(p, v, p2, T)
+        return Uk_poly11(p, v, p2)[1]
     else
-        return Uk_poly_Jn20(p, v, p2, T)
+        return Uk_poly21(p, v, p2)[1]
     end
 end
 
-function Uk_poly_Jn10(p, v, p2, ::Type{T}) where T <: Float64
-    u0 = one(T)
-    u1 = -evalpoly(p2, (0.125, -0.20833333333333334))
+Uk_poly_In(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[1]
+Uk_poly_In(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[1]
+Uk_poly_Kn(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[2]
+Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[2]
+
+@inline function split_evalpoly(x, P)
+    N = length(P)
+    xx = x*x
+
+    out = P[end]
+    out2 = P[end-1]
+
+    for i in N-2:-2:1
+        out = muladd(xx, out, P[i])
+        out2 = muladd(xx, out2, P[i-1])
+    end
+
+    out = x*out
+    return out2 - out, out2 + out
+end
+function Uk_poly3(p, v, p2)
+    u0 = 1.0
+    u1 = evalpoly(p2, (0.125, -0.20833333333333334))
     u2 = evalpoly(p2, (0.0703125, -0.4010416666666667, 0.3342013888888889))
-    u3 = -evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
+    u3 = evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
+
+    Poly = (u0, u1, u2, u3)
+
+   return split_evalpoly(-p/v, Poly)
+end
+function Uk_poly5(p, v, p2)
+    u0 = 1.0
+    u1 = evalpoly(p2, (0.125, -0.20833333333333334))
+    u2 = evalpoly(p2, (0.0703125, -0.4010416666666667, 0.3342013888888889))
+    u3 = evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
     u4 = evalpoly(p2, (0.112152099609375, -2.3640869140625, 8.78912353515625, -11.207002616222994, 4.669584423426247))
-    u5 = -evalpoly(p2, (0.22710800170898438, -7.368794359479632, 42.53499874538846, -91.81824154324002, 84.63621767460073, -28.212072558200244))
+    u5 = evalpoly(p2, (0.22710800170898438, -7.368794359479632, 42.53499874538846, -91.81824154324002, 84.63621767460073, -28.212072558200244))
+    
+    Poly = (u0, u1, u2, u3, u4, u5)
+    return split_evalpoly(-p/v, Poly)
+end
+
+function Uk_poly11(p, v, p2)
+    u0 = 1.0
+    u1 = evalpoly(p2, (0.125, -0.20833333333333334))
+    u2 = evalpoly(p2, (0.0703125, -0.4010416666666667, 0.3342013888888889))
+    u3 = evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
+    u4 = evalpoly(p2, (0.112152099609375, -2.3640869140625, 8.78912353515625, -11.207002616222994, 4.669584423426247))
+    u5 = evalpoly(p2, (0.22710800170898438, -7.368794359479632, 42.53499874538846, -91.81824154324002, 84.63621767460073, -28.212072558200244))
     u6 = evalpoly(p2, (0.5725014209747314, -26.491430486951554, 218.1905117442116, -699.5796273761325, 1059.9904525279999, -765.2524681411817, 212.57013003921713))
-    u7 = -evalpoly(p2, (1.7277275025844574, -108.09091978839466, 1200.9029132163525, -5305.646978613403, 11655.393336864534, -13586.550006434136, 8061.722181737309, -1919.457662318407))
+    u7 = evalpoly(p2, (1.7277275025844574, -108.09091978839466, 1200.9029132163525, -5305.646978613403, 11655.393336864534, -13586.550006434136, 8061.722181737309, -1919.457662318407))
     u8 = evalpoly(p2, (6.074042001273483, -493.91530477308805, 7109.514302489364, -41192.65496889755, 122200.46498301747, -203400.17728041555, 192547.00123253153, -96980.59838863752, 20204.29133096615))
-    u9 = -evalpoly(p2, (24.380529699556064, -2499.8304818112097, 45218.76898136273, -331645.17248456355, 1.2683652733216248e6, -2.8135632265865337e6, 3.763271297656404e6, -2.998015918538107e6, 1.3117636146629772e6, -242919.18790055133))
+    u9 = evalpoly(p2, (24.380529699556064, -2499.8304818112097, 45218.76898136273, -331645.17248456355, 1.2683652733216248e6, -2.8135632265865337e6, 3.763271297656404e6, -2.998015918538107e6, 1.3117636146629772e6, -242919.18790055133))
     u10 = evalpoly(p2, (110.01714026924674, -13886.08975371704, 308186.40461266245, -2.7856181280864547e6, 1.3288767166421818e7, -3.7567176660763346e7, 6.634451227472903e7, -7.410514821153265e7, 5.095260249266464e7, -1.9706819118432228e7, 3.284469853072038e6))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10))
+    u11 = evalpoly(p2, (551.3358961220206, -84005.4336030241, 2.2437681779224495e6, -2.4474062725738727e7, 1.420629077975331e8, -4.9588978427503026e8, 1.1068428168230145e9, -1.6210805521083372e9, 1.5535968995705802e9, -9.394623596815784e8, 3.2557307418576574e8, -4.932925366450996e7))
+
+    Poly = (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11)
+
+    return split_evalpoly(-p/v, Poly)
 end
-function Uk_poly_Jn20(p, v, p2, ::Type{T}) where T <: Float64
-    u0 = one(T)
-    u1 = -evalpoly(p2, (0.125, -0.20833333333333334))
+
+function Uk_poly21(p, v, p2)
+    u0 = 1.0
+    u1 = evalpoly(p2, (0.125, -0.20833333333333334))
     u2 = evalpoly(p2, (0.0703125, -0.4010416666666667, 0.3342013888888889))
-    u3 = -evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
+    u3 = evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173))
     u4 = evalpoly(p2, (0.112152099609375, -2.3640869140625, 8.78912353515625, -11.207002616222994, 4.669584423426247))
-    u5 = -evalpoly(p2, (0.22710800170898438, -7.368794359479632, 42.53499874538846, -91.81824154324002, 84.63621767460073, -28.212072558200244))
+    u5 = evalpoly(p2, (0.22710800170898438, -7.368794359479632, 42.53499874538846, -91.81824154324002, 84.63621767460073, -28.212072558200244))
     u6 = evalpoly(p2, (0.5725014209747314, -26.491430486951554, 218.1905117442116, -699.5796273761325, 1059.9904525279999, -765.2524681411817, 212.57013003921713))
-    u7 = -evalpoly(p2, (1.7277275025844574, -108.09091978839466, 1200.9029132163525, -5305.646978613403, 11655.393336864534, -13586.550006434136, 8061.722181737309, -1919.457662318407))
+    u7 = evalpoly(p2, (1.7277275025844574, -108.09091978839466, 1200.9029132163525, -5305.646978613403, 11655.393336864534, -13586.550006434136, 8061.722181737309, -1919.457662318407))
     u8 = evalpoly(p2, (6.074042001273483, -493.91530477308805, 7109.514302489364, -41192.65496889755, 122200.46498301747, -203400.17728041555, 192547.00123253153, -96980.59838863752, 20204.29133096615))
-    u9 = -evalpoly(p2, (24.380529699556064, -2499.8304818112097, 45218.76898136273, -331645.17248456355, 1.2683652733216248e6, -2.8135632265865337e6, 3.763271297656404e6, -2.998015918538107e6, 1.3117636146629772e6, -242919.18790055133))
+    u9 = evalpoly(p2, (24.380529699556064, -2499.8304818112097, 45218.76898136273, -331645.17248456355, 1.2683652733216248e6, -2.8135632265865337e6, 3.763271297656404e6, -2.998015918538107e6, 1.3117636146629772e6, -242919.18790055133))
     u10 = evalpoly(p2, (110.01714026924674, -13886.08975371704, 308186.40461266245, -2.7856181280864547e6, 1.3288767166421818e7, -3.7567176660763346e7, 6.634451227472903e7, -7.410514821153265e7, 5.095260249266464e7, -1.9706819118432228e7, 3.284469853072038e6))
-    u11 = -evalpoly(p2, (551.3358961220206, -84005.4336030241, 2.2437681779224495e6, -2.4474062725738727e7, 1.420629077975331e8, -4.9588978427503026e8, 1.1068428168230145e9, -1.6210805521083372e9, 1.5535968995705802e9, -9.394623596815784e8, 3.2557307418576574e8, -4.932925366450996e7))
+    u11 = evalpoly(p2, (551.3358961220206, -84005.4336030241, 2.2437681779224495e6, -2.4474062725738727e7, 1.420629077975331e8, -4.9588978427503026e8, 1.1068428168230145e9, -1.6210805521083372e9, 1.5535968995705802e9, -9.394623596815784e8, 3.2557307418576574e8, -4.932925366450996e7))
     u12 = evalpoly(p2, (3038.090510922384, -549842.3275722887, 1.7395107553978164e7, -2.2510566188941526e8, 1.5592798648792574e9, -6.563293792619285e9, 1.79542137311556e10, -3.3026599749800724e10, 4.1280185579753975e10, -3.4632043388158775e10, 1.8688207509295826e10, -5.866481492051847e9, 8.147890961183121e8))
-    u13 = -evalpoly(p2, (18257.755474293175, -3.8718334425726123e6, 1.43157876718889e8, -2.167164983223795e9, 1.763473060683497e10, -8.786707217802327e10, 2.879006499061506e11, -6.453648692453765e11, 1.0081581068653821e12, -1.0983751560812233e12, 8.192186695485773e11, -3.990961752244665e11, 1.144982377320258e11, -1.4679261247695616e10))
+    u13 = evalpoly(p2, (18257.755474293175, -3.8718334425726123e6, 1.43157876718889e8, -2.167164983223795e9, 1.763473060683497e10, -8.786707217802327e10, 2.879006499061506e11, -6.453648692453765e11, 1.0081581068653821e12, -1.0983751560812233e12, 8.192186695485773e11, -3.990961752244665e11, 1.144982377320258e11, -1.4679261247695616e10))
     u14 = evalpoly(p2, (118838.42625678326, -2.9188388122220814e7, 1.2470092935127103e9, -2.1822927757529224e10, 2.0591450323241e11, -1.1965528801961816e12, 4.612725780849132e12, -1.2320491305598287e13, 2.334836404458184e13, -3.166708858478516e13, 3.056512551993532e13, -2.0516899410934438e13, 9.109341185239898e12, -2.406297900028504e12, 2.86464035717679e11))   
-    u15 = -evalpoly(p2, (832859.3040162893, -2.3455796352225152e8, 1.1465754899448236e10, -2.2961937296824646e11, 2.4850009280340854e12, -1.663482472489248e13, 7.437312290867914e13, -2.3260483118893994e14, 5.230548825784446e14, -8.57461032982895e14, 1.0269551960827625e15, -8.894969398810265e14, 5.4273966498765975e14, -2.213496387025252e14, 5.417751075510605e13, -6.019723417234006e12))
+    u15 = evalpoly(p2, (832859.3040162893, -2.3455796352225152e8, 1.1465754899448236e10, -2.2961937296824646e11, 2.4850009280340854e12, -1.663482472489248e13, 7.437312290867914e13, -2.3260483118893994e14, 5.230548825784446e14, -8.57461032982895e14, 1.0269551960827625e15, -8.894969398810265e14, 5.4273966498765975e14, -2.213496387025252e14, 5.417751075510605e13, -6.019723417234006e12))
     u16 = evalpoly(p2, (6.252951493434797e6, -2.0016469281917763e9, 1.1099740513917902e11, -2.5215584749128545e12, 3.100743647289646e13, -2.3665253045164925e14, 1.2126758042503475e15, -4.3793258383640155e15, 1.1486706978449752e16, -2.2268225133911144e16, 3.213827526858624e16, -3.4447226006485144e16, 2.705471130619708e16, -1.5129826322457682e16, 5.705782159023671e15, -1.3010127235496995e15, 1.3552215870309369e14))
-    u17 = -evalpoly(p2, (5.0069589531988926e7, -1.8078220384658062e10, 1.128709145410874e12, -2.886383763141476e13, 4.0004445704303625e14, -3.4503855118462725e15, 2.0064271476309532e16, -8.270945651585064e16, 2.4960365126160426e17, -5.62631788074636e17, 9.575335098169139e17, -1.2336116931960694e18, 1.1961991142756308e18, -8.592577980317548e17, 4.4347954614171904e17, -1.5552983504313904e17, 3.3192764720355224e16, -3.254192619642669e15))
+    u17 = evalpoly(p2, (5.0069589531988926e7, -1.8078220384658062e10, 1.128709145410874e12, -2.886383763141476e13, 4.0004445704303625e14, -3.4503855118462725e15, 2.0064271476309532e16, -8.270945651585064e16, 2.4960365126160426e17, -5.62631788074636e17, 9.575335098169139e17, -1.2336116931960694e18, 1.1961991142756308e18, -8.592577980317548e17, 4.4347954614171904e17, -1.5552983504313904e17, 3.3192764720355224e16, -3.254192619642669e15))
     u18 = evalpoly(p2, (4.259392165047669e8, -1.722832387173505e11, 1.2030115826419191e13, -3.4396530474307594e14, 5.335106978708839e15, -5.1605093193485224e16, 3.37667624979061e17, -1.5736434765189599e18, 5.402894876715982e18, -1.3970803516443374e19, 2.757282981650519e19, -4.178861444656839e19, 4.859942729324836e19, -4.301555703831444e19, 2.846521225167657e19, -1.3639420410571592e19, 4.47020096401231e18, -8.966114215270463e17, 8.30195760673191e16))
-    u19 = -evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
+    u19 = evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
     u20 = evalpoly(p2, (3.646840080706556e10, -1.818726203851104e13, 1.5613123930484672e15, -5.48403360388329e16, 1.0461721131134344e18, -1.2483700995047234e19, 1.0126774169536592e20, -5.8917941350694964e20, 2.548961114664972e21, -8.405915817108351e21, 2.1487414815055883e22, -4.302534303482379e22, 6.783661642951883e22, -8.423222750084323e22, 8.19433100543513e22, -6.173206302884415e22, 3.528435843903409e22, -1.4787743528433614e22, 4.285296082829494e21, -7.671943936729004e20, 6.393286613940837e19))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20))
+    u21 = evalpoly(p2, (3.6490108188498334e11, -2.0052440123627112e14, 1.894406984252143e16, -7.319501491566134e17, 1.5365025218443373e19, -2.0197335419300872e20, 1.8081594057131945e21, -1.1640246461465369e22, 5.591591380366263e22, -2.0566149136271542e23, 5.8965434619782445e23, -1.3337178907798302e24, 2.3967237744351682e24, -3.430872898515746e24, 3.905264103536985e24, -3.511096528332644e24, 2.461506085403875e24, -1.3170969618092387e24, 5.194289094766812e23, -1.4228394823321413e23, 2.417461500896379e22, -1.91862023880665e21))
+    
+    Poly = (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20, u21)
+
+    return split_evalpoly(-p/v, Poly)
 end
 
-function Uk_poly_Jn(p, v, p2, x, ::Type{T}) where T <: BigFloat
+function Uk_poly_Jn(p, v, p2, x::T) where T <: BigFloat
     u0 = one(T)
     u1 = -evalpoly(p2, (3, -5)) / 24
     u2 = evalpoly(p2, (81, -462, 385)) / 1152
@@ -101,22 +122,6 @@ function Uk_poly_Jn(p, v, p2, x, ::Type{T}) where T <: BigFloat
     return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16, u17, u18, u19, u20, u21))
 end
 #=
-function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat
-    u0 = one(T)
-    u1 = 1 / 24 * evalpoly(p2, (3, -5))
-    u2 = 1 / 1152 * evalpoly(p2, (81, -462, 385))
-    u3 = 1 / 414720 * evalpoly(p2, (30375, -369603, 765765, -425425))
-    u4 = 1 / 39813120 * evalpoly(p2, (4465125, -94121676, 349922430, -446185740, 185910725))
-    u5 = 1 / 6688604160 * evalpoly(p2, (1519035525, -49286948607, 284499769554, -614135872350, 566098157625, -188699385875))
-    u6 = 1 / 4815794995200 * evalpoly(p2, (2757049477875, -127577298354750, 1050760774457901, -3369032068261860,5104696716244125, -3685299006138750, 1023694168371875))
-    u7 = 1 / 115579079884800 * evalpoly(p2, (199689155040375, -12493049053044375, 138799253740521843, -613221795981706275, 1347119637570231525, -1570320948552481125, 931766432052080625,  -221849150488590625))
-    u8 = 1 / 22191183337881600 * evalpoly(p2, (134790179652253125, -10960565081605263000, 157768535329832893644, -914113758588905038248, 2711772922412520971550, -4513690624987320777000, 4272845805510421639500, -2152114239059719935000, 448357133137441653125))
-    u9 = 1 / 263631258054033408000 * evalpoly(p2, (6427469716717690265625, -659033454841709672064375, 11921080954211358275362500, -87432034049652400520788332, 334380732677827878090447630, -741743213039573443221773250, 992115946599792610768672500, -790370708270219620781737500, 345821892003106984030190625, -64041091111686478524109375))
-    u10 = 1 / 88580102706155225088000 * evalpoly(p2, (9745329584487361980740625, -1230031256571145165088463750, 27299183373230345667273718125, -246750339886026017414509498824, 1177120360439828012193658602930, -3327704366990695147540934069220, 5876803711285273203043452095250, -6564241639632418015173104205000, 4513386761946134740461797128125, -1745632061522350031610173343750, 290938676920391671935028890625))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10))
-end
-
-
     u0 = one(x)
     u1 = p / 24 * (3 - 5*p^2) * -1 / v
     u2 = p^2 / 1152 * (81 - 462*p^2 + 385*p^4) / v^2
@@ -144,8 +149,6 @@ end
         + 5136561256208409671660362298619778869859994724706875*p^16 - 2284251621937242886581917667066122422330060024456125*p^14 + 730367145705123976114617970888707594104468177381925*p^12 - 163359140754958502896104062604202934925448173291477*p^10
         + 24403480234538299231733883413666768614198435948125*p^8 - 2254933495791765108580529087615802250458013685625*p^6 + 112597271053778779753048514469995937998172890625*p^4 - 2303431987527333128955769182299845911305390625*p^2 + 8178936810213560828419581728001773291015625) * -1 / v^15
 
-
-
 function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
     u0 = one(T)
     u1 = -evalpoly(p2, (0.125, -0.20833333333333334))

src/besselj.jl

@@ -267,7 +267,7 @@ function besselj_debye(v, x)
     p = v / vs
     p2  = v^2 / vmx
 
-    return coef * Uk_poly_Jn(p, v, p2, x, T)
+    return coef * Uk_poly_Jn(p, v, p2, x)
 end
 
 # For 0.0 <= x < 171.5