November 24, 2013
On Sunday, 24 November 2013 at 18:37:48 UTC, Andrei Alexandrescu wrote:
> But that originates as a call to multiplication between two Complex numbers. Can't the problem be addressed at that level?

Don't see why not, but it's going to be an unpleasant mess of if/else unless anyone can think up any clever tricks to avoid it.

BTW is it true that IEEE standards define 0 * inf to be nan? It's counter-intuitive mathematically and I can't find a reference to that effect.
November 24, 2013
On 11/24/13 10:56 AM, Shammah Chancellor wrote:
> On 2013-11-24 18:37:51 +0000, Andrei Alexandrescu said:
>
>> But that originates as a call to multiplication between two Complex
>> numbers. Can't the problem be addressed at that level?
>>
>> Andrei
>
> I don't believe so because IEEE floats define inf*0 to be NaN.   You
> would have to check to see if rhs.re == 0 || lhs.re == 0  and then just
> return zero.  Somewhat unfortunate.  You really do need an imaginary
> type for reasons specified in the original page here:
>
> http://digitalmars.com/d/1.0/cppcomplex.html
>
> -Shammah

Sure. We ain't above defining an imaginary type!

Andrei

November 24, 2013
On 11/24/13 3:24 PM, Joseph Rushton Wakeling wrote:
> On Sunday, 24 November 2013 at 18:37:48 UTC, Andrei Alexandrescu wrote:
>> But that originates as a call to multiplication between two Complex
>> numbers. Can't the problem be addressed at that level?
>
> Don't see why not, but it's going to be an unpleasant mess of if/else
> unless anyone can think up any clever tricks to avoid it.

How many special cases are out there?

> BTW is it true that IEEE standards define 0 * inf to be nan? It's
> counter-intuitive mathematically and I can't find a reference to that
> effect.

It sort of does. If you multiply something that goes to 0 with something that goes to infinity it could go either way.


Andrei

November 24, 2013
On 2013-11-24 23:24:18 +0000, Joseph Rushton Wakeling said:

> BTW is it true that IEEE standards define 0 * inf to be nan? It's counter-intuitive mathematically and I can't find a reference to that effect.

I couldn't find a reference either, but that's certainly how it is handled.   I doubt there is a bug in so many compilers regarding that.  Although, my friend Fred Tydeman makes his living finding these sorts of bugs with compilers: http://www.tybor.com/

November 25, 2013
On 25/11/13 00:35, Andrei Alexandrescu wrote:
> How many special cases are out there?

Well, if you have two complex numbers, z1 and z2, then

   (z1 * z2).re = (z1.re * z2.re) - (z1.im * z2.im)

and

   (z1 * z2).im = (z1.im * z2.re) + (z2.re * z1.im)

So you have to do if's for all four of z1.im, z2.re, z2.re and z1.im, and then you have to decide whether you override the apparent IEEE default just for the case of (re * im) or whether you do it for everything.

I mean, it feels a bit weird if you allow 0 * inf = 0 when it's real part times imaginary part, but not when it's real part times real part.

Do the IEEE standards cover implementation of complex numbers?  I confess complete ignorance here (and the Wikipedia page wasn't any help).

> It sort of does. If you multiply something that goes to 0 with something that
> goes to infinity it could go either way.

I'm not sure I follow.  I mean, if you have two sequences {x_i} --> 0 and {y_i} --> inf, then sure, the sequence of the product {x_i * y_i} can go either way. But if you just think of 0 as a number, then

     0 * inf = 0 * lim{x --> inf} x
             = lim{x --> inf} (0 * x)
             = lim{x --> inf} 0
             = 0

Or am I missing something about how FP numbers are implemented that either makes it convenient to define 0 * inf as not-a-number or that means that there are ambiguities that don't exist for "real" real numbers?
November 25, 2013
On 11/25/13 12:18 AM, Joseph Rushton Wakeling wrote:
> On 25/11/13 00:35, Andrei Alexandrescu wrote:
>> How many special cases are out there?
>
> Well, if you have two complex numbers, z1 and z2, then
>
>     (z1 * z2).re = (z1.re * z2.re) - (z1.im * z2.im)
>
> and
>
>     (z1 * z2).im = (z1.im * z2.re) + (z2.re * z1.im)
>
> So you have to do if's for all four of z1.im, z2.re, z2.re and z1.im,
> and then you have to decide whether you override the apparent IEEE
> default just for the case of (re * im) or whether you do it for everything.
>
> I mean, it feels a bit weird if you allow 0 * inf = 0 when it's real
> part times imaginary part, but not when it's real part times real part.

Doesn't sound all that bad to me. After all the built-in complex must be doing something similar. Of course if a separate imaginary type helps this and other cases, we should define it.

>> It sort of does. If you multiply something that goes to 0 with
>> something that
>> goes to infinity it could go either way.
>
> I'm not sure I follow.  I mean, if you have two sequences {x_i} --> 0
> and {y_i} --> inf, then sure, the sequence of the product {x_i * y_i}
> can go either way. But if you just think of 0 as a number, then
>
>       0 * inf = 0 * lim{x --> inf} x
>               = lim{x --> inf} (0 * x)
>               = lim{x --> inf} 0
>               = 0

Heh, no need for the expansion :o). Zero is zero.

> Or am I missing something about how FP numbers are implemented that
> either makes it convenient to define 0 * inf as not-a-number or that
> means that there are ambiguities that don't exist for "real" real numbers?

Well 0 in FP can always be considered a number so small, 0 was the closest representation. But I'm just speculating as to whether that's the reason for the choice.


Andrei

November 25, 2013
On 25/11/13 18:51, Andrei Alexandrescu wrote:
> Doesn't sound all that bad to me. After all the built-in complex must be doing
> something similar. Of course if a separate imaginary type helps this and other
> cases, we should define it.

Well, if you want it I'm happy to write the patch.  It's just I'm not sure that what is happening with std.complex is actually wrong if it's to be considered correct that 0 * inf = nan.
November 26, 2013
On 11/25/13 10:37 AM, Joseph Rushton Wakeling wrote:
> On 25/11/13 18:51, Andrei Alexandrescu wrote:
>> Doesn't sound all that bad to me. After all the built-in complex must
>> be doing
>> something similar. Of course if a separate imaginary type helps this
>> and other
>> cases, we should define it.
>
> Well, if you want it I'm happy to write the patch.  It's just I'm not
> sure that what is happening with std.complex is actually wrong if it's
> to be considered correct that 0 * inf = nan.

If you got the blessing of some complex number expert, that would be great.

Andrei
November 26, 2013
"Dmitry Olshansky" <dmitry.olsh@gmail.com> wrote in message news:l6tfm8$2hnj$1@digitalmars.com...
>>
>
> Can't it just check for the real part being exactly zero and special- case multiplication for that?
>

There is no such thing as exactly zero in floating point.  Only -0 and +0.


November 26, 2013
On 23/11/13 08:43, Ali Çehreli wrote:
> import std.stdio;
> import std.complex;
>
> void main()
> {
>      writeln(complex(1.0L, -real.infinity) * complex(0.0, 1.0L));
>      writeln((1L - ireal.infinity) * 1i);
> }
>
>
> The output:
>
> inf-nani    <-- "incorrect" according to the quoted page
> inf+1i      <-- correct

But, still operating with builtins,

    writeln((1L - ireal.infinity) * (0 + 1i));

and you get again

    inf-nani

Basically, your nice result with (1L - ireal.infinity) * 1i comes about because in this case, you're not multiplying two complex numbers, but one complex and one imaginary.  In the latter case, there's no 0 to multiply by infinity.