SciFracX · ErikQQY · Dec 23, 2023 · Dec 23, 2023 · Dec 23, 2023 · Dec 23, 2023
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -10,8 +10,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
-          - 'nightly'
+          - '1'
+          - '~1.10.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

diff --git a/Project.toml b/Project.toml
@@ -2,7 +2,7 @@ name = "FractionalCalculus"
 uuid = "638fb199-4bb2-4014-80c8-6dc0d90f156b"
 license = "MIT"
 authors = ["Qingyu Qu <[email protected]>"]
-version = "0.2.10"
+version = "0.3.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -15,13 +15,13 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-InvertedIndices = "1.1"
-QuadGK = "2.4.2"
-SpecialFunctions = "1.6.2, 2"
-SpecialMatrices = "2.0.0"
-SymbolicUtils = "0.19.7"
-Symbolics = "4.4"
-julia = "1.0"
+InvertedIndices = "1.3"
+QuadGK = "2.9"
+SpecialFunctions = "2.3"
+SpecialMatrices = "3"
+SymbolicUtils = "1.5"
+Symbolics = "5.13"
+julia = "1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/docs/make.jl b/docs/make.jl
@@ -6,7 +6,8 @@ makedocs(
     format = Documenter.HTML(
         assets = ["assets/favicon.ico"],
     ),
-
+    warnonly = [:missing_docs, :cross_references],
+    doctest = false,
     pages = [
         "index.md",
         "Fractional Derivative" => [

diff --git a/docs/src/Derivative/derivative.md b/docs/src/Derivative/derivative.md
@@ -4,7 +4,7 @@
 Pages = ["derivative.md"]
 ```
 
-To get start with fractional derivative, you need to know that unlike Newtonian derivatives, fractional derivative is defined via integral.
+To get started with fractional derivatives, you need to know that unlike Newtonian derivatives, fractional derivative is defined via integral.
 
 !!! tip "Non-Local Operators"
 	It is noteworthy that the fractional derivatives are not local operators, which means that we cannot calculate the fractional derivative solely on the basis of function values of $f(x)$ taken from neighborhood of the point $x$. Instead, we have to take the entire history of $f(x)$ (i.e., all function values $f(x)$ for $0<x<a$) into account.
@@ -30,13 +30,13 @@ julia> fracdiff(x->x, 0.5, 1, 0.0001, RLDiffL1())
 
 ## Caputo sense derivative
 
-There many types of definitions of fractional derivative, Caputo is one of these useful definitions. The Caputo fractional derivative is first be proposed in [Michele Caputo's Paper](https://doi.org/10.1111/j.1365-246X.1967.tb02303.x), 
+There are many types of definitions of fractional derivatives, Caputo is one of these useful definitions. The Caputo fractional derivative is first proposed in [Michele Caputo's Paper](https://doi.org/10.1111/j.1365-246X.1967.tb02303.x), 
 
 ```math
 ^CD_t^\alpha f(t) = \frac{1}{\Gamma(n-\alpha)}\int_0^t\frac{f^{(n)}(\tau)d\tau}{(t-\tau)^{\alpha+1-n}}, n=\lceil{\alpha}\rceil
 ```
 
-In **FractionalCalculus.jl**, let's see, if you want to calculate the $0.5$ order fractional derivative of $f(x)=x$ at a $x=1$ with step size $0.0001$, simply typing these:
+In **FractionalCalculus.jl**, let's see, if you want to calculate the $0.5$ order fractional derivative of $f(x)=x$ at a $x=1$ with step size $0.0001$, simply type these:
 
 
 ```julia-repl
@@ -85,10 +85,10 @@ julia> fracdiff(x->x, 0.5, collect(0:0.01:1), 2, GLHighPrecision())
  1.1283791670955126
 ```
 
-Here, we use the high precision algorithm, the fourth parameter means we set the precision order as **p=2**. The returned result means the derivative on the interval $[0, 1]$.
+Here, we use the high-precision algorithm, the fourth parameter means we set the precision order as **p=2**. The returned result means the derivative on the interval $[0, 1]$.
 
 !!! info
-	If the function ``f(t)`` is suitably smooth, then the Grünwald Letnikov sense derivative and the Riemann Liouville sense derivative is equivalent.
+	If the function ``f(t)`` is suitably smooth, then the Grünwald Letnikov sense derivative and the Riemann Liouville sense derivative are equivalent.
 
 ## Riesz sense derivative
 
@@ -179,4 +179,4 @@ julia> fracdiff(x->x, 0.5, 1, 0.01, AtanganaSeda())
 !!! note
 	Here we need to specify the **start point** and **end point**
 
-There are different approximating methods being used in the computing, choose the one you need😉
+There are different approximating methods being used in computing, choose the one you need😉
diff --git a/docs/src/Derivative/derivativeapi.md b/docs/src/Derivative/derivativeapi.md
@@ -1,8 +1,8 @@
 # Fractional derivative API
 
-FractionalCalculus.jl has a disciplinary type system, by specifying which type you want to use, you can use the relating algorithm to compute fractional differentiation and fractional integral. 
+FractionalCalculus.jl has a disciplinary type system, by specifying which type you want to use, you can use the related algorithm to compute fractional differentiation and fractional integral. 
 
-This page contains all of the existing API we can use for computing fractional derivative.
+This page contains all of the existing APIs we can use for computing fractional derivatives.
 
 ```@docs
 FractionalCalculus.fracdiff
@@ -30,10 +30,13 @@ FractionalCalculus.GLLagrangeThreePointInterp
 ## Riemann Liouville sense fractional derivative
 
 ```@docs
-FractionalCalculus.L1
+FractionalCalculus.RLDiffL1
 FractionalCalculus.RLDiffMatrix
+FractionalCalculus.RLLinearSplineInterp
 FractionalCalculus.RLG1
 FractionalCalculus.RLD
+FractionalCalculus.RLDiffL2
+FractionalCalculus.RLDiffL2C
 ```
 
 ## Hadamard sense fractional derivative

diff --git a/docs/src/Derivative/short_memory_effect.md b/docs/src/Derivative/short_memory_effect.md
@@ -8,7 +8,7 @@ Which means in the interval $[t-L, t]$, $L$ is the "memory length".
 _aD^\alpha_t f(t)\approx _{t-L}D^\alpha_t f(t)
 ```
 
-By deploying the **Short Memory Effect**, we can reduce our numerical cost while retain the precision in a way.
+By deploying the **Short Memory Effect**, we can reduce our numerical cost while retaining the precision in a way.
 
 !!! info "Short memory effect in FDE"
-    Want to see how the short memory is used in fractional differential equations? Please see [short memory effect in FDE](http://scifracx.org/FractionalDiffEq.jl/dev/system_of_FDE/#Short-memory-effect-in-FDE)
+    Want to see how short memory is used in fractional differential equations? Please see [the short memory effect in FDE](http://scifracx.org/FractionalDiffEq.jl/dev/system_of_FDE/#Short-memory-effect-in-FDE)
diff --git a/docs/src/Integral/integral.md b/docs/src/Integral/integral.md
@@ -15,7 +15,7 @@ In FractionalCalculus, you can compute the integral of a function with order $\a
 julia> fracint(x->x, 0.5, 0, 1, 0.0001, RLDirect())
 ```
 
-A tuple contains result and estimating error will be returned.
+A tuple contains the result and estimating error will be returned.
 
 ```julia
 julia> fracint(x->x, 0.5, 0, 1, 0.0001, RLDirect())

diff --git a/docs/src/Integral/integralapi.md b/docs/src/Integral/integralapi.md
@@ -1,6 +1,6 @@
 # Fractional integral API
 
-FractionalCalculus.jl has a disciplinary type system, by specifying which type you want to use, you can use the relating algorithm to compute fractional differentiation and fractional integral. 
+FractionalCalculus.jl has a disciplinary type system, by specifying which type you want to use, you can use the related algorithm to compute fractional differentiation and fractional integral. 
 
 This page contains all of the existing API we can use for computing fractional integral.
 

diff --git a/docs/src/Integral/symintegral.md b/docs/src/Integral/symintegral.md
@@ -1,6 +1,6 @@
 # Symbolic Fractional Integration
 
-Similiar with symbolic fractional differentiation in FractionalCalculus.jl, we first define fractional integration rules to implement symbolic fractional integration computing.
+Similiar to symbolic fractional differentiation in FractionalCalculus.jl, we first define fractional integration rules to implement symbolic fractional integration computing.
 
 All we have to do is to call the ```semiint``` function:
 

diff --git a/src/Derivative/RL.jl b/src/Derivative/RL.jl
@@ -96,9 +96,14 @@
 """
 struct RLD <: RLDiff end
 
-
+"""
+RLDiffL2
+"""
 struct RLDiffL2 <: RLDiff end
 
+"""
+RLDiffL2C
+"""
 struct RLDiffL2C <: RLDiff end
 ################################################################
 ###                    Type definition done                  ###
@@ -236,7 +241,7 @@
    return h^(-α)/gamma(-α)*result
 end

 fracdiff(f::FunctionAndNumber, α, end_point, h::Float64, ::RLG1) = fracdiff(f::FunctionAndNumber, α, 0, end_point, h::Float64, RLG1())

 function fracdiff(f::FunctionAndNumber,
                  α::Float64,
@@ -269,34 +274,34 @@
    return n^α/(α*(1-α))*temp
 end

 function fracdiff(f::FunctionAndNumber,
                  α::Float64,
                  point::Real,
                  h::Float64,
                  ::RLDiffL2C)
    N = round(Int, point/h)
    temp = zero(Float64)
    for k=-1:N+1
        temp += Ŵ(k, α, N)*f((N-k)*h)
    end
 end

 function Ŵ(k::Int64, α::Float64, N::Int64)
    expo = 2-α
    if k == -1
        return 1
    elseif k == 0
        return 2^expo-2
    elseif k == 1
        return 3^expo-2^expo
    elseif 2 ≤ k ≤ N-2
        return (k+2)^expo - 2*(k+1)^expo + 2*(k-1)^expo - (k-2)^expo
    elseif k == N-1
        return -N^expo-(N-3)^expo + 2*(N-2)^expo
    elseif k == N
        return -N^expo + 2*(N-1)^expo-(N-2)^expo
    elseif k == N+1
        return N^expo - (N-1)^expo
    end
 end
 #=
@@ -306,30 +311,30 @@
 =#


 function fracdiff(f::FunctionAndNumber,
                  alpha::Float64,
                  point::Real,
                  h::Float64,
                  ::RLDiffL2)

    n = round(Int64, point/h)
    result = zero(Float64)
    for k=-1:n
        result += WK(k, alpha, n)*f((n-k)*h)
    end

    return f(0)*point^(-alpha)/gamma(1-alpha) + derivative(f, 0)*point^(1-alpha)/gamma(2-alpha) + h^(-alpha)/gamma(3-alpha)*result
 end

 function WK(k::Int64, alpha::Float64, n::Int64)
    if k == -1
        return 1
    elseif k == 0
        return 2^(2-alpha)-3
    elseif 1 ≤ k ≤ n-2
        return (k+2)^(2-alpha) - 3*(k+1)^(2-alpha) + 3*k^(2-alpha) - (k-1)^(2-alpha)
    elseif k == n-1
        return -2*n^(2-alpha) + 3*(n-1)^(2-alpha) - (n-2)^(2-alpha)
    elseif k == n
        return n^(2-alpha) - (n-1)^(2-alpha)
    end

diff --git a/src/Derivative/SymbolicDiff.jl b/src/Derivative/SymbolicDiff.jl
@@ -22,92 +22,92 @@
 
 SEMIDIFFRULES = [
 
-        @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys));
+        @rule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys));
         # CONSTANTS AND POWERS
-        @acrule sqrt(~x) => sqrt(pi)/2
-        #@acrule((~x) => 2*sqrt(~x/pi))
-        @acrule (~x)^(~α) => gamma(~α+1)/gamma(~α+1/2)*(~x)^(~α-1/2)
+        #@rule sqrt(~x) => sqrt(pi)/2
+        #@rule((~x) => 2*sqrt(~x/pi))
+        @rule (~x)^(~α) => gamma(~α+1)/gamma(~α+1/2)*(~x)^(~α-1/2)
 
         # BINOMIALS
-        @acrule 1/sqrt(~x) => 0
-        @acrule sqrt(1+~x) => 1/sqrt(pi*~x) + atan(sqrt(~x))/sqrt(pi)
-        @acrule sqrt(1-~x) => 1/sqrt(pi*~x) - atanh(sqrt(~x))/sqrt(pi)
-        @acrule 1/sqrt(1+~x) => 1/(sqrt(pi*~x)*(1+~x))
-        @acrule 1/sqrt(1-~x) => 1/(sqrt(pi*~x)*(1-~x))
-        @acrule 1/(1+~x) => (sqrt(1+~x)-sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1+~x)^(3/2))
-        @acrule 1/(1-~x) => (sqrt(1-~x)+sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1-~x)^(3/2))
-        @acrule 1/(1+~x)^(3/2) => (1-~x)/(sqrt(pi*~x)*(1+~x)^2)
-        @acrule 1/(1-~x)^(3/2) => (1+~x)/(sqrt(pi*~x)*(1-~x)^2)
-        @acrule 1/(1-~x)^(~p) => -Incomplete_beta(~x, -1/2, ~p+1/2)/(2*sqrt(pi)*(1-~x)^(~p+1/2))
-        @acrule 1/sqrt(~x*(1-~x)) => (ellipe(~x)-(1-~x)*ellipk(~x))/(~x*(1-~x)*sqrt(pi))
-        @acrule 1/sqrt(~x*(1+~x)) => (ellipe(~x/(1+~x))-ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
-        @acrule sqrt(~x/(1-~x)) => ellipe(~x)/((1-~x)*sqrt(pi))
+        #@rule 1/sqrt(~x) => 0
+        @rule sqrt(1+~x) => 1/sqrt(pi*~x) + atan(sqrt(~x))/sqrt(pi)
+        @rule sqrt(1-~x) => 1/sqrt(pi*~x) - atanh(sqrt(~x))/sqrt(pi)
+        @rule 1/sqrt(1+~x) => 1/(sqrt(pi*~x)*(1+~x))
+        @rule 1/sqrt(1-~x) => 1/(sqrt(pi*~x)*(1-~x))
+        @rule 1/(1+~x) => (sqrt(1+~x)-sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1+~x)^(3/2))
+        @rule 1/(1-~x) => (sqrt(1-~x)+sqrt(~x)*asinh(sqrt(~x)))/(sqrt(pi*~x)*(1-~x)^(3/2))
+        @rule 1/(1+~x)^(3/2) => (1-~x)/(sqrt(pi*~x)*(1+~x)^2)
+        @rule 1/(1-~x)^(3/2) => (1+~x)/(sqrt(pi*~x)*(1-~x)^2)
+        @rule 1/(1-~x)^(~p) => -Incomplete_beta(~x, -1/2, ~p+1/2)/(2*sqrt(pi)*(1-~x)^(~p+1/2))
+        @rule 1/sqrt(~x*(1-~x)) => (ellipe(~x)-(1-~x)*ellipk(~x))/(~x*(1-~x)*sqrt(pi))
+        @rule 1/sqrt(~x*(1+~x)) => (ellipe(~x/(1+~x))-ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
+        @rule sqrt(~x/(1-~x)) => ellipe(~x)/((1-~x)*sqrt(pi))
         # Doesn't support Legendre
-        @acrule 1/(sqrt(~x)*(1+~x)) => -sqrt(pi)/(2*(1+~x)^(3/2))
-        @acrule 1/(sqrt(~x)*(1-~x)) => sqrt(pi)/(2*(1-~x)^(3/2))
-        @acrule sqrt(~x)/~x => sqrt(pi)/(2*(1+~x)^(3/2))
-        @acrule sqrt(~x)/(-~x) => sqrt(pi)/(2*(1-~x)^(3/2))
-        @acrule sqrt(~x)/(1-~x)^2 => sqrt(pi)*(2+~x)/(4*(1-~x)^(5/2)) 
-        #@acrule((~x)^(~p)/(1-(~x))^((~p)+3/2) => (((~p)+1/2+1/2*(~x))*gamma(~p+1)*(~x)^((~p)-1/2))/(gamma((~p+3/2)*(1-(~x))^(2+(~p))))
-        @acrule sqrt((1-~x)/~x) => (ellipe(~x)-ellipk(~x))/(~x*sqrt(pi))
-        @acrule sqrt((1+~x)/~x) => ((1+~x)*ellipe(~x/(1+~x)) - ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
-        @acrule  (~x)^(~p)/(1+~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, ~x)
-        @acrule  (~x)^(~p)/(1-~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, -~x)
+        @rule 1/(sqrt(~x)*(1+~x)) => -sqrt(pi)/(2*(1+~x)^(3/2))
+        @rule 1/(sqrt(~x)*(1-~x)) => sqrt(pi)/(2*(1-~x)^(3/2))
+        @rule sqrt(~x)/~x => sqrt(pi)/(2*(1+~x)^(3/2))
+        @rule sqrt(~x)/(-~x) => sqrt(pi)/(2*(1-~x)^(3/2))
+        @rule sqrt(~x)/(1-~x)^2 => sqrt(pi)*(2+~x)/(4*(1-~x)^(5/2)) 
+        #@rule((~x)^(~p)/(1-(~x))^((~p)+3/2) => (((~p)+1/2+1/2*(~x))*gamma(~p+1)*(~x)^((~p)-1/2))/(gamma((~p+3/2)*(1-(~x))^(2+(~p))))
+        @rule sqrt((1-~x)/~x) => (ellipe(~x)-ellipk(~x))/(~x*sqrt(pi))
+        @rule sqrt((1+~x)/~x) => ((1+~x)*ellipe(~x/(1+~x)) - ellipk(~x/(1+~x)))/(~x*sqrt(pi*(1+~x)))
+        @rule  (~x)^(~p)/(1+~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, ~x)
+        @rule  (~x)^(~p)/(1-~x)^~r => gamma(~p+1)/gamma(~p+1/2)*(~x)^(~p-1/2)* Hypergeometric1F1(~r, ~p+1, p+1/2, -~x)
 
 
 
         # EXPONENTIAL AND RELATED FUNCTIONS
-        @acrule exp(~x) => 1/sqrt(pi*~x) + exp(~x)*erf(sqrt(~x))
-        @acrule exp(-~x) => 1/sqrt(pi*~x) - 2/sqrt(pi)*dawson(sqrt(-~x))
-        @acrule exp(~x)*erf(sqrt(~x)) => exp(~x)
-        @acrule dawson(sqrt(~x)) => 1/2*sqrt(pi)*exp(-~x)
-        @acrule exp(~x)*erfc(sqrt(~x)) => 1/sqrt(pi*~x) - exp(~x)*erfc(sqrt(~x))
-        @acrule exp(~x)*erfc(-sqrt(~x)) => 1/sqrt(pi*~x) + exp(~x)*erf(-sqrt(~x))
-        @acrule erf(~x) => exp(-1/2*~x)*besseli(0, 1/2*~x)
-        @acrule exp(~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(1/2*~x)*(besseli(1, 1/2*~x)+besseli(0, 1/2*~x))
-        @acrule exp(-~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(-1/2*~x)*(besseli(1, 1/2*~x)-besseli(0, 1/2*~x)) 
-        @acrule cosh(sqrt(~x))/sqrt(~x) => sqrt(pi/~x)*besseli(1, sqrt(~x))/2
-        @acrule  (1-cos(sqrt(~x)))/~x => sqrt(pi)/2*Struve0(sqrt(~x)) 
-        @acrule  (1-cosh(sqrt(~x)))/~x => sqrt(pi)/2*MStruve0(sqrt(~x)) 
-
-        @acrule sin(~x) => sin(~x+pi/4)-sqrt(2)*(AuxiliaryFresnelCos(sqrt(2*~x/π)))
-        @acrule cos(~x) => 1/(sqrt(π*~x))+cos(~x+1/4*π)-sqrt(2)*AuxiliaryFresnelSin(sqrt(2*~x/π))
-        @acrule sinh(~x) => dawson(sqrt(~x))/sqrt(pi)-(exp(~x)*erf(sqrt(~x)))/2
-        @acrule cosh(~x) => 1/sqrt(pi*~x)+exp(~x)*erf(sqrt(~x))/2-dawson(sqrt(~x))/sqrt(pi)
-        @acrule asin(sqrt(~x))/sqrt(1-~x) => sqrt(pi)/(2*(1-~x))
-        @acrule atan(sqrt(~x)) => 1/2*sqrt(pi/(1+~x))
-        @acrule asinh(sqrt(~x))/sqrt(1+~x) => sqrt(pi)/(2*(1+~x))
-        @acrule atanh(sqrt(~x)) => 1/2*sqrt(pi/(1-~x))
+        @rule exp(~x) => 1/sqrt(pi*~x) + exp(~x)*erf(sqrt(~x))
+        @rule exp(-~x) => 1/sqrt(pi*~x) - 2/sqrt(pi)*dawson(sqrt(-~x))
+        @rule exp(~x)*erf(sqrt(~x)) => exp(~x)
+        @rule dawson(sqrt(~x)) => 1/2*sqrt(pi)*exp(-~x)
+        @rule exp(~x)*erfc(sqrt(~x)) => 1/sqrt(pi*~x) - exp(~x)*erfc(sqrt(~x))
+        @rule exp(~x)*erfc(-sqrt(~x)) => 1/sqrt(pi*~x) + exp(~x)*erf(-sqrt(~x))
+        @rule erf(~x) => exp(-1/2*~x)*besseli(0, 1/2*~x)
+        @rule exp(~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(1/2*~x)*(besseli(1, 1/2*~x)+besseli(0, 1/2*~x))
+        @rule exp(-~x)/sqrt(~x) => 1/2*sqrt(pi)*exp(-1/2*~x)*(besseli(1, 1/2*~x)-besseli(0, 1/2*~x)) 
+        @rule cosh(sqrt(~x))/sqrt(~x) => sqrt(pi/~x)*besseli(1, sqrt(~x))/2
+        @rule  (1-cos(sqrt(~x)))/~x => sqrt(pi)/2*Struve0(sqrt(~x)) 
+        @rule  (1-cosh(sqrt(~x)))/~x => sqrt(pi)/2*MStruve0(sqrt(~x)) 
+
+        @rule sin(~x) => sin(~x+pi/4)-sqrt(2)*(AuxiliaryFresnelCos(sqrt(2*~x/π)))
+        @rule cos(~x) => 1/(sqrt(π*~x))+cos(~x+1/4*π)-sqrt(2)*AuxiliaryFresnelSin(sqrt(2*~x/π))
+        @rule sinh(~x) => dawson(sqrt(~x))/sqrt(pi)-(exp(~x)*erf(sqrt(~x)))/2
+        @rule cosh(~x) => 1/sqrt(pi*~x)+exp(~x)*erf(sqrt(~x))/2-dawson(sqrt(~x))/sqrt(pi)
+        @rule asin(sqrt(~x))/sqrt(1-~x) => sqrt(pi)/(2*(1-~x))
+        @rule atan(sqrt(~x)) => 1/2*sqrt(pi/(1+~x))
+        @rule asinh(sqrt(~x))/sqrt(1+~x) => sqrt(pi)/(2*(1+~x))
+        @rule atanh(sqrt(~x)) => 1/2*sqrt(pi/(1-~x))
 
         # BESSEL AND STRUVE FUNCTIONS
-        @acrule besselj(0, sqrt(~x)) => cos(sqrt(~x))/sqrt(pi*(~x))
-        @acrule besselj(1, sqrt(~x))/sqrt(~x) => (cos(sqrt(~x))+sqrt(~x)*sin(sqrt(~x))-1)/(sqrt(~x)*(~x)^(3/2))
-        @acrule  besselj(~v, sqrt(~x))/(~x)^(~v/2) => 1/(2^(~v)*gamma(~v+1)*sqrt(pi*~x))- Struve(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
-        @acrule sqrt(~x)*besselj(1, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi)
-        @acrule (~x)^(~v/2)*besselj(1, sqrt(~x)) => (~x)^((~v)/2-1/4)*besselj(~v-1/2, sqrt(~x))/sqrt(2)
-        @acrule besseli(0, sqrt(~x)) => cosh(sqrt(~x))/sqrt(pi*~x)
-        @acrule besseli(1, sqrt(~x))/sqrt(~x) => (cosh(sqrt(~x))-sqrt(~x)*sinh(sqrt(~x))-1)/(sqrt(pi)*(~x)^(3/2))
-        @acrule  besseli(~v, sqrt(~x))/((~x)^(~v/2)) => 1/(2^~v*gamma(~v+1)*sqrt(pi*(~x))) + Legendre(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
-        @acrule sqrt(~x)*besseli(1, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi)
-        @acrule  (~x)^(~v/2)*besseli(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*besseli(~v-1/2, sqrt(~x))/sqrt(2) 
-        @acrule exp(-~x)*besseli(0, ~x) => exp(-2*~x)/sqrt(pi*(~x))
-        @acrule  Struve(0, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi*(~x))
-        @acrule  Struve(0, sqrt(~x))/(~x) => (sin(sqrt(~x)) - SinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2)) 
-        @acrule sqrt(~x)*Struve(1, sqrt(~x)) =>(1-cos(sqrt(~x)))/sqrt(pi)
-        @acrule (~x)^(~v/2)*Struve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*Struve(~v-1/2, sqrt(~x))/sqrt(2)
-
-        @acrule MStruve(0, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi*(~x))
-        @acrule MStruve(0, sqrt(~x))/(~x) => (sinh(sqrt(~x)) - HyperSinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2))
-        @acrule sqrt(~x)*MStruve(1, sqrt(~x)) => (cosh(sqrt(~x)-1))/sqrt(pi)
-        @acrule  (~x)^(~v/2)*MStruve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*MStruve(~v-1/2, sqrt(~x))/sqrt(2)
+        @rule besselj(0, sqrt(~x)) => cos(sqrt(~x))/sqrt(pi*(~x))
+        @rule besselj(1, sqrt(~x))/sqrt(~x) => (cos(sqrt(~x))+sqrt(~x)*sin(sqrt(~x))-1)/(sqrt(~x)*(~x)^(3/2))
+        @rule  besselj(~v, sqrt(~x))/(~x)^(~v/2) => 1/(2^(~v)*gamma(~v+1)*sqrt(pi*~x))- Struve(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
+        @rule sqrt(~x)*besselj(1, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi)
+        @rule (~x)^(~v/2)*besselj(1, sqrt(~x)) => (~x)^((~v)/2-1/4)*besselj(~v-1/2, sqrt(~x))/sqrt(2)
+        @rule besseli(0, sqrt(~x)) => cosh(sqrt(~x))/sqrt(pi*~x)
+        @rule besseli(1, sqrt(~x))/sqrt(~x) => (cosh(sqrt(~x))-sqrt(~x)*sinh(sqrt(~x))-1)/(sqrt(pi)*(~x)^(3/2))
+        @rule  besseli(~v, sqrt(~x))/((~x)^(~v/2)) => 1/(2^~v*gamma(~v+1)*sqrt(pi*(~x))) + Legendre(~v+1/2, sqrt(~x))/(sqrt(2)*(~x)^(~v/2+1/4))
+        @rule sqrt(~x)*besseli(1, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi)
+        @rule  (~x)^(~v/2)*besseli(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*besseli(~v-1/2, sqrt(~x))/sqrt(2) 
+        @rule exp(-~x)*besseli(0, ~x) => exp(-2*~x)/sqrt(pi*(~x))
+        @rule  Struve(0, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi*(~x))
+        @rule  Struve(0, sqrt(~x))/(~x) => (sin(sqrt(~x)) - SinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2)) 
+        @rule sqrt(~x)*Struve(1, sqrt(~x)) =>(1-cos(sqrt(~x)))/sqrt(pi)
+        @rule (~x)^(~v/2)*Struve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*Struve(~v-1/2, sqrt(~x))/sqrt(2)
+
+        @rule MStruve(0, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi*(~x))
+        @rule MStruve(0, sqrt(~x))/(~x) => (sinh(sqrt(~x)) - HyperSinIntegral(sqrt(~x)))/(sqrt(pi)*(~x)^(3/2))
+        @rule sqrt(~x)*MStruve(1, sqrt(~x)) => (cosh(sqrt(~x)-1))/sqrt(pi)
+        @rule  (~x)^(~v/2)*MStruve(~v, sqrt(~x)) => (~x)^((~v)/2-1/4)*MStruve(~v-1/2, sqrt(~x))/sqrt(2)
 
         # MISCELLANEOUS FUNCTIONS
-        @acrule log(~x) => log(4*~x)/sqrt(pi*~x)
-        @acrule sqrt(~x)*log(~x) => sqrt(pi)/2*(log(~x/4)+2)
-        @acrule log(~x)/sqrt(~x) => sqrt(pi)/~x
+        @rule log(~x) => log(4*~x)/sqrt(pi*~x)
+        @rule sqrt(~x)*log(~x) => sqrt(pi)/2*(log(~x/4)+2)
+        @rule log(~x)/sqrt(~x) => sqrt(pi)/(~x)
 
-        @acrule ellipk(~x) => sqrt(pi)/(2*sqrt(~x*(1-~x)))
-        @acrule ellipe(~x) => 1/2*sqrt(pi*(1-~x)/(~x))
+        @rule ellipk(~x) => sqrt(pi)/(2*sqrt(~x*(1-~x)))
+        @rule ellipe(~x) => 1/2*sqrt(pi*(1-~x)/(~x))
 
 ]