Kahan strongly prefers an Imaginary type .. should Julia?

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
10 messages Options
Reply | Threaded
Open this post in threaded view
|

Kahan strongly prefers an Imaginary type .. should Julia?

Jeffrey Sarnoff
I was rereading How Java's Floating-Point Hurts Everyone Everywhere, and seeing how Julia handled some of the difficulties therein (pp 11-15; this from pg 15).  

A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of ± 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + ı·y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type.  


julia> c = Complex(-1,0);
julia> z = sqrt(conj(c)) - conj(sqrt(c))
0.0 + 2.0im                   # should be 0.0+0.0im
julia> abs(z), z^2
2.0, -4.0+0.0im               # should be 0.0, 0.0+0.0im

Is this considered acceptable?

 


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Stefan Karpinski
We considered this at some point. I don't really recall why we decided against it.

On Thu, Feb 25, 2016 at 8:57 AM, Jeffrey Sarnoff <[hidden email]> wrote:
I was rereading How Java's Floating-Point Hurts Everyone Everywhere, and seeing how Julia handled some of the difficulties therein (pp 11-15; this from pg 15).  

A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of ± 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + ı·y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type.  


julia> c = Complex(-1,0);
julia> z = sqrt(conj(c)) - conj(sqrt(c))
0.0 + 2.0im                   # should be 0.0+0.0im
julia> abs(z), z^2
2.0, -4.0+0.0im               # should be 0.0, 0.0+0.0im

Is this considered acceptable?

 



Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

feza
In reply to this post by Jeffrey Sarnoff
Why does this work for c = Complex(-1.0,0.0) but not when using integers?

On Thursday, February 25, 2016 at 8:57:03 AM UTC-5, Jeffrey Sarnoff wrote:
I was rereading <a href="http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf" target="_blank" rel="nofollow" onmousedown="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;" onclick="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;">How Java's Floating-Point Hurts Everyone Everywhere, and seeing how Julia handled some of the difficulties therein (pp 11-15; this from pg 15).  

A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of ± 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + ı·y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type.  


julia> c = Complex(-1,0);
julia> z = sqrt(conj(c)) - conj(sqrt(c))
0.0 + 2.0im                   # should be 0.0+0.0im
julia> abs(z), z^2
2.0, -4.0+0.0im               # should be 0.0, 0.0+0.0im

Is this considered acceptable?

 


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Simon Byrne
Note that Complex(1,0) is a Complex{Int}, which doesn't have signed zeros. Kahan is referring to floating point pairs, which do have signed zeros, and as freaz points out, do actually behave correctly, so I'm not sure how this is an argument in favour of an imaginary type.

I think a better argument in favour of an imaginary type might be something like
https://github.com/JuliaLang/julia/issues/10000

-Simon


On Thursday, 25 February 2016 18:47:28 UTC, freaz wrote:
Why does this work for c = Complex(-1.0,0.0) but not when using integers?

On Thursday, February 25, 2016 at 8:57:03 AM UTC-5, Jeffrey Sarnoff wrote:
I was rereading <a href="http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf" rel="nofollow" target="_blank" onmousedown="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;" onclick="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;">How Java's Floating-Point Hurts Everyone Everywhere, and seeing how Julia handled some of the difficulties therein (pp 11-15; this from pg 15).  

A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of ± 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + ı·y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type.  


julia> c = Complex(-1,0);
julia> z = sqrt(conj(c)) - conj(sqrt(c))
0.0 + 2.0im                   # should be 0.0+0.0im
julia> abs(z), z^2
2.0, -4.0+0.0im               # should be 0.0, 0.0+0.0im

Is this considered acceptable?

 


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Simon Byrne
Perhaps a better argument is something like:

julia> (1 - Inf*im)*4*im

Inf + 4.0im


julia> (1 - Inf*im)*(4*im)

Inf + NaN*im


(adapted from http://www.digitalmars.com/d/1.0/cppcomplex.html)



On Friday, 26 February 2016 09:44:31 UTC, Simon Byrne wrote:
Note that Complex(1,0) is a Complex{Int}, which doesn't have signed zeros. Kahan is referring to floating point pairs, which do have signed zeros, and as freaz points out, do actually behave correctly, so I'm not sure how this is an argument in favour of an imaginary type.

I think a better argument in favour of an imaginary type might be something like
<a href="https://github.com/JuliaLang/julia/issues/10000" target="_blank" rel="nofollow" onmousedown="this.href=&#39;https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2FJuliaLang%2Fjulia%2Fissues%2F10000\46sa\75D\46sntz\0751\46usg\75AFQjCNF8PhGTW69Aw5_RWHYgaepOHruFrQ&#39;;return true;" onclick="this.href=&#39;https://www.google.com/url?q\75https%3A%2F%2Fgithub.com%2FJuliaLang%2Fjulia%2Fissues%2F10000\46sa\75D\46sntz\0751\46usg\75AFQjCNF8PhGTW69Aw5_RWHYgaepOHruFrQ&#39;;return true;">https://github.com/JuliaLang/julia/issues/10000

-Simon


On Thursday, 25 February 2016 18:47:28 UTC, freaz wrote:
Why does this work for c = Complex(-1.0,0.0) but not when using integers?

On Thursday, February 25, 2016 at 8:57:03 AM UTC-5, Jeffrey Sarnoff wrote:
I was rereading <a href="http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf" rel="nofollow" target="_blank" onmousedown="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;" onclick="this.href=&#39;http://www.google.com/url?q\75http%3A%2F%2Fwww.cs.berkeley.edu%2F~wkahan%2FJAVAhurt.pdf\46sa\75D\46sntz\0751\46usg\75AFQjCNHCZU2FWqL1VnAHzub-hSpjw-A8AQ&#39;;return true;">How Java's Floating-Point Hurts Everyone Everywhere, and seeing how Julia handled some of the difficulties therein (pp 11-15; this from pg 15).  

A streamline goes astray when the complex functions SQRT and LOG are implemented, as is necessary in Fortran and in libraries currently distributed with C/C++ compilers, in a way that disregards the sign of ± 0.0 in IEEE 754 arithmetic and consequently violates identities like SQRT( CONJ( Z ) ) = CONJ( SQRT( Z ) ) and LOG( CONJ( Z ) ) = CONJ( LOG( Z ) ) whenever the COMPLEX variable Z takes negative real values. Such anomalies are unavoidable if Complex Arithmetic operates on pairs (x, y) instead of notional sums x + ı·y of real and imaginary variables. The language of pairs is incorrect for Complex Arithmetic; it needs the Imaginary type.  


julia> c = Complex(-1,0);
julia> z = sqrt(conj(c)) - conj(sqrt(c))
0.0 + 2.0im                   # should be 0.0+0.0im
julia> abs(z), z^2
2.0, -4.0+0.0im               # should be 0.0, 0.0+0.0im

Is this considered acceptable?

 


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Andreas Lobinger
In reply to this post by Stefan Karpinski
On Thursday, February 25, 2016 at 3:14:07 PM UTC+1, Stefan Karpinski wrote:
We considered this at some point. I don't really recall why we decided against it.

Matlab?

>> c = -1+0i;
>> z = sqrt(conj(c)) - conj(sqrt(c))

z =

   0.0000 + 2.0000i

>> abs(z), z^2

ans =

     2


ans =

    -4


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Jason Riedy
In reply to this post by Stefan Karpinski
And Stefan Karpinski writes:
> We considered this at some point. I don't really recall why we
> decided against it.

Apologies for pinging an old thread, but people may be interested
in the tech report by Prof. Kahan and Jim Thomas justifying an
explicit imaginary type:
  https://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6127.html
I'm always amused(?) how functional language folks skip this topic.

Jim Thomas's name may be familiar from the C standard for ongoing
work related to floating point...  (He often keeps the 754 working
group from wandering off into the weeds, too.)

You likely decided against it because there's no sensible way to
go *back* to imaginary when the real part becomes zero.  That
would require distinguishing exact zeros, etc.  Given little use
of explicitly imaginary numbers outside of programming
homework...  The interaction with elementary functions is perhaps
one of the more motivational real-world uses where people copy
formulas from text books / wikipedia.

Imaginary{} building into Complex{} could be a nice, very
well-defined example for Julia's type system, however.  Even
better if someone can demonstrate it outside of the base.

Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Stefan Karpinski
Some history here:
So basically:
  • We originally had an ImaginaryUnit type of which `im` was the singleton instance.
  • People (reasonably) expected functions to work on `im` but they didn't because it wasn't a real numerical type.
  • We had one PR for an Imaginary type and another PR for making `im = Complex(false,true)` and tweaking the promotion rules for Bools to make things work out mostly correctly.
  • The PR for making `im = Complex(false,true)` won out.
I have some doubts about whether this was the right choice:
  • This approach doesn't address all the complex number issues.
  • The special promotion rules for Bools are annoying when adding promotion rules since they always cause ambiguities.
  • If we had an AbstractComplex type and had Imaginary <: AbstractComplex, then most complex functions would just work for imaginary values, although some of them could of course have much more efficient and simple special case implementations. This would also allow a PolarComplex type to be provided in a package.

On Wed, Mar 23, 2016 at 4:34 PM, Jason Riedy <[hidden email]> wrote:
And Stefan Karpinski writes:
> We considered this at some point. I don't really recall why we
> decided against it.

Apologies for pinging an old thread, but people may be interested
in the tech report by Prof. Kahan and Jim Thomas justifying an
explicit imaginary type:
  https://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6127.html
I'm always amused(?) how functional language folks skip this topic.

Jim Thomas's name may be familiar from the C standard for ongoing
work related to floating point...  (He often keeps the 754 working
group from wandering off into the weeds, too.)

You likely decided against it because there's no sensible way to
go *back* to imaginary when the real part becomes zero.  That
would require distinguishing exact zeros, etc.  Given little use
of explicitly imaginary numbers outside of programming
homework...  The interaction with elementary functions is perhaps
one of the more motivational real-world uses where people copy
formulas from text books / wikipedia.

Imaginary{} building into Complex{} could be a nice, very
well-defined example for Julia's type system, however.  Even
better if someone can demonstrate it outside of the base.


Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Erik Schnetter
Could this be a GSOC project?

-erik

On Wed, Mar 23, 2016 at 5:21 PM, Stefan Karpinski <[hidden email]> wrote:

> Some history here:
>
> https://github.com/JuliaLang/julia/pull/4922
> https://github.com/JuliaLang/julia/issues/5295
> https://github.com/JuliaLang/julia/pull/5313
> https://github.com/JuliaLang/julia/pull/5468
>
> So basically:
>
> We originally had an ImaginaryUnit type of which `im` was the singleton
> instance.
> People (reasonably) expected functions to work on `im` but they didn't
> because it wasn't a real numerical type.
> We had one PR for an Imaginary type and another PR for making `im =
> Complex(false,true)` and tweaking the promotion rules for Bools to make
> things work out mostly correctly.
> The PR for making `im = Complex(false,true)` won out.
>
> I have some doubts about whether this was the right choice:
>
> This approach doesn't address all the complex number issues.
> The special promotion rules for Bools are annoying when adding promotion
> rules since they always cause ambiguities.
> If we had an AbstractComplex type and had Imaginary <: AbstractComplex, then
> most complex functions would just work for imaginary values, although some
> of them could of course have much more efficient and simple special case
> implementations. This would also allow a PolarComplex type to be provided in
> a package.
>
>
> On Wed, Mar 23, 2016 at 4:34 PM, Jason Riedy <[hidden email]>
> wrote:
>>
>> And Stefan Karpinski writes:
>> > We considered this at some point. I don't really recall why we
>> > decided against it.
>>
>> Apologies for pinging an old thread, but people may be interested
>> in the tech report by Prof. Kahan and Jim Thomas justifying an
>> explicit imaginary type:
>>   https://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6127.html
>> I'm always amused(?) how functional language folks skip this topic.
>>
>> Jim Thomas's name may be familiar from the C standard for ongoing
>> work related to floating point...  (He often keeps the 754 working
>> group from wandering off into the weeds, too.)
>>
>> You likely decided against it because there's no sensible way to
>> go *back* to imaginary when the real part becomes zero.  That
>> would require distinguishing exact zeros, etc.  Given little use
>> of explicitly imaginary numbers outside of programming
>> homework...  The interaction with elementary functions is perhaps
>> one of the more motivational real-world uses where people copy
>> formulas from text books / wikipedia.
>>
>> Imaginary{} building into Complex{} could be a nice, very
>> well-defined example for Julia's type system, however.  Even
>> better if someone can demonstrate it outside of the base.
>>
>



--
Erik Schnetter <[hidden email]>
http://www.perimeterinstitute.ca/personal/eschnetter/
Reply | Threaded
Open this post in threaded view
|

Re: Kahan strongly prefers an Imaginary type .. should Julia?

Stefan Karpinski
Perhaps, although I had already implemented most of the basic work in the above pull request. But it didn't implement abstract complex types or specialized methods for functions of purely imaginary types. Or a polar complex type.

On Wednesday, March 23, 2016, Erik Schnetter <[hidden email]> wrote:
Could this be a GSOC project?

-erik

On Wed, Mar 23, 2016 at 5:21 PM, Stefan Karpinski <<a href="javascript:;" onclick="_e(event, &#39;cvml&#39;, &#39;stefan@karpinski.org&#39;)">stefan@...> wrote:
> Some history here:
>
> https://github.com/JuliaLang/julia/pull/4922
> https://github.com/JuliaLang/julia/issues/5295
> https://github.com/JuliaLang/julia/pull/5313
> https://github.com/JuliaLang/julia/pull/5468
>
> So basically:
>
> We originally had an ImaginaryUnit type of which `im` was the singleton
> instance.
> People (reasonably) expected functions to work on `im` but they didn't
> because it wasn't a real numerical type.
> We had one PR for an Imaginary type and another PR for making `im =
> Complex(false,true)` and tweaking the promotion rules for Bools to make
> things work out mostly correctly.
> The PR for making `im = Complex(false,true)` won out.
>
> I have some doubts about whether this was the right choice:
>
> This approach doesn't address all the complex number issues.
> The special promotion rules for Bools are annoying when adding promotion
> rules since they always cause ambiguities.
> If we had an AbstractComplex type and had Imaginary <: AbstractComplex, then
> most complex functions would just work for imaginary values, although some
> of them could of course have much more efficient and simple special case
> implementations. This would also allow a PolarComplex type to be provided in
> a package.
>
>
> On Wed, Mar 23, 2016 at 4:34 PM, Jason Riedy <<a href="javascript:;" onclick="_e(event, &#39;cvml&#39;, &#39;jason@lovesgoodfood.com&#39;)">jason@...>
> wrote:
>>
>> And Stefan Karpinski writes:
>> > We considered this at some point. I don't really recall why we
>> > decided against it.
>>
>> Apologies for pinging an old thread, but people may be interested
>> in the tech report by Prof. Kahan and Jim Thomas justifying an
>> explicit imaginary type:
>>   https://www.eecs.berkeley.edu/Pubs/TechRpts/1992/6127.html
>> I'm always amused(?) how functional language folks skip this topic.
>>
>> Jim Thomas's name may be familiar from the C standard for ongoing
>> work related to floating point...  (He often keeps the 754 working
>> group from wandering off into the weeds, too.)
>>
>> You likely decided against it because there's no sensible way to
>> go *back* to imaginary when the real part becomes zero.  That
>> would require distinguishing exact zeros, etc.  Given little use
>> of explicitly imaginary numbers outside of programming
>> homework...  The interaction with elementary functions is perhaps
>> one of the more motivational real-world uses where people copy
>> formulas from text books / wikipedia.
>>
>> Imaginary{} building into Complex{} could be a nice, very
>> well-defined example for Julia's type system, however.  Even
>> better if someone can demonstrate it outside of the base.
>>
>



--
Erik Schnetter <<a href="javascript:;" onclick="_e(event, &#39;cvml&#39;, &#39;schnetter@gmail.com&#39;)">schnetter@...>
http://www.perimeterinstitute.ca/personal/eschnetter/