Tests pass

Author: dlfivefifty

src/ApproxFunFourier.jl

@@ -11,32 +11,33 @@ import FFTW: plan_r2r!, fftwNumber, REDFT10, REDFT01, REDFT00, RODFT00, R2HC, HC
                 r2r!, r2r,  plan_fft, plan_ifft, plan_ifft!, plan_fft!
 
 import ApproxFunBase: normalize!, flipsign, FiniteRange, MatrixFun, UnsetSpace, VFun, RowVector,
-                    UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, 
+                    UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, NoSpace, Space,
                     IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
                     AnyDomain, ZeroSpace, TrivialInterlacer, BlockInterlacer, 
                     AbstractTransformPlan, TransformPlan, ITransformPlan,
                     ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral,
-                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper,
+                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation,
                     Conversion, Multiplication, Derivative, Integral, bandwidths, 
                     ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                     DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral,
-                    ReverseOrientation, Reverse, Dirichlet,
+                    ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven, Dirichlet,
                     TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                     hasfasttransform, canonicalspace, setdomain, prectype, domainscompatible, 
-                    plan_transform, plan_itransform, transform, itransform, hasfasttransform, Integral, 
+                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform, Integral, 
                     domainspace, rangespace, boundary, 
-                    union_rule, conversion_rule, maxspace_rule, conversion_type, hasconversion, points, 
+                    union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points, 
                     rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp, 
                     linesum, differentiate, integrate, linebilinearform, bilinearform, 
                     UnsetNumber, coefficienttimes,
-                    Segment, isambiguous, Vec, eps, 
+                    Segment, isambiguous, Vec, eps, isperiodic,
                     arclength, complexlength,
                     invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
-                    reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients,
+                    reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients, isconvertible,
                     clenshaw, ClenshawPlan, sineshaw,
                     toeplitz_getindex, toeplitz_axpy!, ToeplitzOperator, hankel_getindex, 
                     SpaceOperator, ZeroOperator, InterlaceOperator,
-                    interlace!, reverseeven!, negateeven!, cfstype
+                    interlace!, reverseeven!, negateeven!, cfstype, pad!,
+                    extremal_args, hesseneigvals
 
 import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
             Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,
@@ -66,6 +67,7 @@ import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_che
 
 export Fourier, Taylor, Hardy, CosSpace, SinSpace, Laurent, PeriodicDomain
 
+include("utils.jl")
 include("Domains/Domains.jl")
 
 for T in (:CosSpace,:SinSpace)
@@ -602,5 +604,25 @@ include("specialfunctions.jl")
 include("FourierOperators.jl")
 include("LaurentOperators.jl")
 include("LaurentDirichlet.jl")
+include("roots.jl")
+
+Fun(::typeof(identity), d::Circle) = Fun(Laurent(d),[d.center,0.,d.radius])
+
+Space(d::PeriodicDomain) = Fourier(d)
+Space(d::Circle) = Laurent(d)
+
+
+## Evaluation
+
+Evaluation(d::PeriodicDomain,x::Number,n...) = Evaluation(Laurent(d),complex(x),n...)
+
+## Definite Integral
+
+DefiniteIntegral(d::PeriodicDomain) = DefiniteIntegral(Laurent(d))
+DefiniteLineIntegral(d::PeriodicDomain) = DefiniteLineIntegral(Laurent(d))
+
+## Toeplitz
+union_rule(A::Space{<:PeriodicSegment}, B::Space{<:IntervalOrSegment}) =
+    union(Space(Interval(domain(A))), B)
 
 end #module

src/Domains/Domains.jl

@@ -3,6 +3,7 @@
 
 abstract type PeriodicDomain{T} <: Domain{T} end
 
+isperiodic(::PeriodicDomain) = true
 
 canonicaldomain(::PeriodicDomain) = PeriodicSegment()
 

src/roots.jl

@@ -0,0 +1,73 @@
+
+function companion_matrix(c::Vector{T}) where T
+    n=length(c)-1
+
+    if n==0
+        zeros(T,0,0)
+    else
+        A=zeros(T,n,n)
+        for k=1:n
+            A[k,end]=-c[k]/c[end]
+        end
+        for k=2:n
+            A[k,k-1]=one(T)
+        end
+        A
+    end
+end
+
+
+# if isdir(Pkg.dir("AMVW"))
+#     using AMVW
+#     function complexroots(cfs::Vector)
+#         c=chop(cfs,10eps())
+#
+#         # Only use special routine for large roots
+#         if length(c)≥70
+#             Main.AMVW.rootsAMVW(c)
+#         else
+#             hesseneigvals(companion_matrix(c))
+#         end
+#     end
+# else
+complexroots(cfs::Vector{T}) where {T<:Union{Float64,ComplexF64}} =
+    hesseneigvals(companion_matrix(chop(cfs,10eps())))
+# end
+
+function complexroots(cfs::Vector{T}) where T<:Union{BigFloat,Complex{BigFloat}}
+    a = Fun(Taylor(Circle(BigFloat)),cfs)
+    ap = a'
+    rts = Array{Complex{BigFloat}}(complexroots(Vector{ComplexF64}(cfs)))
+    # Do 3 Newton steps
+    for _ = 1:3
+        rts .-= a.(rts)./ap.(rts)
+    end
+    rts
+end
+
+complexroots(neg::Vector, pos::Vector) =
+    complexroots([reverse(chop(neg,10eps()), dims=1);pos])
+complexroots(f::Fun{Laurent{DD,RR}}) where {DD,RR} =
+    mappoint.(Ref(Circle()), Ref(domain(f)),
+        complexroots(f.coefficients[2:2:end],f.coefficients[1:2:end]))
+complexroots(f::Fun{Taylor{DD,RR}}) where {DD,RR} =
+    mappoint.(Ref(Circle()), Ref(domain(f)), complexroots(f.coefficients))
+
+
+
+function roots(f::Fun{Laurent{DD,RR}}) where {DD,RR}
+    irts=filter!(z->in(z,Circle()),complexroots(Fun(Laurent(Circle()),f.coefficients)))
+    if length(irts)==0
+        Complex{Float64}[]
+    else
+        rts=fromcanonical.(f, tocanonical.(Ref(Circle()), irts))
+        if isa(domain(f),PeriodicSegment)
+            sort!(real(rts))  # Make type safe?
+        else
+            rts
+        end
+    end
+end
+
+
+roots(f::Fun{Fourier{D,R}}) where {D,R} = roots(Fun(f,Laurent))

src/specialfunctions.jl

@@ -42,3 +42,21 @@ function Base.imag(f::Fun{Laurent{DD,RR}}) where {DD,RR}
 
     Fun(Fourier(domain(f)),ret)
 end
+
+
+function log(f::Fun{Fourier{D,R},T}) where {T<:Real,D,R}
+    if isreal(domain(f))
+        cumsum(differentiate(f)/f)+log(first(f))
+    else
+        # this makes sure differentiate doesn't
+        # make the function complex
+        g=log(setdomain(f,PeriodicSegment()))
+        setdomain(g,domain(f))
+    end
+end
+
+extremal_args(f::Fun{<:Space{<:PeriodicSegment}}) = roots(differentiate(f))
+function extremal_args(f::Fun{<:Space{<:PeriodicDomain}})
+    S=typeof(space(f))
+    fromcanonical.(f,extremal_args(setcanonicaldomain(f)))
+end

src/utils.jl

@@ -0,0 +1,84 @@
+# diff from CosSpace -> SinSpace
+
+function cosspacediff(v::AbstractVector{T}) where T<:Number
+    if length(v)==1
+        w = zeros(T,1)
+    else
+        w = Array{T}(undef, length(v)-1)
+        for k=1:length(v)-1
+            @inbounds w[k] = -k*v[k+1]
+        end
+    end
+
+    w
+end
+
+# diff from SinSpace -> CosSpace
+
+function sinspacediff(v::AbstractVector{T}) where T<:Number
+    w = Array{T}(undef, length(v)+1)
+    w[1] = zero(T)
+    for k=1:length(v)
+        @inbounds w[k+1] = k*v[k]
+    end
+
+    w
+end
+
+# diff from Fourier -> Fourier
+
+function fourierdiff(v::AbstractVector{T}) where T<:Number
+    n = 2(length(v)÷2)+1
+    w = Array{T}(undef, n)
+    w[1] = zero(T)
+    n > 1 && (w[n-1] = zero(T))
+    for k=1:n÷2-1
+        @inbounds w[2k] = -k*v[2k+1]
+        @inbounds w[2k+1] = k*v[2k]
+    end
+    n > 1 && (w[n] = (n÷2)*v[n-1]; n == length(v) && (w[n-1] = -(n÷2)*v[n]))
+
+    w
+end
+
+# diff from Taylor -> Taylor
+
+function taylor_diff(v::AbstractVector{T}) where T<:Number
+    w = Array{T}(undef, length(v))
+    for k=1:length(v)
+        @inbounds w[k] = (k-1)*v[k]
+    end
+
+    w
+end
+
+# diff from Hardy{false} -> Hardy{false}
+
+function hardyfalse_diff(v::AbstractVector{T}) where T<:Number
+    w = Array{T}(undef, length(v))
+    for k=1:length(v)
+        @inbounds w[k] = -k*v[k]
+    end
+
+    w
+end
+
+# diff from Laurent -> Laurent
+
+function laurentdiff(v::AbstractVector{T}) where T<:Number
+    n = length(v)
+    w = Array{T}(undef, n)
+    w[1] = zero(T)
+    n=length(v)
+
+    for k=1:(isodd(n) ? n÷2 : n÷2-1)
+        @inbounds w[2k] = -k*v[2k]
+        @inbounds w[2k+1] = k*v[2k+1]
+    end
+
+    if iseven(n)
+        @inbounds w[n] = -(n÷2)*v[n]
+    end
+
+    w
+end

test/runtests.jl

@@ -1,6 +1,6 @@
 using ApproxFunFourier, ApproxFunBase, Test, BlockArrays, BlockBandedMatrices, SpecialFunctions, LinearAlgebra
     import ApproxFunBase: testspace, testtransforms, testmultiplication,
-                      testbandedoperator, testcalculus, Block
+                      testbandedoperator, testcalculus, Block, Vec, testfunctional
     import SpecialFunctions: factorial
 
 @testset "Periodic Domains" begin    
@@ -114,8 +114,6 @@ end
     v = randn(4) .+ im.*randn(4)
     @test P*v == P*real(v) + im*(P*imag(v))
 
-    @test norm(Fun(x->Fun(cos,Fourier(-π .. π),20)(x),20)-Fun(cos,20)) <100eps()
-    @test norm(Fun(x->Fun(cos,Fourier(-π .. π))(x))-Fun(cos)) <100eps()
     @test norm(Fun(cos,Fourier)'+Fun(sin,Fourier)) < 100eps()
 
     f=Fun(x->exp(-10sin((x-.1)/2)^2),Fourier)
@@ -156,8 +154,8 @@ end
     A = Multiplication(mySin,Fourier())
     @test A.op[1,1] == 0
 
-    @test norm(ApproxFun.Reverse(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(-t-0.2)-0.1),Fourier())) < 10eps()
-    @test norm(ApproxFun.ReverseOrientation(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(t-0.2)-0.1),Fourier(PeriodicSegment(2π,0)))) < 10eps()
+    @test norm(ApproxFunBase.Reverse(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(-t-0.2)-0.1),Fourier())) < 10eps()
+    @test norm(ApproxFunBase.ReverseOrientation(Fourier())*Fun(t->cos(cos(t-0.2)-0.1),Fourier()) - Fun(t->cos(cos(t-0.2)-0.1),Fourier(PeriodicSegment(2π,0)))) < 10eps()
 end
 
 @testset "Laurent" begin
@@ -169,7 +167,6 @@ end
     @test Fun(Laurent(-1..1),[1,1.,1.])(0.1) ≈ 1+2cos(π*(0.1+1))
     @test Fun(Laurent(0..1),[1,1.,1.])(0.1) ≈ 1+2cos(2π*0.1)
 
-    @test norm(Fun(x->Fun(cos,Laurent)(x))-Fun(cos)) <100eps()
     @test norm(Fun(cos,Laurent)'+Fun(sin,Laurent)) < 100eps()
 
     B=Evaluation(Laurent(0..2π),0,1)
@@ -367,20 +364,11 @@ end
 end
 
 @testset "Integral equations" begin    
-    @time for S in (Fourier(Circle()),Laurent(Circle()),Taylor(Circle()),
-        CosSpace(Circle()),JacobiWeight(-0.5,-0.5,Chebyshev()),
-        JacobiWeight(-0.5,-0.5,Chebyshev(Segment(1.0,2.0+im))),
-        JacobiWeight(0.5,0.5,Ultraspherical(1,Segment(1.0,2.0+im))))
+    @time for S in (Fourier(Circle()),Laurent(Circle()),Taylor(Circle()),CosSpace(Circle()))
         testfunctional(DefiniteLineIntegral(S))
     end
 
     Σ = DefiniteIntegral()
-
-    f1 = Fun(t->cos(cos(t)),-π..π)
-    f = Fun(t->cos(cos(t)),Laurent(-π..π))
-
-    @test sum(f1) ≈ Σ*f
-
     f1 = Fun(t->cos(cos(t))/t,Laurent(Circle()))
     f2 = Fun(t->cos(cos(t))/t,Fourier(Circle()))
     @test Σ*f1 ≈ Σ*f2