remove restriction on same band limit for m as l

issue #16

Author: MikaelSlevinsky

src/SphericalHarmonics/SphericalHarmonics.jl

@@ -16,7 +16,7 @@ include("thinplan.jl")
 
 function plan_sph2fourier(A::AbstractMatrix; opts...)
     M, N = size(A)
-    if M ≤ 1022
+    if M ≤ 1023
         SlowSphericalHarmonicPlan(A)
     else
         ThinSphericalHarmonicPlan(A; opts...)

src/SphericalHarmonics/fastplan.jl

@@ -10,7 +10,6 @@ end
 
 function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
-    @assert N == 2M-1
     n = (N+1)÷2
     RP = RotationPlan(T, n-1)
     p1 = plan_normleg2cheb(A)

src/SphericalHarmonics/sphfunctions.jl

@@ -15,7 +15,7 @@ function sphrand{T}(::Type{T}, m::Int, n::Int)
     for i = 1:m
         A[i,1] = rand(T)
     end
-    for j = 1:n
+    for j = 1:n-1
         for i = 1:m-j
             A[i,2j] = rand(T)
             A[i,2j+1] = rand(T)
@@ -29,7 +29,7 @@ function sphrandn{T}(::Type{T}, m::Int, n::Int)
     for i = 1:m
         A[i,1] = randn(T)
     end
-    for j = 1:n
+    for j = 1:n-1
         for i = 1:m-j
             A[i,2j] = randn(T)
             A[i,2j+1] = randn(T)
@@ -43,7 +43,7 @@ function sphones{T}(::Type{T}, m::Int, n::Int)
     for i = 1:m
         A[i,1] = one(T)
     end
-    for j = 1:n
+    for j = 1:n-1
         for i = 1:m-j
             A[i,2j] = one(T)
             A[i,2j+1] = one(T)

src/SphericalHarmonics/thinplan.jl

@@ -14,7 +14,6 @@ end
 
 function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
-    @assert N == 2M-1
     n = (N+1)÷2
     RP = RotationPlan(T, n-1)
     p1 = plan_normleg2cheb(A)

test/sphericalharmonictests.jl

@@ -26,31 +26,31 @@ println()
 
 import FastTransforms: sphevaluatepi, normalizecolumns!, maxcolnorm
 
-n = 255
-A = sphrandn(Float64, n+1, n+1);
+n = 256
+A = sphrandn(Float64, n, n);
 normalizecolumns!(A);
 B = sph2fourier(A);
 
 for θ in (0.123, 0.456)
-    S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,1] for ℓ in 1:n+1)
-    SA = sum(sphevaluatepi(θ,ℓ-1,0)*A[ℓ,1] for ℓ in 1:n+1)
+    S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,1] for ℓ in 1:n)
+    SA = sum(sphevaluatepi(θ,ℓ-1,0)*A[ℓ,1] for ℓ in 1:n)
     @test norm(S0-SA) < 1000eps()
-    for m in 3:2:n+1
-        S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-2] for ℓ in 1:n+1)
-        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n+1)
+    for m in 3:2:n
+        S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-2] for ℓ in 1:n)
+        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n)
         @test norm(S0-SA) < 1000eps()
-        S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-1] for ℓ in 1:n+1)
-        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n+1)
+        S0 = sum(cospi((ℓ-1)*θ)*B[ℓ,2m-1] for ℓ in 1:n)
+        SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n)
         @test norm(S0-SA) < 1000eps()
     end
 end
 
-for θ in (0.123, 0.456), m in 2:2:n+1
-    S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-2] for ℓ in 1:n+1)
-    SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n+1)
+for θ in (0.123, 0.456), m in 2:2:n
+    S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-2] for ℓ in 1:n)
+    SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-2] for ℓ in 1:n)
     @test norm(S1-SA) < 1000eps()
-    S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-1] for ℓ in 1:n+1)
-    SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n+1)
+    S1 = sum(sinpi(ℓ*θ)*B[ℓ,2m-1] for ℓ in 1:n)
+    SA = sum(sphevaluatepi(θ,ℓ+m-2,m-1)*A[ℓ,2m-1] for ℓ in 1:n)
     @test norm(S1-SA) < 1000eps()
 end
 
@@ -60,8 +60,8 @@ println()
 
 import FastTransforms: normalizecolumns!, maxcolnorm
 
-n = VERSION < v"0.6.0-" ? 255 : 511
-A = sphrandn(Float64, n+1, n+1);
+n = VERSION < v"0.6.0-" ? 256 : 512
+A = sphrandn(Float64, n, n);
 normalizecolumns!(A);
 
 B = sph2fourier(A)
@@ -75,8 +75,8 @@ C = P\B
 println("The backward difference between slow plan and original: ", maxcolnorm(A-C))
 
 if VERSION ≥ v"0.6.0-"
-    n = 1023
-    A = sphrandn(Float64, n+1, n+1);
+    n = 1024
+    A = sphrandn(Float64, n, n);
     normalizecolumns!(A);
 
     B = sph2fourier(A; sketch = :none)

test/test_fastplan.jl

@@ -1,8 +1,8 @@
 import FastTransforms: allranks, normalizecolumns!, maxcolnorm
 
-n = 362
+m, n = 362, 300
 
-A = sphrandn(Float64, n+1, n+1);
+A = sphrandn(Float64, m, n);
 normalizecolumns!(A);
 B = zero(A);
 C = zero(A);

test/test_thinplan.jl

@@ -1,8 +1,8 @@
 import FastTransforms: allranks, normalizecolumns!, maxcolnorm
 
-n = VERSION < v"0.6.0-" ? 362 : 724
+m, n = VERSION < v"0.6.0-" ? (362, 300) : (724, 600)
 
-A = sphrandn(Float64, n+1, n+1);
+A = sphrandn(Float64, m, n);
 normalizecolumns!(A);
 B = zero(A);
 C = zero(A);