sincos is broken for nonstandard types

[dpsanders]

Due to

julia/base/math.jl

Line 428 in 3e3f966

res = Base.FastMath.sincos_fast(x)

:

 res = Base.FastMath.sincos_fast(x)

This breaks sin, cos, exp etc. for Complex{T} with non-standard T (i.e. not Float64, Float32), that previously worked fine (but apparently weren't tested).

Presumably it's sufficient to restrict the argument type of the above sincos definition and add a generic fallback? Should that be

sincos(x) = (sin(x), cos(x))

or is there a better solution?

[yuyichao]

What's an example of non-standard T and what's the error message?

The fallback method is defined for sincos_fast and I don't see why that line will cause any error.

[dpsanders]

julia> using IntervalArithmetic

julia> a = (1..1) + (2..2)*im
[1, 1] + [2, 2]*im

julia> typeof(a)
Complex{IntervalArithmetic.Interval{Float64}}

julia> sin(a)
ERROR: TypeError: sin: in typeassert, expected AbstractFloat, got IntervalArithmetic.Interval{Float64}
Stacktrace:
 [1] sincos at ./math.jl:428 [inlined]
 [2] sin(::Complex{IntervalArithmetic.Interval{Float64}}) at ./complex.jl:701

[yuyichao]

k. so it's throwing an error due to

julia/base/fastmath.jl

Line 313 in 3e3f966

sincos_fast(v::Real) = sincos_fast(float(v)::AbstractFloat)

since float doesn't return an AbstractFloat. The assertion is to avoid infinite recursion so it shouldn't be removed. The alternative way I've thought about but decide not to do before seeing a float that triggers the assertion is

function sincos_fast(v::Real)
    vf = float(v)
    if isa(vf, AbstractFloat)
        sincos_fast(vf)
    else
        sin_fast(vf), cos_fast(vf)
    end
end

Another way is to restrict the input type to Integer

sincos_fast(v::Integer) = sincos_fast(float(v))

And rely on float(::Integer) to not return another integer type...

[dpsanders]

Ah, I see, thanks. And you can't just have two different fallbacks, one for v::AbstractFloat and the other for v::Real?

[yuyichao]

There are two fallbacks

julia/base/fastmath.jl

Lines 312 to 314 in 3e3f966

	sincos_fast(v::AbstractFloat) = (sin_fast(v), cos_fast(v))
	sincos_fast(v::Real) = sincos_fast(float(v)::AbstractFloat)
	sincos_fast(v) = (sin_fast(v), cos_fast(v))

and the isa(vf, AbstractFloat) method is basically implementing an fallback for Real.

[giordano]

There are other functions that don't work with Complex{IntervalArithmetic.Interval{Float64}} and aren't related with sincos, like log1p, exp2, exp10.

[dpsanders]

Also log:

julia> using IntervalArithmetic

julia> z = (1.1..1.2) + (2.1..2.2)im
[1.09999, 1.20001] + [2.09999, 2.20001]*im

julia> log(z)
ERROR: StackOverflowError:
Stacktrace:
 [1] log(::Complex{IntervalArithmetic.Interval{Float64}}) at ./complex.jl:470 (repeats 80000 times)

[giordano]

...and also sqrt, log10, log2, asin, acos, atan, and all hyperbolic functions. And sinpi, cospi, sinc, cosc, depending on the branch they enter:

julia> sinpi((1..2) + (3..4)im)
ERROR: StackOverflowError:
Stacktrace:
 [1] sinpi(::IntervalArithmetic.Interval{Float64}) at ./special/trig.jl:220 (repeats 80000 times)

julia> cospi((1..2) + (3..4)im)
ERROR: StackOverflowError:
Stacktrace:
 [1] cospi(::IntervalArithmetic.Interval{Float64}) at ./special/trig.jl:221 (repeats 80000 times)

julia> sinc((1..2) + (3..4)im)
ERROR: StackOverflowError:
Stacktrace:
 [1] sinc(::Complex{IntervalArithmetic.Interval{Float64}}) at ./special/trig.jl:298 (repeats 80000 times)

julia> cosc((1..2) + (3..4)im)
ERROR: StackOverflowError:
Stacktrace:
 [1] cosc(::Complex{IntervalArithmetic.Interval{Float64}}) at ./special/trig.jl:310 (repeats 80000 times)

[dpsanders]

Thanks for the list.

[dpsanders]

@yuyichao What is special about AbstractFloat here?

[yuyichao]

where?

[dpsanders]

In sincos_fast.

[stevengj]

There are lots of cases where float(x) does not produce an AbstractFloat, e.g. float(1+0im), so that assertion doesn't seem quite right to me. It is perfectly reasonable for float on an Interval type to produce another Interval type, for example.

One option would be to simply let it throw a StackOverflow error if float(v) === v here. That seems simplest to me; it's not like you really get "infinite recursion", right?

Another option would be to use something like sincos_fast(v::Real) = sincos_fast(forcefloat(v)) where

function forcefloat(v)
    f = float(v)
    typeof(f) === typeof(v) && error("float(::",typeof(v),") produces a result of the same type")
    return f
end

(The compiler optimizes out the type check in cases where inference succeeds.)

[yuyichao]

In sincos_fast.

#22095 (comment)

e.g. float(1+0im),

This case is irrelevant since it's not ::Real

One option would be to simply let it throw a StackOverflow error if float(v) === v here. That seems simplest to me; it's not like you really get "infinite recursion", right?

I don't think that fixes any problem? You will get a infinite recursion.

Another option would be to use something like sincos_fast(v::Real) = sincos_fast(forcefloat(v)) where

It's basically what the current code is doing. forcefloat implemented this way won't actually prevent a infinite recursion.

[pkofod]

I've thought about but decide not to do

I figure there's an obvious reason, but why?

[yuyichao]

I did not know any exception so I have no way to back it up in case anyone asks (and I'm 100% sure someone would). Especially since the infinite recursion in other functions clearly shows that it is not currently done in any other part in Base.

[giordano]

Especially since the infinite recursion in other functions clearly shows that it is not currently done in any other part in Base.

Yes, I believe there are several places in Base where it's implicitly assumed that float(x::Real) returns an AbstractFloat. I wouldn't say that's a too bad assumption.

[dpsanders]

My question was more along the lines of "why is AbstractFloat singled out?". I think the answer is simply because there are specialised versions of sincos defined in Base for each of the subtypes of AbstractFloat.

It is definitely a bug that all these complex functions that have relatively easy fallback definitions in terms of their real counterparts (ignoring subtleties of branch cuts etc) do not work for such "non-standard" complex types.

[stevengj]

@yuyichao, "infinite recursion" is not infinite because it eventually throws a StackOverflowError. So, it's basically a choice between throwing one exception or another, no?

The point of my forcefloat approach is that it catches recursion loops, but still works for types in which float(::T) does not produce an abstractfloat as long as there is an appropriate method.

[stevengj]

How about

function sincos(x::Real)
   f = float(x)
   if typeof(f) === typeof(x)
      sincos(f)
   else
      (sin(f), cos(f))
   end
end

[yuyichao]

I think the answer is simply because there are specialised versions of sincos defined in Base for each of the subtypes of AbstractFloat.

No. It's because the fallback is defined for AbstractFloat.

"infinite recursion" is not infinite because it eventually throws a StackOverflowError. So, it's basically a choice between throwing one exception or another, no?

..... Errr, well, StackOverflowError IS infinite recursion, in that real infinite recursion (without overflowing the stack) will not happen currently. If you don't like how it is described I can use "useless/unhelpful error". Note that this is also a case where if we disable certain codegen llvm options can actually infinitely recurse without raising any error.

The point of my forcefloat approach is that it catches recursion loops, but still works for types in which float(::T) does not produce an abstractfloat as long as there is an appropriate method.

There could be cases where this is the case. Certainly not for IntervalArithmetic.

How about

I personally prefer not using typeof(before) === typeof(after) as the check since it doesn't theoritically avoid infinite recursion (reads: StackOverflowError). Given that the intent of the float call is to covert to AbstractFloat I think the original version I proposed above is clearer.