re-add constants

Author: oscardssmith

src/bessely.jl

@@ -15,17 +15,17 @@
 #    Hankel's asymptotic expansion is used
 #    where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
 #    See section 4 of [3] for more details and [1] for coefficients of polynomials
-# 
+#
 #   Branch 3: x >= 25.0
 #              bessely0 = sqrt(2/(pi*x))*beta(x)*(sin(x - pi/4 - alpha(x))
 #   See modified expansions given in [3]. Exact coefficients are used.
 #
 #   Calculation of bessely1 is done in a similar way as bessely0.
 #   See [3] for details on similarities.
-# 
+#
 # [1] https://github.com/deepmind/torch-cephes
 # [2] Cephes Math Library Release 2.8:  June, 2000 by Stephen L. Moshier
-# [3] Harrison, John. "Fast and accurate Bessel function computation." 
+# [3] Harrison, John. "Fast and accurate Bessel function computation."
 #     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
 #
 
@@ -63,19 +63,22 @@ function _bessely0_compute(x::Float64)
         p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
+        xinv = inv(x)
+        x2 = xinv*xinv
         if x < 120.0
-            p = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
-            q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
+            p1 = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
+            q1 = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
+            p = evalpoly(x2, p1)
+            q = evalpoly(x2, q1)
         else
-            p = (one(T), -1/16, 53/512, -4447/8192)
-            q = (-1/8, 25/384, -1073/5120, 375733/229376)
+            p2 = (one(T), -1/16, 53/512, -4447/8192)
+            q2 = (-1/8, 25/384, -1073/5120, 375733/229376)
+            p = evalpoly(x2, p2)
+            evalpoly(x2, q2)
         end
-        xinv = inv(x)
-        x2 = xinv*xinv
 
-        p = evalpoly(x2, p)
         a = SQ2OPI(T) * sqrt(xinv) * p
-        xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))
+        xn = muladd(xinv, q, - PIO4(T))
 
         # the following computes b = sin(x + xn) more accurately
         # see src/misc.jl
@@ -136,19 +139,22 @@ function _bessely1_compute(x::Float64)
         p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
+        xinv = inv(x)
+        x2 = xinv*xinv
         if x < 120.0
-            p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296) 
-            q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
+            p1 = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
+            q1 = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
+            p = evalpoly(x2, p1)
+            q = evalpoly(x2, q1)
         else
-            p = (one(T), 3/16, -99/512, 6597/8192)
-            q = (3/8, -21/128, 1899/5120, -543483/229376)
+            p2 = (one(T), 3/16, -99/512, 6597/8192)
+            q2 = (3/8, -21/128, 1899/5120, -543483/229376)
+            p = evalpoly(x2, q2)
+            q = evalpoly(x2, q2)
         end
-        xinv = inv(x)
-        x2 = xinv*xinv
 
-        p = evalpoly(x2, p)
         a = SQ2OPI(T) * sqrt(xinv) * p
-        xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
+        xn = muladd(xinv, q, - 3 * PIO4(T))
 
         # the following computes b = sin(x + xn) more accurately
         # see src/misc.jl

src/constants.jl

@@ -28,6 +28,30 @@ const QQ_y0(::Type{Float64}) = (
     7.24046774195652478189E3, 3.88240183605401609683E3, 8.56430025976980587198E2,
     6.43178256118178023184E1, 1.00000000000000000000E0
 )
+const PP_j1(::Type{Float64}) = (
+    1.00000000000000000254E0, 5.21451598682361504063E0,
+    8.42404590141772420927E0, 5.11207951146807644818E0,
+    1.12719608129684925192E0, 7.31397056940917570436E-2,
+    7.62125616208173112003E-4
+)
+const PQ_j1(::Type{Float64}) = (
+    9.99999999999999997461E-1, 5.20982848682361821619E0,
+    8.39985554327604159757E0, 5.07386386128601488557E0,
+    1.10514232634061696926E0, 6.88455908754495404082E-2,
+    5.71323128072548699714E-4
+)
+const QP_j1(::Type{Float64}) = (
+    2.52070205858023719784E1, 2.11688757100572135698E2,
+    5.97489612400613639965E2, 7.10856304998926107277E2,
+    3.66779609360150777800E2, 7.58238284132545283818E1,
+    4.98213872951233449420E0, 5.10862594750176621635E-2
+)
+const QQ_j1(::Type{Float64}) = (
+    3.36093607810698293419E2, 2.82619278517639096600E3,
+    7.99704160447350683650E3, 9.56231892404756170795E3,
+    4.98641058337653607651E3, 1.05644886038262816351E3,
+    7.42373277035675149943E1, 1.00000000000000000000E0
+)
 const YP_y1(::Type{Float64}) = (
     -7.78877196265950026825E17, 2.02439475713594898196E17,
     -8.12770255501325109621E15, 1.14509511541823727583E14,
@@ -50,7 +74,7 @@ const besseli0_small_coefs(::Type{Float64}) = (
     0.001736111111149161, 6.94444444107536e-5, 1.9290123635366806e-6,
     3.9367592765038015e-8, 6.151201574092085e-10, 7.593827956729909e-12,
     7.596677643342155e-14, 6.255282299620455e-16, 4.470993793303175e-18,
-    2.1859737023077178e-20, 2.0941557335286373e-22    
+    2.1859737023077178e-20, 2.0941557335286373e-22
 )
 const besseli0_med_coefs(::Type{Float32}) = (
     0.3989423094522082f0, 0.049857687814784155f0, 0.028606362029333827f0,