Create a more unified syntax for requesting a specific return type of functions?

round(1.1) --> 1.0
round(Int, 1.1) --> 1

works, but not

sqrt(Complex, -1)
acos(Complex, 1.1)
...

and other inverse trigonometric functions where the real domain is limited. In these cases the desired result is obtainable with

sqrt(-1 + 0m)  --> 0.0 + 1.0im
acos(1.1 + 0im) -->  0.0 - 0.4435682543851153im

which is concise, but feels a bit “hackish” or FORTRAN-like, or

sqrt(Complex(-1))  --> 0.0 + 1.0im
acos(Complex(1.1)) -->  0.0 - 0.4435682543851153im

where the intention is perhaps clearer but “object-oriented” syntax rather than "functional’.

Another candidate is

s=[2^x for x in 58:62]  # <--- some large integers
sum(x -> x*x, s)      # --> 0, incorrect result, instead
sum(BigInt, x -> x*x, s)  # would give the right result, but does not work
sum(x -> BigInt(x)*x, s) # is correct and equivalent, but less "nice"

So generally an optional first type argument would indicate that the return should be of this type, but also that the function will try to “do the right thing” to get the correct result. What exactly this means depends on the function.

Is this good or bad? Are there more candidate functions which usefully could be updated to support an optional first type argument?

2 Likes

I thought that maybe init= would work, but the only a bit nicer way I could find, is to convert s as in sum(x -> x*x, BigInt.(s))

As to whether I find that useful, I am not 100% sure, I think providing the inputs of right type to me sounds a bit more reasonable, like converting s instead?

1 Like

For my own internal functions, I’ve taken a liking to making a parametric callable type where the parameter is the return type.

For example

struct mysum{R} end

function (::mysum{R})(f, args...)::R where {R}
    x = zero(R)
    for vals in zip(args...)
        x+= convert(R, f(vals...))
    end
    x
end

mysum{BigInt}()(x -> x*x, s)
3 Likes

Providing the input as Complex and not Real is the present solution for functions like sqrt, acos, … In principle acos(Complex, 1.1) could be optimized a bit for Real input (but Complex output), which presently internally requires “ugly” tests whether the imaginary input is zero. For sum I cannot think of such a benefit that `sum(T, f, itr; … ) would have. Just it would perhaps be more consistent if many/most functions have this optional 1st type argument.

This would prevent us from using do block syntax with these higher order functions.

No idea whether this was the original intention, but from a math perspective, the current situation with sqrt etc. seems most intuitive to me. And I’m not sure why it is more or less object oriented than the alternatives.

round(Int, 1.1) seems to do the opposite (in some sense) to the proposed sqrt(Complex, -1) no? Rounding with the optional first argument restricts the output to a “narrower” type, whereas the square root would lift the input to a “wider” type before processing it. So I’m not sure if it’s good to mix these two meanings (although either one might have been chosen as convenience functions I guess).

With

Base.round(func::Function, args...; kwargs...) = func(round(args...; kwargs...))

do block syntax is possible?:

round(1.1) do x println(x) end        # --> 1.0
round(Int, 1.1) do x println(x) end   # --> 1

With sqrt(Complex(x)) the compiler is told how to compute the square root of x, namely by first converting x to a complex number. I think this is called “imperative” programming.

With sqrt(Complex, x) the compiler were told that Complex output is desired, but not how this should be achieved. (In this case the sqrt function could choose between a simplified algorithm if x is Real, but possibly negative, or a fully complex method for complex x, or the input could be converted to Complex anyway as above.) I think that this is called “declarative” or functional (?) programming (I’m not a CS.)

To more consistently embrace the “declarative” style, my proposal is that methods like sqrt(T, x), acos(T, x), sum(T, f, itr; …) could be made available.

Not all callable objects are a subtype of Function. One simple example is Complex which is a subtype of Number, yet even in your example, we can call Complex(x).

In my opinion, this is a weird comparison to make. For the round example, the Int argument is indeed controlling the output type you’re requesting, but for sqrt, it’s not actually about the output type but the input domain.

For example,

julia> sqrt(Complex{Int}(1))
1.0 + 0.0im

is not returning a Complex{Int}, it’s returning a Complex{Float64}. The reason sqrt(-1) throws an error is not that the return type isn’t what we want, it’s that the input is not in the domain of the function.

So I see no reason why this comparison should bother us or seem inconsistent, and I certainly don’t think it has anything to do with “object oriented” versus “functional”.

3 Likes

sum((x -> x*x) ∘ BigInt, s) or sum(abs2 ∘ BigInt, s) does this. For any function accepting a predicate, or a function like sqrt where Complex can be applied (for free or nearly so) to the input to achieve the desired effect, there just doesn’t seem to be much reason for this.

The reason round(Int, x) exists is because round(Int(x)) will not work (for non-integers). I would say that round(Int, x) is more akin to convert(Int, x), the only difference being that it allows inexact conversions by specifying the approximation method.

4 Likes

With

Base.sqrt(::Type{Complex{T}}, x::Real) where {T} = x>=zero(T) ? Complex{T}(√(x)) : Complex{T}(zero(T), √(-x))
Base.sqrt(::Type{Complex{T}}, x::Complex) where {T} = Complex{T}(√(x))
Base.sqrt(::Type{T}, x) where {T<:Real} = T(√(x))

I get:

sqrt(Float32, 2)  # --> 1.4142135f0
sqrt(Int, 1)  # --> 1
sqrt(Int, 4)  # --> 2
sqrt(Int, 2)  # --> ERROR: InexactError: Int64(1.4142135623730951)
sqrt(Complex{Float64}, 1)  # --> 1.0 + 0.0im
sqrt(Complex{Float64}, -1) # --> 0.0 + 1.0im
sqrt(Complex{Int}, 1)  # --> 1 + 0im
sqrt(Complex{Int}, -1)  # --> 0 + 1im
sqrt(Complex{Int}, 1+im)  # --> ERROR: InexactError: Int64(1.09868411346781)

The methods seem to extend the sqrt function consistently, amazing how easy it is with Julia!

BTW, is there any input z=Complex{Int}(x, y) with both x,y≠0 where sqrt(Complex{Int}, z) as defined above does not throw an InexactError (number theoretical problem)?

Personally not a fan because it’s a lot of manual forwarding in new barebones methods with no consistency on what happens.

round is a special case because you very often want an integer result but you may also want to retain the input type for further work, so the output type does become a worthy input. But look at the forwarding work it does besides the actual rounding:

round(::Type{T}, x) where T = round(T, x, RoundNearest)

round(::Type{T}, x, r::RoundingMode) where T = _round_convert(T, round(x, r), x, r)
_round_convert(::Type{T}, x_integer, x, r) where T = convert(T, x_integer)

It complicates the method table to essentially mirror convert, and for the general case, it’s much simpler to do it ourselves when functions expect particular input types or when we want to process an output type, with no loss in clarity.

If all the BigInt input would do is convert the output’s type like round does, then it would not change the incorrect result’s value BigInt(sum(x -> x*x, s)). You can already specify a BigInt input, and it gives the wrong results because it only affects the accumulated value, not the overflow while processing the input iterable:

julia> sum(x -> x*x, s; init=BigInt(0))
0

Changing the input iterable BigInt.(s) or the elementwise function x -> BigInt(x)*x is necessary for evading overflow, and a BigInt input that does either automatically adds ambiguity about what happens. A user may even assume it’s only converting the output value’s type like convert or round does, with frustrating outcomes when they start tweaking code. Considering that the clearer code is the same amount of work, it’s worth doing in my opinion.

1 Like

I would argue that it doesn’t depend on the desired output whether the square root of x can be computed, but only on what x itself actually is.

It just seems to be the most practical solution to distinguish between reals and complex numbers, since many (most?) use cases use strictly real arithmetic and don’t (want to) assume that every number is actually a complex number and can be treated as such.

But, granted, where to draw the line does seem to be a bit arbitrary. E.g. sqrt(2) doesn’t throw an error, but returns the square root of the real number 2 … so integers are sometimes (always?) implicitly assumed to be real numbers (which is also reflected in their type relation Int <: Real), although in some situations it might be necessary to stay within the integers. Again, it might simply come down to what is most useful in most of the use cases, which is presumably real arithmetic.

I don’t think this is more consistent. Something like sqrt(Int, 4) which would return 2 (as Int) would have the same meaning as round(Int, 1.1) to my understanding (“do the operation and convert whatever comes out to the desired type if possible”). Mixing this with “interpret the input in whatever way that gives me the desired output” seems like a very different thing to do, and is sometimes already done, e.g. in sqrt(2).

Marked as solution :slight_smile: . I agree that sum([T, ] ...) is not a good example for my point, it is unclear what this would mean.

However, for sqrt([T, } ... and other functions which are normally used in the real domain, but can also work in the complex one, I still think that it is a good idea. But this could of course be implemented from case to case, perhaps also in a package.

Yeah, for various reasons, we ended up with the situation where we consider the integers and rationals to be embedded in the reals, but we do not consider the reals to be embedded in the complex numbers.

The performance gain for a sqrt(::Type{T}, ...) compared to the standard sqrt is larger than expected. My benchmark assumes that I need to compute the square root of a large number of Floats, which are both positive and negative. Therefore the output needs to be Complex anyway. The standard sqrt forces to convert also the input to Complex and do the sqrt computation with Complex machinery, which is unnecessary in this case and can be avoided with

Base.sqrt(::Type{Complex{T}}, x::Real) where {T} = x>=zero(T) ? Complex{T}(√(x)) : Complex{T}(zero(T), √(-x))

Results:

ta=randn(100_000);
base_sqrt(ta) = begin                                                                                                                                                                                                                 
▌     foreach(ta) do x                                                                                                                                                                                                                  
▌         s = sqrt(Complex(x))                                                                                                                                                                                                          
▌     end                                                                                                                                                                                                                               
▌ end
my_sqrt(ta) = begin                                                                                                                                                                                                                   
▌     foreach(ta) do x                                                                                                                                                                                                                  
▌         s = sqrt(Complex{Float64}, x)                                                                                                                                                                                                 
▌     end                                                                                                                                                                                                                               
▌ end
@btime base_sqrt($ta)  # --> 1.320 ms (0 allocations: 0 bytes)
@btime my_sqrt($ta)   # --> 403.834 μs (0 allocations: 0 bytes)

The “shortcut” for complex output, but real input is ~three times faster!

With

Base.sqrt(::Type{Complex{T}}, x::Complex) where {T} = Complex{T}(√(x))

some fun stuff in the Complex plane of integers is revealed, to my surprise:

sqrt(Complex{Int}, 1+im) gives InexactError because there is no integer square root of 1+im. Perhaps not completely unexpected, a solution is given for Pythagorean triplets:

sqrt(Complex{Int}, 3+4im)  # --> 2 + 1im

but sqrt(Complex{Int}, 4+3im) does not work. Sometimes the larger of the integers needs to be on the real, not imagninary axis:

sqrt(Complex{Int}, 21+20im)  # --> 5 + 2im

and sqrt(Complex{Int}, 20+21im) does not work. …

I’m sure that there are mathematical explanations for these observations and all this is investigated and published somewhere. For me it is suprising and fun that a one-line method definition of sqrt in Julia enables to play with this.

1 Like

Indeed, there are situations where this can be useful. But no other method in Base (and probably most of the ecosystem) supports a ::Type first argument (and there’s been little demand for such), so adding a Base.sqrt method buys you virtually nothing (and is type piracy except as a direct change to Base). It seems you’d do just as well to define your own local function

csqrt(x::Real) = x < 0 ? Complex(zero(x), sqrt(-x)) : Complex(sqrt(x), zero(x))
csqrt(x::Number) = sqrt(x)

This should give equal performance to your definition with no need for piracy. Personally, I find round(Complex{Int}, csqrt(x), RoundNearest #= or other RoundingMode =#) to be more clear and useful than sqrt(Complex{Int}, x) could ever be. It’s a rare situation that one wants a Complex{Int} from sqrt, since it’s so often inexact, but when one does I suspect one would usually be content with rounding (in which case one wants to be able to control the RoundingMode).

As for nonnegative integer sqrt, there is the isqrt which is basically round(Int, sqrt(x), RoundDown). I don’t think there exists an agreeable definition of isqrt(::Complex) (except when the result is exact), so one should probably use a round of sqrt for that.