fix format

Author: heltonmc

src/besseli.jl

@@ -135,30 +135,30 @@ end
 #####
 
 function besseli0(z::ComplexF64)
-    check = false
+    isconj = false
     if real(z) < 0.0
         z = -z
     end
     if imag(z) < 0.0
         z = conj(z)
-        check = true
+        isconj = true
     end
 
     if abs(z) > 17.5
         zinv = 1 / z
         e = exp(z)
-        s = sqrt(2 * z * π)
+        sinv =sqrt(zinv)/sqrt(2*pi)
         p = evalpoly(zinv*zinv, (1.0, 0.0703125, 0.112152099609375, 0.5725014209747314, 6.074042001273483, 110.01714026924674, 3038.090510922384, 118838.42625678325, 6.252951493434797e6, 4.259392165047669e8, 3.646840080706556e10, 3.8335346613939443e12, 4.8540146868529006e14, 7.286857349377656e16))
         p2 = zinv*evalpoly(zinv*zinv, (0.125, 0.0732421875, 0.22710800170898438, 1.7277275025844574, 24.380529699556064, 551.3358961220206, 18257.755474293175, 832859.3040162893, 5.0069589531988926e7, 3.8362551802304335e9, 3.6490108188498334e11, 4.218971570284097e13, 5.827244631566907e15, 9.47628809926011e17))
-        r = e / s * (p + p2) + im * (p - p2) / (e * s)
+        r = e * sinv * (p + p2) + im * (p - p2) * sinv / e
     elseif abs(z) < 5.2
         r = evalpoly(z*z, (1.0, 0.25, 0.015625, 0.00043402777777777775, 6.781684027777777e-6, 6.781684027777778e-8, 4.709502797067901e-10, 2.4028075495244395e-12, 9.385966990329842e-15, 2.896903392077112e-17, 7.242258480192779e-20, 1.4963343967340453e-22, 2.5978027721077174e-25, 3.842903509035085e-28, 4.9016626390753635e-31, 5.4462918211948485e-34, 5.318644356635594e-37))
     else
         if angle(z) <= π / 4.4
             zz = z*z
-            zz2 = zz*zz
-            r = evalpoly(zz2, (1.0, 0.015625, 6.781684027777777e-6, 4.709502797067901e-10, 9.385966990329842e-15, 7.242258480192779e-20, 2.5978027721077174e-25, 4.9016626390753635e-31, 5.318644356635594e-37, 3.5500798014623073e-43, 1.5365650110207356e-49, 4.499321282809354e-56, 9.228877211181495e-63, 1.3652185223641266e-69, 1.4929270885431173e-76, 1.232764473958843e-83, 7.829549665715613e-91, 3.887148093924665e-98))
-            r += zz*evalpoly(zz2, (0.25, 0.00043402777777777775, 6.781684027777778e-8, 2.4028075495244395e-12, 2.896903392077112e-17, 1.4963343967340453e-22, 3.842903509035085e-28, 5.4462918211948485e-34, 4.60090342269515e-40, 2.458504017633177e-46, 8.71068600351891e-53, 2.1263333094562166e-59, 3.691550884472598e-66, 4.681819349671216e-73, 4.4379521062518346e-80, 3.206983543077115e-87, 1.7974172786307652e-94, 7.932955293723807e-102))
+            z4 = zz*zz
+            r = evalpoly(z4, (1.0, 0.015625, 6.781684027777777e-6, 4.709502797067901e-10, 9.385966990329842e-15, 7.242258480192779e-20, 2.5978027721077174e-25, 4.9016626390753635e-31, 5.318644356635594e-37, 3.5500798014623073e-43, 1.5365650110207356e-49, 4.499321282809354e-56, 9.228877211181495e-63, 1.3652185223641266e-69, 1.4929270885431173e-76, 1.232764473958843e-83, 7.829549665715613e-91, 3.887148093924665e-98))
+            r += zz*evalpoly(z4, (0.25, 0.00043402777777777775, 6.781684027777778e-8, 2.4028075495244395e-12, 2.896903392077112e-17, 1.4963343967340453e-22, 3.842903509035085e-28, 5.4462918211948485e-34, 4.60090342269515e-40, 2.458504017633177e-46, 8.71068600351891e-53, 2.1263333094562166e-59, 3.691550884472598e-66, 4.681819349671216e-73, 4.4379521062518346e-80, 3.206983543077115e-87, 1.7974172786307652e-94, 7.932955293723807e-102))
         else
             zz = (-1.3144540816651462 + 20.054693865816724im, -1.2152720220184496 + 18.54146805523844im, 
                   -1.0758582968968085 + 16.41442564500392im, -0.9186088936209557 + 14.015263371275323im, 
@@ -211,35 +211,35 @@ function besseli0(z::ComplexF64)
         r = s1 / (s2 * sqrt(z) * exp(-z))
         end
     end
-    check && (r = conj(r))
+    isconj && (r = conj(r))
     return r
 end
 
 function besseli0(z::ComplexF32)
     z = ComplexF64(z)
-    check = false
+    isconj = false
     if real(z) < 0.0
         z = -z
     end
     if imag(z) < 0.0
         z = conj(z)
-        check = true
+        isconj = true
     end
 
     if abs(z) > 10.0
         zinv = 1 / z
         e = exp(z)
-        s = sqrt(2 * z * π)
+        sinv = sqrt(zinv) * SQ1O2PI(Float64)
         p = evalpoly(zinv*zinv, (1.0, 0.0703125, 0.112152099609375, 0.5725014209747314, 6.074042001273483, 110.01714026924674))
         p2 = zinv*evalpoly(zinv*zinv, (0.125, 0.0732421875, 0.22710800170898438, 1.7277275025844574, 24.380529699556064, 551.3358961220206))
-        r = e / s * (p + p2) + im * (p - p2) / (e * s)
+        r = e * sinv * (p + p2) + im * sinv * (p - p2) / e
     else
         zz = z*z
-        zz2 = zz*zz
-        r = evalpoly(zz2, (1.0, 0.015625, 6.781684027777777e-6, 4.709502797067901e-10, 9.385966990329842e-15, 7.242258480192779e-20, 2.5978027721077174e-25, 4.9016626390753635e-31, 5.318644356635594e-37, 3.5500798014623073e-43))
-        r += zz*evalpoly(zz2, (0.25, 0.00043402777777777775, 6.781684027777778e-8, 2.4028075495244395e-12, 2.896903392077112e-17, 1.4963343967340453e-22, 3.842903509035085e-28, 5.4462918211948485e-34, 4.60090342269515e-40, 2.458504017633177e-46))
+        z4 = zz*zz
+        r = evalpoly(z4, (1.0, 0.015625, 6.781684027777777e-6, 4.709502797067901e-10, 9.385966990329842e-15, 7.242258480192779e-20, 2.5978027721077174e-25, 4.9016626390753635e-31, 5.318644356635594e-37, 3.5500798014623073e-43))
+        r += zz*evalpoly(z4, (0.25, 0.00043402777777777775, 6.781684027777778e-8, 2.4028075495244395e-12, 2.896903392077112e-17, 1.4963343967340453e-22, 3.842903509035085e-28, 5.4462918211948485e-34, 4.60090342269515e-40, 2.458504017633177e-46))
     end
-    check && (r = conj(r))
+    isconj && (r = conj(r))
     return ComplexF32(r)
 end
 
@@ -252,26 +252,26 @@ function besseli1(z::ComplexF64)
     end
     if imag(z) < 0.0
         z = conj(z)
-        check = true
+        isconj = true
     else
-        check = false
+        isconj = false
     end
 
     if abs(z) > 17.5
         zinv = 1 / z
         e = exp(z)
-        s = sqrt(2 * z * π)
+        sinv = sqrt(zinv) * SQ1O2PI(Float64)
         p = evalpoly(zinv*zinv, (1.0, -0.1171875, -0.144195556640625, -0.6765925884246826, -6.883914268109947, -121.59789187653587, -3302.2722944808525, -127641.2726461746, -6.656367718817688e6, -4.502786003050393e8, -3.8338575207427895e10, -4.0118385991331978e12, -5.060568503314726e14, -7.572616461117957e16, -1.3262572853205555e19))
         p2 = zinv*evalpoly(zinv*zinv, (0.375, 0.1025390625, 0.2775764465332031, 1.993531733751297, 27.248827311268542, 603.8440767050702, 19718.37591223663, 890297.8767070678, 5.310411010968523e7, 4.043620325107754e9, 3.827011346598606e11, 4.406481417852279e13, 6.065091351222699e15, 9.83388387659068e17))
-        r = e / s * (p - p2) - im * (p + p2) / (e * s)
+        r = e * sinv * (p - p2) - im * (p + p2) * sinv / e
     elseif abs(z) < 5.2
         r = z*evalpoly(z*z, (0.5, 0.0625, 0.0026041666666666665, 5.425347222222222e-5, 6.781684027777778e-7, 5.651403356481481e-9, 3.363930569334215e-11, 1.5017547184527747e-13, 5.214426105738801e-16, 1.4484516960385557e-18, 3.2919356728148996e-21, 6.234726653058522e-24, 9.991549123491221e-27, 1.3724655389411017e-29, 1.6338875463584545e-32, 1.7019661941233902e-35, 1.564307163716351e-38))
     else
         if angle(z) <= π / 4.4
             zz = z*z
-            zz2 = zz*zz
-            r = evalpoly(zz2, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45, 3.658488121477942e-51, 9.781133223498595e-58, 1.845775442236299e-64, 2.5281824488224564e-71, 2.574012221626064e-78, 1.9883297967078112e-85, 1.1862954038963049e-92, 5.553068705606664e-100))
-            r += zz*evalpoly(zz2, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48, 1.9797013644361157e-54, 4.4298610613671176e-61, 7.099136316293458e-68, 8.360391695841456e-75, 7.396586843753058e-82, 5.010911786057992e-89, 2.6432607038687723e-96))
+            z4 = zz*zz
+            r = evalpoly(z4, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45, 3.658488121477942e-51, 9.781133223498595e-58, 1.845775442236299e-64, 2.5281824488224564e-71, 2.574012221626064e-78, 1.9883297967078112e-85, 1.1862954038963049e-92, 5.553068705606664e-100))
+            r += zz*evalpoly(z4, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48, 1.9797013644361157e-54, 4.4298610613671176e-61, 7.099136316293458e-68, 8.360391695841456e-75, 7.396586843753058e-82, 5.010911786057992e-89, 2.6432607038687723e-96))
             r *= z    
         else
             zz = (-1.3144285719254842 + 20.054304662400146im, -1.2128746095534062 + 18.50489060934118im,
@@ -325,7 +325,7 @@ function besseli1(z::ComplexF64)
         r = s1 / (s2 * sqrt(z) * exp(-z))
         end
     end
-    check && (r = conj(r))
+    isconj && (r = conj(r))
     return r*c
 end
 
@@ -339,26 +339,26 @@ function besseli1(z::ComplexF32)
     end
     if imag(z) < 0.0
         z = conj(z)
-        check = true
+        isconj = true
     else
-        check = false
+        isconj = false
     end
 
     if abs(z) > 10.0
         zinv = 1 / z
         e = exp(z)
-        s = sqrt(2 * z * π)
+        sinv = sqrt(zinv) * SQ1O2PI(Float64)
         p = evalpoly(zinv*zinv, (1.0, -0.1171875, -0.144195556640625, -0.6765925884246826, -6.883914268109947, -121.59789187653587, -3302.2722944808525))
         p2 = zinv*evalpoly(zinv*zinv, (0.375, 0.1025390625, 0.2775764465332031, 1.993531733751297, 27.248827311268542, 603.8440767050702, 19718.37591223663))
-        r = e / s * (p - p2) - im * (p + p2) / (e * s)
+        r = e * sinv * (p - p2) - im * sinv * (p + p2) / e
     else
         zz = z*z
-        zz2 = zz*zz
-        r = evalpoly(zz2, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45))
-        r += zz*evalpoly(zz2, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48))
+        z4 = zz*zz
+        r = evalpoly(z4, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45))
+        r += zz*evalpoly(z4, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48))
         r *= z
     end
-    check && (r = conj(r))
+    isconj && (r = conj(r))
     return ComplexF32(r*c)
 end
 

test/besseli_test.jl

@@ -68,7 +68,7 @@ i1x_32 = besseli1x.(Float32.(x32))
 # test complex implementation
 
 for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 15.0, 17.0, 18.0, 20.0, 25.0, 50.0, 150.0], a in 0:pi/12:2pi
-    z = x*exp(im*a)
+    z = x*cis(a)
     @test isapprox(besseli0(z), SpecialFunctions.besseli(0, z), rtol=2e-14)
     @test isapprox(besseli0(ComplexF32(z)), ComplexF32(SpecialFunctions.besseli(0, ComplexF32(z))), rtol=1e-7)
     @test isapprox(besseli1(z), SpecialFunctions.besseli(1, z), rtol=2e-14)

test/besselj_test.jl

@@ -80,7 +80,7 @@ x = [1e12, 5e12, 1e13, 5e13, 1e14, 5e14, 1e15, 5e15, 1e16, 5e16, 1e17, 5e17, 1e1
 # test complex implementation
 
 for x in [0.0, 0.01, 0.5, 1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 15.0, 17.0, 18.0, 20.0, 25.0, 50.0, 150.0], a in 0:pi/12:2pi
-    z = x*exp(im*a)
+    z = x*cis(a)
     @test isapprox(besselj0(z), SpecialFunctions.besselj(0, z), rtol=2e-14)
     @test isapprox(besselj0(ComplexF32(z)), ComplexF32(SpecialFunctions.besselj(0, ComplexF32(z))), rtol=1e-7)
     @test isapprox(besselj1(z), SpecialFunctions.besselj(1, z), rtol=3e-14)