change file structure

Author: heltonmc

src/Bessels.jl

@@ -4,13 +4,42 @@ export besselj0
 export besselj1
 export bessely0
 export bessely1
+export besseli0
+export besseli0x
+export besseli1
+export besseli1x
+export besselk0
+export besselk0x
+export besselk1
+export besselk1x
 
+
+include("Float32/besseli.jl")
+include("Float32/besselj.jl")
+include("Float32/besselk.jl")
+include("Float32/bessely.jl")
+include("Float32/constants.jl")
+
+include("Float64/besseli.jl")
+include("Float64/besselj.jl")
+include("Float64/besselk.jl")
+include("Float64/bessely.jl")
+include("Float64/constants.jl")
+
+include("Float128/besseli.jl")
+include("Float128/besselj.jl")
+include("Float128/besselk.jl")
+include("Float128/bessely.jl")
+include("Float128/constants.jl")
+
+
+include("chebyshev.jl")
 include("constants.jl")
-include("j0_y0_constants.jl")
-include("j0.jl")
-include("y0.jl")
-include("j1.jl")
+#include("j0_y0_constants.jl")
+#include("j0.jl")
+#include("y0.jl")
+#include("j1.jl")
 
-include("parse.jl")
+#include("parse.jl")
 
 end

src/Float128/besseli.jl

@@ -0,0 +1,69 @@
+function besselj0(x::BigFloat)
+    x = abs(x)
+    T = eltype(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 2.0
+        z = x * x
+        p = z * z * evalpoly(z, J0_2N) / evalpoly(z, J0_2D)
+        p -= 0.25 * z
+        p += one(T)
+        return p
+    else
+        xinv = inv(x)
+        z  = xinv * xinv
+        if xinv <= 0.25
+            if xinv <= 0.125
+                if xinv <= 0.0625
+                    p = evalpoly(z, P16_IN) / evalpoly(z, P16_ID)
+                    q = evalpoly(z, Q16_IN) / evalpoly(z, Q16_ID)
+                else
+                    p = evalpoly(z, P8_16N) / evalpoly(z, P8_16D)
+                    q = evalpoly(z, Q8_16N) /  evalpoly(z, Q8_16D)
+                end
+            elseif xinv <= 0.1875
+                p = evalpoly(z, P5_8N) / evalpoly(z, P5_8D)
+	            q = evalpoly(z, Q5_8N) / evalpoly(z, Q5_8D)
+            else
+                p = evalpoly(z, P4_5N) / evalpoly(z, P4_5D)
+	            q = evalpoly(z, Q4_5N) / evalpoly(z, Q4_5D)
+            end
+        elseif xinv <= 0.5
+            if xinv <= 0.375
+                if xinv <= 0.3125
+                    p = evalpoly(z, P3r2_4N) / evalpoly(z, P3r2_4D)
+	                q = evalpoly(z, Q3r2_4N) / evalpoly(z, Q3r2_4D)
+                else
+                    p = evalpoly(z, P2r7_3r2N) / evalpoly(z, P2r7_3r2D)
+                    q = evalpoly(z, Q2r7_3r2N) / evalpoly(z, Q2r7_3r2D)
+                end
+            elseif xinv <= 0.4375
+                p = evalpoly(z, P2r3_2r7N) / evalpoly(z, P2r3_2r7D)
+	            q = evalpoly(z, Q2r3_2r7N) / evalpoly(z, Q2r3_2r7D)
+            else
+                p = evalpoly(z, P2_2r3N) / evalpoly(z, P2_2r3D)
+	            q = evalpoly(z, Q2_2r3N) / evalpoly(z, Q2_2r3D)
+            end
+        end
+        
+        p = 1.0 + z * p
+        q = z * xinv * q
+        q = q - 0.125 * xinv
+        c = cos(x)
+        s = sin(x)
+        ss = s - c
+        cc = s + c
+        z = -cos(x + x)
+        if (s * c) < 0
+            cc = z / ss;
+        else
+            ss = z / cc;
+        end
+        z = ONEOSQPI(T) * (p * cc - q * ss) / sqrt(x)
+        return z
+    end
+end

src/Float128/besselj.jl

@@ -0,0 +1,72 @@
+function bessely0(x::BigFloat)
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
+        end
+    elseif isinf(x)
+        return zero(x)
+    end
+  
+    xx = abs(x)
+
+    if xx <= 2.0
+        z = xx * xx
+        p = evalpoly(z, Y0_2N) / evalpoly(z, Y0_2D)
+        p = TWOOPI(BigFloat) * log(x) * besselj0(x) + p
+        return p
+    end
+
+    xinv = inv(xx)
+    z = xinv * xinv
+    if xinv <= 0.25
+        if xinv <= 0.125
+            if xinv <= 0.0625
+                p = evalpoly(z, P16_IN) / evalpoly(z, P16_ID)
+                q = evalpoly(z, Q16_IN) / evalpoly(z, Q16_ID)
+            else
+                p = evalpoly(z, P8_16N) / evalpoly(z, P8_16D)
+                q = evalpoly(z, Q8_16N) / evalpoly(z, Q8_16D)
+            end
+        elseif xinv <= 0.1875
+            p = evalpoly(z, P5_8N) / evalpoly(z, P5_8D)
+	        q = evalpoly(z, Q5_8N) / evalpoly(z, Q5_8D)
+        else
+            p = evalpoly(z, P4_5N) / evalpoly(z, P4_5D)
+	        q = evalpoly(z, Q4_5N) / evalpoly(z, Q4_5D)
+        end
+    else
+        if xinv <= 0.375
+            if xinv <= 0.3125
+                p = evalpoly(z, P3r2_4N) / evalpoly(z, P3r2_4D)
+	            q = evalpoly(z, Q3r2_4N) / evalpoly(z, Q3r2_4D)
+            else 
+                p = evalpoly(z, P2r7_3r2N) / evalpoly(z, P2r7_3r2D)
+	            q = evalpoly(z, Q2r7_3r2N) / evalpoly(z, Q2r7_3r2D)
+            end
+        elseif xinv <= 0.4375
+            p = evalpoly(z, P2r3_2r7N) / evalpoly(z, P2r3_2r7D)
+	        q = evalpoly(z, Q2r3_2r7N) / evalpoly(z, Q2r3_2r7D)
+        else
+            p = evalpoly(z, P2_2r3N) / evalpoly(z, P2_2r3D)
+	        q = evalpoly(z, Q2_2r3N) / evalpoly(z, Q2_2r3D)
+        end
+    end
+
+    p = 1.0 + z * p
+    q = z * xinv * q
+    q = q - 0.125 * xinv
+    c = cos(x)
+    s = sin(x)
+    ss = s - c
+    cc = s + c
+    z = -cos(x + x)
+    if s * c < zero(x)
+        cc = z / ss
+    else 
+        ss = z / cc
+    end
+    z = ONEOSQPI(BigFloat) * (p * ss + q * cc) / sqrt(x)
+    return z
+end

src/Float128/besselk.jl

@@ -0,0 +1,52 @@
+function besseli0(x::Float32)
+    T = Float32
+    if x < zero(x)
+        x = -x
+    end
+    if x <= 8.0f0
+        y = 0.5f0 * x - 2.0f0
+        return exp(x) * chbevl(y, A_i0(T))
+    else
+        return exp(x) * chbevl(32.0f0 / x - 2.0f0, B_i0(T)) / sqrt(x)
+    end
+end
+function besseli0x(x::Float32)
+    T = Float32
+    if x < zero(x)
+        x = -x
+    end
+    if x <= 8.0f0
+        y = 0.5f0 * x - 2.0f0
+        return chbevl(y, A_i0(T))
+    else
+        return chbevl(32.0f0 / x - 2.0f0, B_i0(T)) / sqrt(x)
+    end
+end
+function besseli1(x::Float32)
+    T = Float32
+    z = abs(x)
+    if z <= 8.0f0
+        y = 0.5f0 * z - 2.0f0
+        z = chbevl(y, A_i1(T)) * z * exp(z)
+    else
+        z = exp(z) * chbevl(32.0f0 / z - 2.0f0, B_i1(T)) / sqrt(z)
+    end
+    if x < zero(x)
+        z = -z
+    end
+    return z
+end
+function besseli1x(x::Float32)
+    T = Float32
+    z = abs(x)
+    if z <= 8.0f0
+        y = 0.5f0 * z - 2.0f0
+        z = chbevl(y, A_i1(T)) * z
+    else
+        z = chbevl(32.0f0 / z - 2.0f0, B_i1(T)) / sqrt(z)
+    end
+    if x < zero(x)
+        z = -z
+    end
+    return z
+end

src/Float128/bessely.jl

@@ -0,0 +1,51 @@
+function besselj0(x::Float32)
+    T = Float32
+    x = abs(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 2.0f0
+        z = x * x
+        if x < 1.0f-3
+            return 1.0f0 - 0.25f0 * z
+        end
+        DR1 = 5.78318596294678452118f0
+        p = (z - DR1) * evalpoly(z, JP_j0(T))
+        return p
+    else
+        q = inv(x)
+        w = sqrt(q)
+        p = w * evalpoly(q, MO_j0(T))
+        w = q * q
+        xn = q * evalpoly(w, PH_j0(T)) - PIO4(Float32)
+        p = p * cos(xn + x)
+        return p
+    end
+end
+
+function besselj1(x::Float32)
+    x = abs(x)
+    if iszero(x)
+        return zero(x)
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 2.0f0
+        z = x * x
+        Z1 = 1.46819706421238932572f1
+        p = (z - Z1) * x * evalpoly(z, JP32)
+        return p
+    else
+        q = inv(x)
+        w = sqrt(q)
+        p = w * evalpoly(q, MO132)
+        w = q * q
+        xn = q * evalpoly(w, PH132) - THPIO4(Float32)
+        p = p * cos(xn + x)
+        return p
+    end
+end

src/j0_y0_constants.jl renamed to src/Float128/constants.jl

@@ -0,0 +1,68 @@
+function besselk0(x::Float32)
+    T = Float32
+    if x <= zero(x)
+        return throw(DomainError(x, "NaN result for non-NaN input."))
+    end
+
+    if x <= 2.0f0
+        y = x * x - 2.0f0
+        y = chbevl(y, A_k0(T)) - log(0.5f0 * x) * besseli0(x)
+        return y
+    else
+        z = 8.0f0 / x - 2.0f0
+        y = exp(-x) * chbevl(z, B_k0(T)) / sqrt(x)
+        return y
+    end
+end
+function besselk0x(x::Float32)
+    T = Float32
+    if x <= zero(x)
+        return throw(DomainError(x, "NaN result for non-NaN input."))
+    end
+
+    if x <= 2.0f0
+        y = x * x - 2.0f0
+        y = chbevl(y, A_k0(T)) - log(0.5f0 * x) * besseli0(x)
+        return y * exp(x)
+    else
+        z = 8.0f0 / x - 2.0f0
+        y = chbevl(z, B_k0(T)) / sqrt(x)
+        return y
+    end
+end
+
+function besselk1(x::Float32)
+    T = Float32
+    if x <= zero(x)
+        if iszero(x)
+            return Inf32
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
+        end
+    end
+    if x <= 2.0f0
+        y = x * x - 2.0f0
+        y = log(0.5f0 * x) * besseli1(x) + chbevl(y, A_k1(T)) / x
+        return y
+    else
+        return exp(-x) * chbevl(8.0f0 / x - 2.0f0, B_k1(T)) / sqrt(x)
+    end
+end
+
+function besselk1x(x::Float32)
+    T = Float32
+    if x <= zero(x)
+        if iszero(x)
+            return Inf32
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
+        end
+    end
+    if x <= 2.0f0
+        y = x * x - 2.0f0
+        y = log(0.5f0 * x) * besseli1(x) + chbevl(y, A_k1(T)) / x
+        return y * exp(x)
+    else
+        return chbevl(8.0f0 / x - 2.0f0, B_k1(T)) / sqrt(x)
+    end
+end