Problem with `roots` and `BigFloat`

[lcw]

Hi,

I am having an issue with roots and BigFloat. Please see below for minimum example reproducing the error. I am happy to debug further but I am not sure what to do next.

Thanks for looking at it.

Best, Lucas

❯ julia
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.2 (2017-12-13 18:08 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-apple-darwin14.5.0

julia> Pkg.status("ApproxFun")
 - ApproxFun                     0.8.0-             master

julia> using ApproxFun

julia> x = Fun(identity,BigFloat(0)..BigFloat(10))
Fun(Chebyshev(【0.000000000000000000000000000000000000000000000000000000000000000000000000000000,1.000000000000000000000000000000000000000000000000000000000000000000000000000000e+01】),BigFloat[5.000000000000000000000000000000000000000000000000000000000000000000000000000000, 5.000000000000000000000000000000000000000000000000000000000000000000000000000069])

julia> roots(abs(x^2))
ERROR: MethodError: no method matching plan_ichebyshevtransform(::SubArray{BigFloat,1,Array{BigFloat,1},Tuple{StepRange{Int64,Int64}},true})
Closest candidates are:
  plan_ichebyshevtransform(::Array{T<:Union{BigFloat, Complex{BigFloat}},1}; kind) where T<:Union{BigFloat, Complex{BigFloat}} at /Users/lucas/.julia/v0.6/ApproxFun/src/Extras/fftGeneric.jl:51
  plan_ichebyshevtransform(::Array{D<:DualNumbers.Dual,1}; kind) where D<:DualNumbers.Dual at /Users/lucas/.julia/v0.6/ApproxFun/src/Extras/dualnumbers.jl:46
  plan_ichebyshevtransform(::AbstractArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},1}; kind) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/lucas/.julia/v0.6/ApproxFun/src/LinearAlgebra/chebyshevtransform.jl:93
Stacktrace:
 [1] roots(::ApproxFun.Fun{ApproxFun.Chebyshev{ApproxFun.Segment{BigFloat},BigFloat},BigFloat,SubArray{BigFloat,1,Array{BigFloat,1},Tuple{StepRange{Int64,Int64}},true}}) at /Users/lucas/.julia/v0.6/ApproxFun/src/Extras/roots.jl:75
 [2] _mapreduce(::ApproxFun.#roots, ::Base.#vcat, ::IndexLinear, ::Array{ApproxFun.Fun,1}) at ./reduce.jl:273
 [3] roots(::ApproxFun.Fun{ApproxFun.PiecewiseSpace{Tuple{ApproxFun.Chebyshev{ApproxFun.Segment{BigFloat},BigFloat},ApproxFun.Chebyshev{ApproxFun.Segment{BigFloat},BigFloat},ApproxFun.Chebyshev{ApproxFun.Segment{BigFloat},BigFloat}},ApproxFun.UnionDomain{Tuple{ApproxFun.Segment{BigFloat},ApproxFun.Segment{BigFloat},ApproxFun.Segment{BigFloat}},BigFloat},BigFloat},BigFloat,Array{BigFloat,1}}) at /Users/lucas/.julia/v0.6/ApproxFun/src/Extras/roots.jl:416

[MikaelSlevinsky]

This reveals a couple issues. Let x = Fun(identity,BigFloat(0)..BigFloat(10)).

• f = abs(x^2) returns a domain with three pieces, one piece is actually a point,
• values(f) fails because the transform isn't expecting subarrays with step ranges,
• roots(abs2(x)) returns an unacceptable answer.

Note that if x = Fun(identity, 0..10), then roots(abs(x^2)) ≈ [7.45058e-8] and roots(abs2(x)) ≈ [0.0, 7.45058e-8] which is reasonable.

[dlfivefifty]

I fixed the transform on master. But now it returns BigFloat(NaN)

[dlfivefifty]

This works however: roots(abs(x^BigFloat(2)))

This is because x^2.0 returns a Fun with JacobiWeight space.

[dlfivefifty]

The reason roots(abs2(x)) fails is because roots for BigFloat calls standard roots and then does a small number of Newton iterations:

https://github.com/JuliaApproximation/ApproxFun.jl/blob/72d88dfe546687a0c1e7e4aaf039bf9997bf11e0/src/Extras/roots.jl#L89

I think this is misguided, as one reason you use BigFloat is to deal with cases where Float64 is not accurate enough, in which case the Float64 roots might be missing some of the BigFloat roots.

The solution is to rewrite the rest of the roots commands to work with BigFloat. But this takes some work:

• [ ] Replace colleague_eigvals with a routine that supports general types, maybe AMVW.jl?
• [ ] Change splitPoint to be type dependent.

[MikaelSlevinsky]

I thought this was a while loop. Someone must've changed it to only do 3 iterations.

There are always going to be corner cases with roots, especially multiple roots, but somehow about 8 digits (sqrt(eps())) seems to me like it should be sufficient to starting Newton in BigFloat.

So far, AMWV.jl is only for monomial bases?

[dlfivefifty]

Newton iteration defeats the point of allow BigFloat, which is to ramp up the accuracy to get better results.

Take for example:

1E-15sin(1000x) + sin(x)

If you use Newton iteration you will only see 1 root, when there are hundreds of roots.

It's also possible for Newton iteration to diverge.

[MikaelSlevinsky]

Can you show me some of the nontrivial roots? I'm looking on google's plotter and only see the obvious ones. (In fact, for a prefactor smaller than something like 4.5e-3, say 4.5e-3sin(1000x) + sin(x), there might not be any others beyond piZ)

[dlfivefifty]

Sorry, meant Sin(100000000000 x) 😛