add in/de-crementαβ!

and toline! and fromline!

Author: MikaelSlevinsky

src/cjt.jl

@@ -9,10 +9,7 @@ function cjt(c::AbstractVector,plan::ChebyshevJacobiPlan)
         elseif α == 0.5 && β == -0.5
             decrementα!(ret,α,β)
         elseif α == 0.5 && β == 0.5
-            N > 1 && (ret[N] *= (4N-2)/N)
-            N > 2 && (ret[N-1] *= (4N-6)/(N-1))
-            @inbounds for i=N-2:-1:2 ret[i] = (4i-2)/i*(ret[i]+(2i+1)/(8i+8)*ret[i+2]) end
-            N > 2 && (ret[1] += 0.1875ret[3])
+            decrementαβ!(ret,α,β)
         end
         for i=1:N ret[i] *= Cx(i-1.0)/sqrtpi end
         return ret
@@ -30,17 +27,13 @@ function cjt(c::AbstractVector,plan::ChebyshevUltrasphericalPlan)
     N ≤ 1 && return c
     if λ == 0 || λ == 1
         ret = copy(c)
-        if λ == 1
-            N > 1 && (ret[N] *= (4N-2)/N)
-            N > 2 && (ret[N-1] *= (4N-6)/(N-1))
-            @inbounds for i=N-2:-1:2 ret[i] = (4i-2)/i*(ret[i]+(2i+1)/(8i+8)*ret[i+2]) end
-            N > 2 && (ret[1] += 0.1875ret[3])
-        end
+        λ == 1 && decrementαβ!(ret,λ-one(λ)/2,λ-one(λ)/2)
         for i=1:N ret[i] *= Cx(i-1.0)/sqrtpi end
         return ret
     else
         # Ultraspherical line
-        ret = ultra2cheb(c,λ,plan)
+        ret = toline!(copy(c),λ-one(λ)/2,λ-one(λ)/2)
+        ret = ultra2cheb(ret,modλ(λ),plan)
         return ret
     end
 end
@@ -59,10 +52,7 @@ function icjt(c::AbstractVector,plan::ChebyshevJacobiPlan)
             incrementα!(ret,α-1,β)
             return ret
         elseif α == 0.5 && β == 0.5
-            N > 2 && (ret[1] -= 0.1875ret[3])
-            @inbounds for i=2:N-2 ret[i] = i/(4i-2)*ret[i] - (2i+1)/(8i+8)*ret[i+2] end
-            N > 2 && (ret[N-1] *= (N-1)/(4N-6))
-            N > 1 && (ret[N] *= N/(4N-2))
+            incrementαβ!(ret,α-1,β-1)
             return ret
         else
             return ret
@@ -82,18 +72,12 @@ function icjt(c::AbstractVector,plan::ChebyshevUltrasphericalPlan)
     if λ == 0 || λ == 1
         ret = copy(c)
         for i=1:N ret[i] *= sqrtpi/Cx(i-1.0) end
-        if λ == 1
-            N > 2 && (ret[1] -= 0.1875ret[3])
-            @inbounds for i=2:N-2 ret[i] = i/(4i-2)*ret[i] - (2i+1)/(8i+8)*ret[i+2] end
-            N > 2 && (ret[N-1] *= (N-1)/(4N-6))
-            N > 1 && (ret[N] *= N/(4N-2))
-            return ret
-        else
-            return ret
-        end
+        λ == 1 && incrementαβ!(ret,λ-3one(λ)/2,λ-3one(λ)/2)
+        return ret
     else
         # Ultraspherical line
-        ret = cheb2ultra(c,λ,plan)
+        ret = cheb2ultra(c,modλ(λ),plan)
+        fromline!(ret,λ-one(λ)/2,λ-one(λ)/2)
         return ret
     end
 end
@@ -114,12 +98,12 @@ function plan_icjt(c::AbstractVector,α,β;M::Int=7)
 end
 
 function plan_cjt(c::AbstractVector,λ;M::Int=7)
-    P = ForwardChebyshevUltrasphericalPlan(c,λ,M)
+    P = ForwardChebyshevUltrasphericalPlan(c,modλ(λ),M)
     P.CUC.λ = λ
     P
 end
 function plan_icjt(c::AbstractVector,λ;M::Int=7)
-    P = BackwardChebyshevUltrasphericalPlan(c,λ,M)
+    P = BackwardChebyshevUltrasphericalPlan(c,modλ(λ),M)
     P.CUC.λ = λ
     P
 end

src/specialfunctions.jl

@@ -405,16 +405,17 @@ function incrementβ!(c::AbstractVector,α,β)
     N > 1 && (c[N] *= (αβ+N)/(αβ+2N-1))
     c
 end
-#=
+
 function incrementαβ!(c::AbstractVector,α,β)
     @assert α == β
     N = length(c)
-    @inbounds for i=1:N-2 c[i] = (2α+i-2)*(2α+i-1)/(2α+2i-1)/(2α+2i)*c[i] - (α+i+1)/(4α+4i+2)*c[i+2] end
-    N > 1 && (c[N-1] *= (2α+N-3)*(2α+N-2)/(2α+2N-3)/(2α+2N-2))
-    N > 2 && (c[N] *= (2α+N-2)*(2α+N-1)/(2α+2N-1)/(2α+2N))
+    N > 2 && (c[1] -= (α+2)/(4α+10)*c[3])
+    @inbounds for i=2:N-2 c[i] = (2α+i)*(2α+i+1)/(2α+2i-1)/(2α+2i)*c[i] - (α+i+1)/(4α+4i+6)*c[i+2] end
+    N > 1 && (c[N-1] *= (2α+N-1)*(2α+N)/(2α+2N-3)/(2α+2N-2))
+    N > 0 && (c[N] *= (2α+N)*(2α+N+1)/(2α+2N-1)/(2α+2N))
     c
 end
-=#
+
 function decrementα!(c::AbstractVector,α,β)
     αβ,N = α+β,length(c)
     N > 1 && (c[N] *= (αβ+2N-2)/(αβ+N-1))
@@ -430,16 +431,17 @@ function decrementβ!(c::AbstractVector,α,β)
     N > 1 && (c[1] -= (α+1)/(αβ+2)*c[2])
     c
 end
-#=
+
 function decrementαβ!(c::AbstractVector,α,β)
     @assert α == β
     N = length(c)
-    N > 1 && (c[N-1] *= (2α+2N-5)*(2α+2N-4)/(2α+N-5)/(2α+N-4))
-    N > 2 && (c[N] *= (2α+2N-3)*(2α+2N-2)/(2α+N-4)/(2α+N-3))
-    @inbounds for i=N-2:-1:1 c[i] = (2α+2i-3)*(2α+2i-2)/(2α+i-4)/(2α+i-3)*(c[i] + (α+i)/(4α+4i-2)*c[i+2])  end
+    N > 0 && (c[N] *= (2α+2N-3)*(2α+2N-2)/(2α+N-2)/(2α+N-1))
+    N > 1 && (c[N-1] *= (2α+2N-5)*(2α+2N-4)/(2α+N-3)/(2α+N-2))
+    @inbounds for i=N-2:-1:2 c[i] = (2α+2i-3)*(2α+2i-2)/(2α+i-2)/(2α+i-1)*(c[i] + (α+i)/(4α+4i+2)*c[i+2]) end
+    N > 2 && (c[1] += (α+1)/(4α+6)*c[3])
     c
 end
-=#
+
 
 function modαβ(α)
     if -0.5 < α ≤ 0.5
@@ -452,6 +454,16 @@ function modαβ(α)
     a
 end
 
+function modλ(λ)
+    if 0 ≤ λ < 1
+        l = λ
+    else
+        l = λ%1
+        l < 0 && (l+=1)
+    end
+    l
+end
+
 function tosquare!(ret::AbstractVector,α,β)
     a,b = modαβ(α),modαβ(β)
     A,B = α-a,β-b
@@ -505,3 +517,32 @@ function fromsquare!(ret::AbstractVector,α,β)
     end
     ret
 end
+
+
+function toline!(ret::AbstractVector,α,β)
+    @assert α == β
+    a,b = modαβ(α),modαβ(β)
+    A,B = α-a,β-b
+    if α ≤ -0.5
+        incrementαβ!(ret,α,β)
+    else
+        for i=A:-1:1
+            decrementαβ!(ret,i+a,i+a)
+        end
+    end
+    ret
+end
+
+function fromline!(ret::AbstractVector,α,β)
+    @assert α == β
+    a,b = modαβ(α),modαβ(β)
+    A,B = α-a,β-b
+    if α ≤ -0.5
+        decrementαβ!(ret,a,b)
+    else
+        for i=0:A-1
+            incrementαβ!(ret,i+a,i+a)
+        end
+    end
+    ret
+end