Added GCD methods for complex numbers #35374
Conversation
Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>
Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>
Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>
* faster BigInt hashing
When the input type is a dispatch tuple, we know that a subtype lookup
is sufficient, so we can do a slightly more direct lookup and avoid
considering `Union{}`-declared methods as possible results.
This also fixes an issue with the special `jl_matching_methods` query
being done by `jl_insert_backedges`: we would set the `done` flag when a
method "fully covered" the result, even if that method was deleted.
Rather than continue to abuse the meaning of `world == 0` (which also
has performance implications), we add a specific flag for this case.
…ang#34861) The function was moved by 65b8e7e
…uliaLang#34764)" (JuliaLang#34878) This reverts commit 4ab8a9d.
…mate (JuliaLang#34859) equality since BLAS and Julia might use square roots that round differently. This has occasionally caused CI failures.
This allows julia code to more easily consume the output of jl-parse-all (see JuliaLang#34715) and is consistent with the way line number nodes have file name information in all non-toplevel constructs.
This changes unnamed type variables to print as valid syntax, via the
`var""` printing macro.
For example, before this commit, `Vector{<:Bool}` printed as a
non-syntactic string:
```julia
julia> Vector{<:Bool}
Array{#s16,1} where #s16<:Bool
```
And now it prints as a valid, syntactic expression:
```julia
julia> Vector{<:Bool}
Array{var"#s46",1} where var"#s46"<:Bool
```
We have deleteat!(vec::Vector, itr), but it returns the modified vector. If you instead want to obtain the delted elements, there is splice!(::Vector, ::UnitRange) as a special case of the three-argument version that inserts additional elements where the old ones were deleted. However, these is no two argument splice! for arbitrary iterables, so probably the best way to write that currently is: ``` inds = collect(itr) vals = A[inds] deleteat!(A, inds) vals ``` which is both less efficient and more verbose than ideal. I'm wondering if perhaps deleteat! should have been returning the deleted elements, but that's not a change we can make in 1.x. Instead, I'm proposing here to extend the 2 argument (but not the 3 argument for obvious reasons) variant of splice! to arbitrary iterables.
* More tests for Triangular, div and mult * Long form function defs * avoid useless unwrapping-wrapping in triangular multiplication Co-authored-by: Daniel Karrasch <Daniel.Karrasch@gmx.de>
When we added the check to prevent inlining through ambiguous methods, we failed exclude this check for calls to `invoke` (which skip ambiguous methods, since the method is explicitly selected). Fixes JuliaLang#34900.
Here we implement a syntax for ccall with Julia-like type annotations on the arguments. Compared to ccall: * The new syntax is more familiar to Julia programmers * The new syntax gives a notation for calling C varargs functions * Support for specifying calling convention is not yet implemented as that will require another syntax discussion * The semantics for interpolating "function like things" is much simplified (only function pointers are allowed)
* Tranpose for non-numeric eltype bidiag tests * More tests for bidiag mul
This fixes the case in JuliaLang#34482 and adds some more robustness for similar cases, though I wouldn't be surprised if there were more dragons hiding here. Intrinsics are a bit special and from the start, we sort of expected them to only ever be called correctly under pain of segfaults or other undefined behavior. We've been gradually making these more robust, but fundamentally, they were never intended to be used by users directly, only through the type-validating wrappers in Base. What has changed to cause the recent influx of issues in this area is that people now like to do compiler transforms that happily recurse through these wrappers and perform transforms that are not always legal on Intrinsics. This should help catch a number of them, but this interface is still not very thoroughly validated, and I would be surprised to see crashes or other errors stemming from incorrect usage here.
Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>
|
Would be good to get someone with complex number expertise to review. |
Co-Authored-By: Stefan Karpinski <stefan@karpinski.org>
simeonschaub
left a comment
simeonschaub
left a comment
There was a problem hiding this comment.
Overall looks good, but see my comments about checking for overflow. I think the div and rem should also have their own docstring explaining how they work.
base/complex.jl
Outdated
| a, b = Complex{R}(a), Complex{R}(b) | ||
| b̅ = conj(b) | ||
| # TODO: Handle overflow when calculating a*b̅ | ||
| t = a*b̅ |
There was a problem hiding this comment.
Since div can't overflow for real numbers, I think, we should try to guarantee that div doesn't overflow for complex numbers either. I think this safety is worth the extra performance cost in most uses of this function. Implementing Base.checked_mul(::Complex{<:Integer}, ::Complex{<:Integer}) should actually be fairly straightforward to do.
There was a problem hiding this comment.
That's true, but an issue with the current division algorithm (a/b = a*b̅/abs2(b)), is that due to abs2 overflow, we cannot ensure that a/b never overflows. We can indeed try to catch and handle all the possible overflows, but that will be it.
We can always find some inputs which have b::Complex{T} and abs2(b)>typemax(T), which will make div(a, b) overflow.
E.g. a = Complex{Int32}(1, 1), b= Complex{Int32}(10^5, 10^5). Here, a/b will overflow due to the current algorithm.
There was a problem hiding this comment.
I'm not sure, what you mean. I mean, we check in abs2 and in *, that nothing overflows and then you're just calling div on real integers, which can't overflow, so we should catch all cases of overflow. I'm wondering though, if there aren't any better algorithms for division of complex integers, because checking in each step will prevent overflow in the result, but we will also throw errors for cases where either a*b' or abs2(b) overflows, but a/b should actually still be representable as the complex integer type of a and b.
There was a problem hiding this comment.
we will also throw errors for cases where either
a*b'orabs2(b)overflows
This is what I meant. It is possible that abs2(b) overflows, but a/b is still representable as a complex integer. E.g. div(Complex{Int32}(1, 1), Complex{Int32}(10^5, 10^5)).
Also, should checked_mul(::Complex, ::Complex) be in complex.jl or in checked.jl ?
There was a problem hiding this comment.
As soon as one of abs(b) or a*b' overflows, isn't the end result also wrong with your current algorithm, even if the result might still be representable?
There was a problem hiding this comment.
@ravibitsgoa Could you explain your usecase for this? That would probably help to understand some of the constraints
There was a problem hiding this comment.
#35125
This PR is based on this issue. After quite a long discussion there, we decided we should avoid using floating points or rationals or wider type internally. Although those three are quite decent solutions IMO, but they will create a problem when we want to find exact integer answers of div, rem and gcd.
There was a problem hiding this comment.
I really cannot imagine those overflow-checked stuff (if you get it right) will ever be more efficient than just calculating with widened precision.
At least I would like to see a performance comparison.
There was a problem hiding this comment.
Can someone help me find good test cases for performance testing for these methods?
Since I am new to Julia and do not have much experience with performance testing. I am unable to construct tests for which gcd would take more iterations.
I tried running a few tests for the gcd method, which take like 0.000001 seconds to complete.
There was a problem hiding this comment.
For Euclidean GCD algorithm, the worst case occurs when arguments are consecutive Fibonacci numbers.
Using this, I tested gcd(f1 + 0im, f2+ 0im), for f1 = 1836311903, f2= 2971215073, since these are the largest Int64 arguments for which gcd does not overflow.
Benchmarking with @time gives 0.000009 seconds (3 allocations: 96 bytes) for the current overflow-checked implementation.
I compared it against another implementation, which calculates gcd by widening arguments into Complex{BigInt} first. For such a complex gcd method, @time gives around 0.000150 seconds (2.23 k allocations: 46.305 KiB).
Hence, we can safely conclude that for arguments below up to f1, f2, calculating gcd without widening into Complex{BigInt} might be a better decision.
If someone finds better benchmarks, I am happy to try them as well.
…for complex integers. Added docstrings for div, rem, divrem methods.
… Added tests for it.
…for complex integers. Improved docstrings for div, rem, divrem methods.
|
Please fix this ambiguity error: |
This comment has been minimized.
This comment has been minimized.
|
Hi everyone, |
|
@simonbyrne, @stevengj, any chance either of you could review this and merge if it looks good? |
|
Hi everyone, |
| resulting in the return of a negative value of imaginary part. This overflow | ||
| occurs only when `conj` is applied to the minimum representable value of a | ||
| signed integer. That is, when `x == typemin(typeof(x))`, | ||
| `imag(conj(x * im)) == x < 0`, not `-x` as might be expected. |
There was a problem hiding this comment.
This warning seems unnecessary. We don't have similar warnings for - or other basic arithmetic operators.
| 1 - 3im | ||
| julia> conj(1 + typemin(Int64)im) | ||
| 1 - 9223372036854775808im |
There was a problem hiding this comment.
same … this doesn't seem the place to talk about integer overflow
There was a problem hiding this comment.
I just added it because the expected behaviour would be 1 - typemin(Int64) im.
65 participants
No description provided.