www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - std.complex

reply "Craig Dillabaugh" <cdillaba cg.scs.carleton.ca> writes:
The complex type in std.complex restricts the real/imaginary
parts of the number to be float/double/real.  I am curious to
know if there is a reason why integral types are not permitted. I
checked the C++ equivalent and it does not have the same
requirement.

I mention this because some of my work is done with radar
satellite images. All pixels in such an images are stored as
complex numbers, but in all cases I am aware of they are stored a
short int values. Most software that operates on the images uses
complex<short> (most of it is C++).

Is there any reason why complex numbers in D's standard lib must
be of non-integral types?  I believe in C++ the type is optimized
for floating point values, but allows other types.
Nov 18 2013
next sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 11/18/13 5:44 PM, Craig Dillabaugh wrote:
 The complex type in std.complex restricts the real/imaginary
 parts of the number to be float/double/real.  I am curious to
 know if there is a reason why integral types are not permitted. I
 checked the C++ equivalent and it does not have the same
 requirement.

 I mention this because some of my work is done with radar
 satellite images. All pixels in such an images are stored as
 complex numbers, but in all cases I am aware of they are stored a
 short int values. Most software that operates on the images uses
 complex<short> (most of it is C++).

 Is there any reason why complex numbers in D's standard lib must
 be of non-integral types?  I believe in C++ the type is optimized
 for floating point values, but allows other types.

The simple reason is we couldn't find appropriate applications. If you make a good argument, we'll include integral types as well. Submit an enhancement request on bugzilla including your example and let's take it from there. Andrei
Nov 18 2013
parent reply Stewart Gordon <smjg_1998 yahoo.com> writes:
On 19/11/2013 02:03, Andrei Alexandrescu wrote:> On 11/18/13 5:44 PM, 
Craig Dillabaugh wrote:
<snip>
 Is there any reason why complex numbers in D's standard lib must
 be of non-integral types?  I believe in C++ the type is optimized
 for floating point values, but allows other types.

The simple reason is we couldn't find appropriate applications.

I don't understand. At the moment Complex appears to me to be type-agnostic - as long as a type supports the standard arithmetic operators and assignment of the value 0 to it, it will work. The only thing preventing it from working at the moment is this line struct Complex(T) if (isFloatingPoint!T) So why do you need an appropriate application in order not to have this arbitrary restriction? Or have I missed something? It isn't just integer types. Somebody might want to use complex with a library (fixed-point, arbitrary precision, decimal, etc.) numeric type. Fractal generators, for example, are likely to use this a lot. Or even more exotically, use Complex!(Complex!real) to implement hypercomplex numbers. Stewart.
Jan 01 2014
next sibling parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 1/1/14 3:12 PM, Lars T. Kyllingstad wrote:
 On Wednesday, 1 January 2014 at 12:29:35 UTC, Stewart Gordon wrote:
  [...]

 Or even more exotically, use Complex!(Complex!real) to implement
 hypercomplex numbers.

This is an extremely marginal use case. Currently Complex!(Complex!T) folds to Complex!T.

What was the motivation? Andrei
Jan 01 2014
prev sibling next sibling parent reply Stewart Gordon <smjg_1998 yahoo.com> writes:
On 01/01/2014 19:55, Joseph Rushton Wakeling wrote:
<snip>
 There are binary operations on complex numbers where the only sensible
 outcome seems to be non-integral real and imaginary parts. Addition,
 subtraction and multiplication are OK with integral types, but division
 really seems unpleasant to implement absent floating point,

Then why not just disable division if it's a non-float type, rather than preventing the whole complex template from being used with that type? This is like cutting off somebody's arm because they have a sore thumb. Moreover, we have no way in the general case of determining whether T is an integral type, a library float-esque type, or (for example) a Galois field type. So disabling it _just in case_ division doesn't work is crazy. There must be better ways to do it.
 exponentiation even more so.

Exponentiation by a non-negative integer is straightforward. So we should at least support this case for Gaussian integers. <snip>
 However, I think relaxing the template constraints like this would best
 be done in the context of a library float-esque type (e.g. BigFloat)
 being implemented in Phobos, which could then be used to provide both
 proof-of-concept and the primary test case.

What do you mean by "in the context of", exactly? Restricting it to some float-esque type that is in Phobos would still be overly restrictive.
 Or even more exotically, use Complex!(Complex!real) to implement
 hypercomplex numbers.

Perhaps best implemented as a type in its own right?

It would, but removing the restriction would simplify the implementation of Hypercomplex(T) by enabling it to be a wrapper for Complex!(Complex!T). Stewart.
Jan 02 2014
parent reply Stewart Gordon <smjg_1998 yahoo.com> writes:
On 02/01/2014 19:32, Joseph Rushton Wakeling wrote:
 On Thursday, 2 January 2014 at 18:12:56 UTC, Stewart Gordon wrote:
 Then why not just disable division if it's a non-float type, rather
 than preventing the whole complex template from being used with that
 type? This is like cutting off somebody's arm because they have a sore
 thumb.

Because that also seems to me to be an unpleasant option. A complex number implementation that fails to support ordinary arithmetic operations in all circumstances is pretty non-intuitive and will probably lead to unpleasant bugs in users' code.

The compiler rejecting the code is the most pleasant bug that's possible IMO.
 Moreover, we have no way in the general case of determining whether T
 is an integral type, a library float-esque type, or (for example) a
 Galois field type.  So disabling it _just in case_ division doesn't
 work is crazy.  There must be better ways to do it.

I don't follow your point here. We can constrain T however we see fit. The point isn't to have some perfect representation of every mathematical possibility, it's to have useful code that serves a good range of use-cases while being robust and easy to maintain.

Not being overly restrictive in what types you will allow it to be used with is an important part of serving a good range of use cases. <snip>
 What do you mean by "in the context of", exactly?  Restricting it to
 some float-esque type that is in Phobos would still be overly
 restrictive.

No, I mean that until you have at least one actual float-esque type to test with, it is probably unwise to relax the template constraints that currently mandate built-in FP types.

I'm sure we have a small handful already. We just need to find them. For instance, given the time I could probably dig up my rational number implementation and update it to current D.
 It would, but removing the restriction would simplify the
 implementation of Hypercomplex(T) by enabling it to be a wrapper for
 Complex!(Complex!T).

And complicate the implementation of Complex itself, for the sake of a likely very marginal special interest that could be supported quite well by an independent type.

How would it complicate the implementation? Removing the undocumented rule whereby Complex!(Complex!T) folds to Complex!T would be a slight simplification. Maybe the right course of action is to have a parameter to the template that suppresses the type restriction and the folding rule, so that people who want to use it on a type that might not work properly can do so. This would be a relatively small change. Stewart.
Jan 02 2014
parent reply Stewart Gordon <smjg_1998 yahoo.com> writes:
On 02/01/2014 21:10, Joseph Rushton Wakeling wrote:
 On Thursday, 2 January 2014 at 20:23:40 UTC, Stewart Gordon wrote:

 How would it complicate the implementation?  Removing the undocumented
 rule whereby Complex!(Complex!T) folds to Complex!T would be a slight
 simplification.

It would involve non-trivial revisions to the internals of Complex.

Please be specific. <snip>
 Maybe the right course of action is to have a parameter to the
 template that suppresses the type restriction and the folding rule, so
 that people who want to use it on a type that might not work properly
 can do so.  This would be a relatively small change.

I do not know how familiar you are with the internals of std.complex.Complex but I think it would be a good idea to go through them in some detail before asserting that changes like this would be "relatively small".

Oh yes, some of the templated functions that take or return a complex would need this extra parameter adding. Another way it could be done is to have a Complex template that implements these rules and which programmers would normally use, and have this aliasing ComplexImpl which actually provides the implementation and which programmers can use directly if they want to bypass the restrictions. The difficulty is documenting it in a way that would make sense to normal users.... Stewart.
Jan 02 2014
parent reply Stewart Gordon <smjg_1998 yahoo.com> writes:
On 03/01/2014 00:03, Joseph Rushton Wakeling wrote:
 On 03/01/14 00:33, Stewart Gordon wrote:
 Please be specific.

You know, I'm starting to find your tone irritating. You are the one who's asking for functionality that goes beyond any Complex implementation that I'm aware of in any other language, and claiming that these things would be trivial to implement.

I wasn't asking for it to go beyond the existing complex implementation or any other. I was proposing that the arbitrary restriction be removed so that the implementation we already have would work on them. As I said originally:
 I don't understand. At the moment Complex appears to me to be
 type-agnostic - as long as a type supports the standard arithmetic
 operators and assignment of the value 0 to it, it will work. The
 only thing preventing it from working at the moment is this line

     struct Complex(T)  if (isFloatingPoint!T)

 So why do you need an appropriate application in order not to have
 this arbitrary restriction? Or have I missed something?

OK, so division has now been mentioned. And you have now mentioned the use of FPTemporary. However....
 I would expect a person who claims with confidence that something is
 trivial, to actually know the internals of the code well enough to
 understand what parts of it would need to be modified. On the basis
 of what you've written, I have no reason to believe that you do.

I had read the code before I made that comment, and from reading it I did believe at the time that the implementation was ready for types other than float, real and double. So the use of FPTemporary is something I missed, but it's taken you this long to point it out to me. However....
 There are numerous places inside current std.complex.Complex where
 temporary values are used mid-calculation.  Those are all of type
 FPTemporary (which in practice means real).So, to handle library
 types (whether library floating-point types such as a BigFloat
 implementation, or a Complex!T so as to support hypercomplex numbers)
 you'd either have to special-case those functions or you'd have to
 provide an alternative Temporary template to handle the temporary
 internal values in the different cases.

FPTemporary is a template. At the moment it's defined only for the built-in floating point types. So what we really need is to define FPTemporary for other types. For int and smaller integral types, we can define it as real just as we do for float/double/real. Whether it's adequate for long would be platform-dependent. For other types, I suppose we can reference a property of the type, or just use the type itself if such a property isn't present. One possible idea (untested): ----- template FPTemporary(F) if (is(typeof(-F.init + F.init * F.init - F.init))) { static if (isFloatingPoint!F || isIntegral!F) { alias real FPTemporary; } else static if (is(F.FPTemporary)) { alias F.FPTemporary FPTemporary; } else { alias F FPTemporary; } } ----- Of course, this isn't a completely general solution, and I can see now that whatever type we use would need to have trigonometric functions in order to support exponentiation. So unless we conditionalise the inclusion of this operation....
 You'd also need to address questions of closure under operations
 (already an issue for the Imaginary type), and precision-related issues
 -- see e.g. the remarks by Craig Dillabaugh and Ola Fosheim Grøstad
 elsewhere in this thread.

Oh yes, and it would be crazy to try and make it work for unsigned integer types. But even if we were to resolve the issues with FPTemporary and that, it would still fall under my earlier suggestion of making it so that if people want to use Complex on an unsupported type then they can explicitly suppress the type restriction, but should understand that it might not work properly. <snip>
 I don't personally see any rationale for implementing hypercomplex
 numbers as a specialization of Complex given that they can be just as
 well implemented as a type in their own right.

Really, I was just thinking that somebody who wants a quick-and-dirty hypercomplex number implementation for some app might try to do it that way. Stewart.
Jan 03 2014
parent Stewart Gordon <smjg_1998 yahoo.com> writes:
On 03/01/2014 17:04, Joseph Rushton Wakeling wrote:
 On 03/01/14 14:32, Stewart Gordon wrote:
 I wasn't asking for it to go beyond the existing complex
 implementation or any other. I was proposing that the arbitrary
 restriction be removed so that the implementation we already have
 would work on them.

Yes, but it isn't an arbitrary restriction. Template constraints are fundamentally a promise to users about what can be expected to work. Integral types, or library types, won't work without significant modifications to the internals of the code. It would be a false promise to relax those constraints.

But I believed the restriction to be arbitrary at the time I made my original point, hence my point. <snip>
 Yes, it ought to be possible to redefine FPTemporary (or define an
 alternative) to determine proper internal temporaries for any
 "float-esque" case.  I was toying with something along the lines of,

      template FPTemporary(F)
          if (isNumeric!F || isFloatLike!F)
      {
          alias typeof(real.init * F.init) FPTemporary;
      }

 ... where isFloatLike would test for appropriate floating-point-like
 properties of F -- although this is probably far too simplistic.

How can isFloatLike be implemented? And how can we test for bigint or Galois field types?
 E.g. how do you handle the case of a float-like library type
 implemented as a class, not a struct?

What's the difficulty here? <snip>
 It doesn't matter if you document it as "This might not work", by
 providing the option you are still essentially saying, "This is an OK
 way to use the type."

 I think that's essentially an encouragement of bad code and a violation
 of the design principle that the easy thing to do should be the right
 thing to do.

There are already violations of this principle in D, such as being able to cast away const.
 Really, I was just thinking that somebody who wants a quick-and-dirty
 hypercomplex number implementation for some app might try to do it
 that way.

I understand that, but quick-and-dirty solutions are often bad ones, and in this case, it just wouldn't work given the internals of Complex.

Addition, subtraction and multiplication would work. So the programmer could just copy the code and reimplement division and exponentiation so that they work (or just get rid of them if they aren't needed). Stewart.
Jan 03 2014
prev sibling next sibling parent Stewart Gordon <smjg_1998 yahoo.com> writes:
On 01/01/2014 23:12, Lars T. Kyllingstad wrote:
<snip>
 /*  Makes Complex!(Complex!T) fold to Complex!T.

      The rationale for this is that just like the real line is a
      subspace of the complex plane, the complex plane is a subspace
      of itself. Example of usage:

This doesn't seem to make sense - _any_ set is a subspace of itself. But I suppose one way to look at it is that Complex!T means the minimal vector space containing x and i*x for all x in T.
      ---
      Complex!T addI(T)(T x)
      {
          return x + Complex!T(0.0, 1.0);
      }
      ---

So the point is to enable templated functions to accept a real or complex argument and return a complex number.
      The above will work if T is both real and complex.
 */

T is "both real and complex"? How is this possible? Stewart.
Jan 02 2014
prev sibling parent Stewart Gordon <smjg_1998 yahoo.com> writes:
On 02/01/2014 15:38, Simen Kjærås wrote:
 On 2014-01-01 12:29, Stewart Gordon wrote:
 Or even more exotically, use Complex!(Complex!real) to implement
 hypercomplex numbers.

For a more generic solution to this, see Cayley-Dickson construction[1] and my implementation of such:

Hypercomplex numbers are outside the scope of Cayley-Dickson construction. The latter defines quaternions, octonions and so on. Though it is confusing as there are at least two definitions of "hypercomplex number": (a) a element of any number system that extends the complex numbers (b) an element of a particular such number system, which is what I was using it to mean. The hypercomplex numbers to which I was referring are as Fractint uses the term, and briefly described at http://mathworld.wolfram.com/HypercomplexNumber.html ij = k, i^2 = j^2 = -1, k = 1. Multiplication is commutative. But the distinctive thing about these is that they can be represented as an ordered pair of complex numbers, and then the complex multiplication formula works on them, hence my suggestion.
 https://github.com/Biotronic/Collectanea/blob/master/biotronic/CayleyDickson2.d


 It's nowhere near finished, has some bugs, and I haven't worked on it
 for at least a year, but it's an interesting way to create complex,
 hypercomplex, dual, and split-complex numbers (and combinations thereof).

So it's a generalisation of Cayley-Dickson construction to produce a range of 2^n-dimensional algebras over the reals. I'll have to look at it in more detail when I've more time. Stewart.
Jan 02 2014
prev sibling next sibling parent "Craig Dillabaugh" <cdillaba cg.scs.carleton.ca> writes:
On Tuesday, 19 November 2013 at 02:03:05 UTC, Andrei Alexandrescu
wrote:
 On 11/18/13 5:44 PM, Craig Dillabaugh wrote:
 The complex type in std.complex restricts the real/imaginary
 parts of the number to be float/double/real.  I am curious to
 know if there is a reason why integral types are not 
 permitted. I
 checked the C++ equivalent and it does not have the same
 requirement.

 I mention this because some of my work is done with radar
 satellite images. All pixels in such an images are stored as
 complex numbers, but in all cases I am aware of they are 
 stored a
 short int values. Most software that operates on the images 
 uses
 complex<short> (most of it is C++).

 Is there any reason why complex numbers in D's standard lib 
 must
 be of non-integral types?  I believe in C++ the type is 
 optimized
 for floating point values, but allows other types.

The simple reason is we couldn't find appropriate applications. If you make a good argument, we'll include integral types as well. Submit an enhancement request on bugzilla including your example and let's take it from there. Andrei

Regards, Craig
Nov 18 2013
prev sibling next sibling parent reply Shammah Chancellor <anonymous coward.com> writes:
//Hijack

http://digitalmars.com/d/1.0/cppcomplex.html

 • Consider the formula (1 - infinity*i) * i which should produce 
 (infinity + i). However, if instead the second factor is (0 + i) rather 
 than just i, the result is (infinity + NaN*i), a spurious NaN was 
 generated.
 • A distinct imaginary type preserves the sign of 0, necessary for 
 calculations involving branch cuts.

Is this stuff no longer an issue? -Shammah
Nov 22 2013
next sibling parent reply =?UTF-8?B?QWxpIMOHZWhyZWxp?= <acehreli yahoo.com> writes:
On 11/22/2013 09:22 PM, Craig Dillabaugh wrote:

 On Saturday, 23 November 2013 at 04:37:19 UTC, Shammah Chancellor wrote:
 //Hijack

 http://digitalmars.com/d/1.0/cppcomplex.html

 • Consider the formula (1 - infinity*i) * i which should produce
 (infinity + i). However, if instead the second factor is (0 + i)
 rather than just i, the result is (infinity + NaN*i), a spurious NaN
 was generated.
 • A distinct imaginary type preserves the sign of 0, necessary for
 calculations involving branch cuts.

Is this stuff no longer an issue? -Shammah

I believe D used to have builtin complex types, back in the old days. They have been removed (deprecated?)

It still compiles.
 and replaced by the library type
 std.complex.

 At least that is my understanding.

And it makes what Shammah Chancellor quoted even more interesting. cdouble and idouble still work correctly but std.complex produces "incorrect" result: 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
 Craig

Ali
Nov 22 2013
parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 11/24/13 9:54 AM, Joseph Rushton Wakeling wrote:
 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

It's because 0.0L * (-real.infinity) evaluates to nan.

Has this been submitted as a bug report? Andrei
Nov 24 2013
parent reply Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-24 18:00:45 +0000, Andrei Alexandrescu said:

 On 11/24/13 9:54 AM, Joseph Rushton Wakeling wrote:
 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

It's because 0.0L * (-real.infinity) evaluates to nan.

Has this been submitted as a bug report? Andrei

It's more a fundamental problem with a complex type in general. C++ has this issue as well. You need a purely imaginary type with the appropiate operations between Complex and Imaginary defined.
Nov 24 2013
next sibling parent reply Dmitry Olshansky <dmitry.olsh gmail.com> writes:
24-Nov-2013 22:03, Shammah Chancellor пишет:
 On 2013-11-24 18:00:45 +0000, Andrei Alexandrescu said:

 On 11/24/13 9:54 AM, Joseph Rushton Wakeling wrote:
 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

It's because 0.0L * (-real.infinity) evaluates to nan.

Has this been submitted as a bug report? Andrei

It's more a fundamental problem with a complex type in general. C++ has this issue as well. You need a purely imaginary type with the appropiate operations between Complex and Imaginary defined.

Can't it just check for the real part being exactly zero and special- case multiplication for that? -- Dmitry Olshansky
Nov 24 2013
parent reply "Daniel Murphy" <yebblies nospamgmail.com> writes:
"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.
Nov 25 2013
parent reply Dmitry Olshansky <dmitry.olsh gmail.com> writes:
26-Nov-2013 09:06, Daniel Murphy пишет:
 "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.

Well, let it be magnitude, "exactly" implies as good zero test as we need. I'm still of the opinion that a few predictable branches like `if(rhs.re.isZero)` won't kill it. Regardless, the point is that built-ins hardly help here at all as you may just as well get an "exactly" 0+1i for some complex computation, and the result won't suddenly change type to i{real,double,float}. -- Dmitry Olshansky
Nov 26 2013
parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 11/26/13 7:13 AM, Dmitry Olshansky wrote:
 26-Nov-2013 09:06, Daniel Murphy пишет:
 "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.

Well, let it be magnitude, "exactly" implies as good zero test as we need. I'm still of the opinion that a few predictable branches like `if(rhs.re.isZero)` won't kill it. Regardless, the point is that built-ins hardly help here at all as you may just as well get an "exactly" 0+1i for some complex computation, and the result won't suddenly change type to i{real,double,float}.

I think the problem is we currently don't have a means to distinguish between two cases: (a) "There is no real part here at all" (b) "There is a real part that happens to be zero, or really as close to zero as it can ever get in a discrete representation" It seems the two cases behave differently in a few corner cases, and that is reasonable. Andrei
Nov 26 2013
prev sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 11/24/13 10:03 AM, Shammah Chancellor wrote:
 On 2013-11-24 18:00:45 +0000, Andrei Alexandrescu said:

 On 11/24/13 9:54 AM, Joseph Rushton Wakeling wrote:
 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

It's because 0.0L * (-real.infinity) evaluates to nan.

Has this been submitted as a bug report? Andrei

It's more a fundamental problem with a complex type in general. C++ has this issue as well. You need a purely imaginary type with the appropiate operations between Complex and Imaginary defined.

But that originates as a call to multiplication between two Complex numbers. Can't the problem be addressed at that level? Andrei
Nov 24 2013
next sibling parent reply Shammah Chancellor <anonymous coward.com> writes:
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
Nov 24 2013
parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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
Nov 24 2013
prev sibling next sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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
Nov 24 2013
next sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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
Nov 25 2013
parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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
Nov 25 2013
next sibling parent reply Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-26 01:02:12 +0000, Andrei Alexandrescu said:

 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

With any luck, I'll be working on my PhD in computational physics using D toward the end of next year. I'll definitely be giving feedback :) As of right now, I have been out of the computational physics field for too long to give good feedback. I can't wait to make a version of the concurrency library that works across machines :D
Nov 26 2013
parent reply Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 12:52, Shammah Chancellor wrote:
 With any luck, I'll be working on my PhD in computational physics using D
toward
 the end of next year.    I'll definitely be giving feedback :) As of right now,
 I have been out of the computational physics field for too long to give good
 feedback.  I can't wait to make a version of the concurrency library that works
 across machines :D

Sounds very cool. What kind of stuff will you be working on?
Nov 26 2013
parent Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-26 12:55:26 +0000, Joseph Rushton Wakeling said:

 Sounds very cool.  What kind of stuff will you be working on?

Who knows. Haven't even finished my apps yet. The deadline is in a couple weeks :)
Nov 26 2013
prev sibling next sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this. Thanks, Andrei
Nov 26 2013
parent Shammah Chancellor <anonymous coward.com> writes:
On 2013-12-20 16:25:53 +0000, Joseph Rushton Wakeling said:

 On 26/11/13 17:28, Andrei Alexandrescu wrote:
 On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

Provisional code: https://github.com/D-Programming-Language/phobos/pull/1797 I also realized there was no explicit issue report, so created one: https://d.puremagic.com/issues/show_bug.cgi?id=11787

Awesome. Thank you! -Shammah
Dec 26 2013
prev sibling parent Walter Bright <newshound2 digitalmars.com> writes:
On 11/26/2013 3:22 AM, Joseph Rushton Wakeling wrote:
 (It also leaves std.complex' behaviour incompatible with std::complex and other
 library implementations, though I guess that's less of a concern so long as we
 think what std.complex does is right.)

Making it behave differently than others would be a disaster.
Jan 01 2014
prev sibling parent Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-27 09:14:04 +0000, Joseph Rushton Wakeling said:

 I'm sure you've heard that old anecdote of the professor back in the 
 1950s, or was it the 1920s, who, on hearing a student say "infinity", 
 said: "I won't have bad language in class!" :-)

Yes. I think that's part of the reason 0 • inf = NaN in IEEE. The values are taken to be limits of unknown functions. Thus 0 * inf is uncalculable without knowing those functions. There is no such thing as the value infinity.
Nov 27 2013
prev sibling next sibling parent Shammah Chancellor <anonymous coward.com> writes:
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/
Nov 24 2013
prev sibling parent "H. S. Teoh" <hsteoh quickfur.ath.cx> writes:
On Wed, Nov 27, 2013 at 11:29:55PM +0100, Joseph Rushton Wakeling wrote:
 On 26/11/13 17:28, Andrei Alexandrescu wrote:
On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
So, as other people have suggested, really the only thing we can
reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

I've started work, code is here: https://github.com/WebDrake/phobos/tree/imaginary Feedback, forks/pull requests, etc. welcome.

For the benefit of others: If you already have a Phobos fork on github, forking this repo won't do anything. What you need to do in that case is to add this repo as a remote, then fetch the 'imaginary' branch and check it out, like this: # Add this repo as a new remote named 'webdrake' git remote add webdrake https://github.com/WebDrake/phobos.git # Fetch the branch named 'imaginary' into webdrake/imaginary git fetch webdrake imaginary # Checkout the fetched branch into a local branch named # 'imaginary' (you can name this whatever you like, it doesn't # matter) git checkout -b imaginary webdrake/imaginary Now you can review / edit the code, push to origin (github), and then submit pull requests (make sure to target the WebDrake repo when submitting the pull, since otherwise it will go to Phobos master by default, which is probably not what you want). T -- Debugging is twice as hard as writing the code in the first place. Therefore, if you write the code as cleverly as possible, you are, by definition, not smart enough to debug it. -- Brian W. Kernighan
Nov 28 2013
prev sibling parent reply "Daniel Murphy" <yebblies nospamgmail.com> writes:
"Craig Dillabaugh" <craig.dillabaugh gmail.com> wrote in message 
news:izowcplookyzxrpzmoci forum.dlang.org...
 On Saturday, 23 November 2013 at 04:37:19 UTC, Shammah Chancellor wrote:
 //Hijack

 http://digitalmars.com/d/1.0/cppcomplex.html

 . Consider the formula (1 - infinity*i) * i which should produce 
 (infinity + i). However, if instead the second factor is (0 + i) rather 
 than just i, the result is (infinity + NaN*i), a spurious NaN was 
 generated.
 . A distinct imaginary type preserves the sign of 0, necessary for 
 calculations involving branch cuts.

Is this stuff no longer an issue? -Shammah

I believe D used to have builtin complex types, back in the old days. They have been removed (deprecated?) and replaced by the library type std.complex. At least that is my understanding. Craig

They are not deprecated yet, but it has been 'planned' for a long time. It's one of those things where not deprecating them doesn't hurt anyone, so it isn't high priority.
Nov 23 2013
parent reply "Daniel Murphy" <yebblies nospamgmail.com> writes:
"Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> wrote in message 
news:pyjudfnteduztpporndj forum.dlang.org...
 On Saturday, 23 November 2013 at 08:18:35 UTC, Daniel Murphy wrote:
 They are not deprecated yet, but it has been 'planned' for a long time.
 It's one of those things where not deprecating them doesn't hurt anyone, 
 so
 it isn't high priority.

Must say that, whatever the behind-the-scenes of the implementation, it seems a shame to lose the nice z = x + y*i notation. OTOH I guess that could lead to ambiguous code, e.g. int i = 4; complex z = x + y*i; // what does i mean here?

I feel mostly the same way - except - I've never actually used them in my own code, and when I've built test cases for them while fixing other bugs I've found some horrific bugs in dmd. This makes me think that _nobody_ is using them. And if that's true, we might as well get rid of them. eg https://d.puremagic.com/issues/show_bug.cgi?id=7594
Nov 23 2013
parent reply Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-23 14:02:38 +0000, Daniel Murphy said:

 "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> wrote in 
 message news:pyjudfnteduztpporndj forum.dlang.org...
 On Saturday, 23 November 2013 at 08:18:35 UTC, Daniel Murphy wrote:
 They are not deprecated yet, but it has been 'planned' for a long time.
 It's one of those things where not deprecating them doesn't hurt anyone, so
 it isn't high priority.

Must say that, whatever the behind-the-scenes of the implementation, it seems a shame to lose the nice z = x + y*i notation. OTOH I guess that could lead to ambiguous code, e.g. int i = 4; complex z = x + y*i; // what does i mean here?

I feel mostly the same way - except - I've never actually used them in my own code, and when I've built test cases for them while fixing other bugs I've found some horrific bugs in dmd. This makes me think that _nobody_ is using them. And if that's true, we might as well get rid of them. eg https://d.puremagic.com/issues/show_bug.cgi?id=7594

I disagree. I was using them for physics simulations. They are very useful for the computational physics community. Just because most people are still using FORTRAN does not mean they won't switch eventually. -Shammah
Nov 23 2013
parent Shammah Chancellor <anonymous coward.com> writes:
On 2013-11-24 15:50:46 +0000, Joseph Rushton Wakeling said:

 On Saturday, 23 November 2013 at 15:13:22 UTC, Shammah Chancellor wrote:
 I disagree.  I was using them for physics simulations.   They are very 
 useful for the computational physics community.   Just because most 
 people are still using FORTRAN does not mean they won't switch 
 eventually.

Would it cause you any particular disadvantage to use the library std.complex rather than the built-in complex type?

It would if the they don't work correctly. There needs to be an Imaginary type and some proper operations between complex and imaginary types. That doesn't seem to be the case currently. I personally think having the built-in type is very helpful. However, I can understand from a language perspective that having "i" around is hard for the parser. Also, the argument "If complex/imaginary is built-in, why not have quaterions also" seems to imply that it should be a library type. -Shammah
Nov 24 2013
prev sibling next sibling parent "Craig Dillabaugh" <craig.dillabaugh gmail.com> writes:
On Saturday, 23 November 2013 at 04:37:19 UTC, Shammah Chancellor 
wrote:
 //Hijack

 http://digitalmars.com/d/1.0/cppcomplex.html

 • Consider the formula (1 - infinity*i) * i which should 
 produce (infinity + i). However, if instead the second factor 
 is (0 + i) rather than just i, the result is (infinity + 
 NaN*i), a spurious NaN was generated.
 • A distinct imaginary type preserves the sign of 0, necessary 
 for calculations involving branch cuts.

Is this stuff no longer an issue? -Shammah

I believe D used to have builtin complex types, back in the old days. They have been removed (deprecated?) and replaced by the library type std.complex. At least that is my understanding. Craig
Nov 22 2013
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Saturday, 23 November 2013 at 08:18:35 UTC, Daniel Murphy 
wrote:
 They are not deprecated yet, but it has been 'planned' for a 
 long time.
 It's one of those things where not deprecating them doesn't 
 hurt anyone, so
 it isn't high priority.

Must say that, whatever the behind-the-scenes of the implementation, it seems a shame to lose the nice z = x + y*i notation. OTOH I guess that could lead to ambiguous code, e.g. int i = 4; complex z = x + y*i; // what does i mean here?
Nov 23 2013
prev sibling next sibling parent "bearophile" <bearophileHUGS lycos.com> writes:
Shammah Chancellor:
 //Hijack

More hijack, see also: http://forum.dlang.org/thread/lxuhrynxujrujoksrvdi forum.dlang.org Bye, bearophile
Nov 23 2013
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Saturday, 23 November 2013 at 15:13:22 UTC, Shammah Chancellor 
wrote:
 I disagree.  I was using them for physics simulations.   They 
 are very useful for the computational physics community.   Just 
 because most people are still using FORTRAN does not mean they 
 won't switch eventually.

Would it cause you any particular disadvantage to use the library std.complex rather than the built-in complex type?
Nov 24 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
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

It's because 0.0L * (-real.infinity) evaluates to nan.
Nov 24 2013
prev sibling next sibling parent "bearophile" <bearophileHUGS lycos.com> writes:
Shammah Chancellor:

 It's more a fundamental problem with a complex type in general.
   C++ has this issue as well.   You need a purely imaginary 
 type with the appropiate operations between Complex and 
 Imaginary defined.

Can't you add a new name to std.complex to implement the purely imaginary type? Bye, bearophile
Nov 24 2013
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
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.
Nov 24 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
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?
Nov 25 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
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.
Nov 25 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
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.
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 23/11/13 10:42, Joseph Rushton Wakeling wrote:
 Must say that, whatever the behind-the-scenes of the implementation, it seems a
 shame to lose the nice z = x + y*i notation.

 OTOH I guess that could lead to ambiguous code, e.g.

      int i = 4;
      complex z = x + y*i;     // what does i mean here?

I was wrong about this -- you have to write complex z = x + y * 1i ... which is unambiguous though I guess potentially still prone to typos :-)
Nov 26 2013
prev sibling next sibling parent Iain Buclaw <ibuclaw gdcproject.org> writes:
On 19 November 2013 02:03, Andrei Alexandrescu
<SeeWebsiteForEmail erdani.org> wrote:
 On 11/18/13 5:44 PM, Craig Dillabaugh wrote:
 The complex type in std.complex restricts the real/imaginary
 parts of the number to be float/double/real.  I am curious to
 know if there is a reason why integral types are not permitted. I
 checked the C++ equivalent and it does not have the same
 requirement.

 I mention this because some of my work is done with radar
 satellite images. All pixels in such an images are stored as
 complex numbers, but in all cases I am aware of they are stored a
 short int values. Most software that operates on the images uses
 complex<short> (most of it is C++).

 Is there any reason why complex numbers in D's standard lib must
 be of non-integral types?  I believe in C++ the type is optimized
 for floating point values, but allows other types.

The simple reason is we couldn't find appropriate applications. If you make a good argument, we'll include integral types as well. Submit an enhancement request on bugzilla including your example and let's take it from there.

Gaussian integers / Graph plotting? IMO I don't see a reason why not to support in the library if it requires no change on the compiler side. Regards Iain,
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 02:02, Andrei Alexandrescu wrote:
 If you got the blessing of some complex number expert, that would be great.

I'm in Italy, I could ask the pope ... :-) As far as I can see the behaviour of std.complex.Complex is consistent both with C++'s std::complex and with D's own internal cfloat/cdouble/creal types. The only case we don't cover with std.complex is that of multiplying a complex number with a purely imaginary one (the equivalent of multiplying, say, a creal and an ireal; multiplying a creal and a creal, even if the real part of the latter is 0, will give the same results as std.complex). We _could_ tweak Complex internally so that in the event that its real part is 0, it behaves like a purely imaginary number, and if its imaginary part is 0, it behaves like a purely real number; but that's problematic for reasons already discussed, because it doesn't distinguish between a purely imaginary number versus one where the real part is so vanishingly small that FP considers it to be 0. (It also leaves std.complex' behaviour incompatible with std::complex and other library implementations, though I guess that's less of a concern so long as we think what std.complex does is right.) So, as other people have suggested, really the only thing we can reasonably do is to define a separate Imaginary type -- or get round the problems with the built-in complex and imaginary types. It is really a shame if they are as problematic as I've heard. :-(
Nov 26 2013
prev sibling next sibling parent "John Colvin" <john.loughran.colvin gmail.com> writes:
On Tuesday, 26 November 2013 at 11:22:38 UTC, Joseph Rushton 
Wakeling wrote:
 On 26/11/13 02:02, Andrei Alexandrescu wrote:
 If you got the blessing of some complex number expert, that 
 would be great.

I'm in Italy, I could ask the pope ... :-) As far as I can see the behaviour of std.complex.Complex is consistent both with C++'s std::complex and with D's own internal cfloat/cdouble/creal types. The only case we don't cover with std.complex is that of multiplying a complex number with a purely imaginary one (the equivalent of multiplying, say, a creal and an ireal; multiplying a creal and a creal, even if the real part of the latter is 0, will give the same results as std.complex). We _could_ tweak Complex internally so that in the event that its real part is 0, it behaves like a purely imaginary number, and if its imaginary part is 0, it behaves like a purely real number; but that's problematic for reasons already discussed, because it doesn't distinguish between a purely imaginary number versus one where the real part is so vanishingly small that FP considers it to be 0. (It also leaves std.complex' behaviour incompatible with std::complex and other library implementations, though I guess that's less of a concern so long as we think what std.complex does is right.) So, as other people have suggested, really the only thing we can reasonably do is to define a separate Imaginary type -- or get round the problems with the built-in complex and imaginary types. It is really a shame if they are as problematic as I've heard. :-(

It seems to me that one really shouldn't special case for absolute values when dealing with floating point. A floating Imaginary type and an integer Complex type - both in std.complex - are the correct solutions IMO.
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 13:52, John Colvin wrote:
 It seems to me that one really shouldn't special case for absolute values when
 dealing with floating point. A floating Imaginary type and an integer Complex
 type - both in std.complex - are the correct solutions IMO.

Yup, agree.
Nov 26 2013
prev sibling next sibling parent "John Colvin" <john.loughran.colvin gmail.com> writes:
On Tuesday, 26 November 2013 at 13:20:39 UTC, Joseph Rushton 
Wakeling wrote:
 On 26/11/13 13:52, John Colvin wrote:
 It seems to me that one really shouldn't special case for 
 absolute values when
 dealing with floating point. A floating Imaginary type and an 
 integer Complex
 type - both in std.complex - are the correct solutions IMO.

Yup, agree.

Just had a chat with one of our local numerics experts: He agreed that hidden special casing that breaks standards compliance is bad and that whatever the solution was it should be explicit. However, the following situation arises: Complex!double(1, -double.inf) * Complex!int(0, i); What does IEEE say about that? 0 * inf == nan?
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 14:32, John Colvin wrote:
 What does IEEE say about that? 0 * inf == nan?

Can't speak for the standard, but in both D and C++, 0 * inf evaluates to nan.
Nov 26 2013
prev sibling next sibling parent "David Nadlinger" <code klickverbot.at> writes:
On Monday, 25 November 2013 at 08:18:43 UTC, Joseph Rushton 
Wakeling wrote:
 But if you just think of 0 as a number, then

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

This is where your argument falls apart, as mathematically, you can't do that unless lim{x --> inf} x is well-defined. See also: http://en.wikipedia.org/wiki/Riemann_sphere David
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 17:11, David Nadlinger wrote:
 On Monday, 25 November 2013 at 08:18:43 UTC, Joseph Rushton Wakeling wrote:
 But if you just think of 0 as a number, then

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

This is where your argument falls apart, as mathematically, you can't do that unless lim{x --> inf} x is well-defined.

I was using a very lazy shorthand there, I'm glad someone thought to call me on it. Can we take it as read that I was basically thinking of a sequence {x_n} such that for every K there is an N such that for n > N, x_n > K ... ? :-) And then you have 0 * lim{n --> inf} x_n = ... etc. The fun stuff must surely arrive when you want to show this kind of stuff in the context of real numbers being defined as equivalence classes of infinite sequences of rationals, à la Cauchy ...
 See also: http://en.wikipedia.org/wiki/Riemann_sphere

I don't recall ever actually studying the Riemann sphere, which really seems to me like a gap in my education :-\
Nov 26 2013
prev sibling next sibling parent "David Nadlinger" <code klickverbot.at> writes:
On Tuesday, 26 November 2013 at 16:30:30 UTC, Joseph Rushton 
Wakeling wrote:
 On 26/11/13 17:11, David Nadlinger wrote:
 On Monday, 25 November 2013 at 08:18:43 UTC, Joseph Rushton 
 Wakeling wrote:
 But if you just think of 0 as a number, then

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

This is where your argument falls apart, as mathematically, you can't do that unless lim{x --> inf} x is well-defined.

I was using a very lazy shorthand there, I'm glad someone thought to call me on it. Can we take it as read that I was basically thinking of a sequence {x_n} such that for every K there is an N such that for n > N, x_n > K ... ? :-) And then you have 0 * lim{n --> inf} x_n = ... etc.

x_n = n actually fulfils that property, and I think most people would understand the limit notation for real numbers exactly the way you intended. But that was not my point. To be able to write "lim{n --> inf} x_n" in a meaningful way (and consequently »pull in« the multiplication), a (or as it turns out, the) limit must exist in the metric space you are working in. If your metric space contains ∞, and it has the property that 0 . ∞ = 0, then your argument is correct. But such a symbol ∞ does not exist in the real numbers. I guess it might help to think back to your first university-level analysis courses, where I'm sure these distinctions were discussed many times in proofs concerning the existence of limits, e.g. integrability of certain functions, …
 See also: http://en.wikipedia.org/wiki/Riemann_sphere

I don't recall ever actually studying the Riemann sphere, which really seems to me like a gap in my education :-\

Well, the reason I bring this up is that what the »right« behavior is all comes down to the definition of your numerical system. IEEE 754 includes infinity as an actual value, contrary to the usual definition of real numbers in mathematics. However, it also distinguishes between +∞ and -∞, so it can't model the Riemann sphere, which is one of the most straightforward ways to perform the extension of the complex plane with a concept infinity in mathematics. David
Nov 26 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 22:11, David Nadlinger wrote:
 x_n = n actually fulfils that property, and I think most people would
understand
 the limit notation for real numbers exactly the way you intended.

 But that was not my point. To be able to write "lim{n --> inf} x_n" in a
 meaningful way (and consequently »pull in« the multiplication), a (or as it
 turns out, the) limit must exist in the metric space you are working in.

 If your metric space contains ∞, and it has the property that 0 . ∞ = 0,
then
 your argument is correct. But such a symbol ∞ does not exist in the real
numbers.

 I guess it might help to think back to your first university-level analysis
 courses, where I'm sure these distinctions were discussed many times in proofs
 concerning the existence of limits, e.g. integrability of certain functions,
…

Well, it has been 12+ years ... :-P Still, it's very, very annoying when one misplaces fundamental stuff like that. Thank you for reminding me :-) Now I need to dust off my copy of "What is Mathematics?" ...
 Well, the reason I bring this up is that what the »right« behavior is all
comes
 down to the definition of your numerical system.

 IEEE 754 includes infinity as an actual value, contrary to the usual definition
 of real numbers in mathematics. However, it also distinguishes between +∞ and
 -∞, so it can't model the Riemann sphere, which is one of the most
 straightforward ways to perform the extension of the complex plane with a
 concept infinity in mathematics.

I'm sure you've heard that old anecdote of the professor back in the 1950s, or was it the 1920s, who, on hearing a student say "infinity", said: "I won't have bad language in class!" :-)
Nov 27 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 17:28, Andrei Alexandrescu wrote:
 On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

I've started work, code is here: https://github.com/WebDrake/phobos/tree/imaginary Feedback, forks/pull requests, etc. welcome.
Nov 27 2013
prev sibling next sibling parent Iain Buclaw <ibuclaw gdcproject.org> writes:
On 27 November 2013 22:29, Joseph Rushton Wakeling
<joseph.wakeling webdrake.net> wrote:
 On 26/11/13 17:28, Andrei Alexandrescu wrote:
 On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

I've started work, code is here: https://github.com/WebDrake/phobos/tree/imaginary Feedback, forks/pull requests, etc. welcome.

That repo doesn't seem to exist (must by my imagination).
Nov 28 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 29/11/13 08:50, Iain Buclaw wrote:
 That repo doesn't seem to exist (must by my imagination).

It's just the branch "imaginary" in my regular Phobos repo (I always use feature branches for new stuff). The link works for me, what goes wrong for you?
Nov 29 2013
prev sibling next sibling parent Iain Buclaw <ibuclaw gdcproject.org> writes:
On 29 November 2013 08:57, Joseph Rushton Wakeling
<joseph.wakeling webdrake.net> wrote:
 On 29/11/13 08:50, Iain Buclaw wrote:
 That repo doesn't seem to exist (must by my imagination).

It's just the branch "imaginary" in my regular Phobos repo (I always use feature branches for new stuff). The link works for me, what goes wrong for you?

Turned out the RJ45 cable in my computer was imaginary. :o)
Nov 29 2013
prev sibling next sibling parent Iain Buclaw <ibuclaw gdcproject.org> writes:
On 29 November 2013 11:50, Iain Buclaw <ibuclaw gdcproject.org> wrote:
 On 29 November 2013 08:57, Joseph Rushton Wakeling
 <joseph.wakeling webdrake.net> wrote:
 On 29/11/13 08:50, Iain Buclaw wrote:
 That repo doesn't seem to exist (must by my imagination).

It's just the branch "imaginary" in my regular Phobos repo (I always use feature branches for new stuff). The link works for me, what goes wrong for you?

Turned out the RJ45 cable in my computer was imaginary. :o)

Fixed by switching to real. :-P
Nov 29 2013
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 26/11/13 17:28, Andrei Alexandrescu wrote:
 On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

Provisional code: https://github.com/D-Programming-Language/phobos/pull/1797 I also realized there was no explicit issue report, so created one: https://d.puremagic.com/issues/show_bug.cgi?id=11787
Dec 20 2013
prev sibling next sibling parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Tuesday, 26 November 2013 at 21:11:25 UTC, David Nadlinger 
wrote:
 IEEE 754 includes infinity as an actual value, contrary to the 
 usual definition of real numbers in mathematics.

A bit late, but AFAIK inf represents either overflow or N/0... And 0 represents either underflow or zero. Computations on those values should be conservative unless you trap overflow/underflow exceptions and handle those as special cases. You want to preserve overflow... Then you have the denormal numbers (underflow where you retain some digits). It is tempting to think of FP as real numbers, but they are not, of course, so libraries should IMO be conservative.
Jan 01 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Wednesday, 1 January 2014 at 12:29:35 UTC, Stewart Gordon 
wrote:
 I don't understand.  At the moment Complex appears to me to be 
 type-agnostic - as long as a type supports the standard 
 arithmetic operators and assignment of the value 0 to it, it 
 will work.  The only thing preventing it from working at the 
 moment is this line

     struct Complex(T)  if (isFloatingPoint!T)

 So why do you need an appropriate application in order not to 
 have this arbitrary restriction?  Or have I missed something?

There are binary operations on complex numbers where the only sensible outcome seems to be non-integral real and imaginary parts. Addition, subtraction and multiplication are OK with integral types, but division really seems unpleasant to implement absent floating point, exponentiation even more so. I imagine there are ways to resolve this, but it certainly simplifies implementation to assume floating-point, and absent a compelling application there is not much reason to avoid that simplification.
 It isn't just integer types.  Somebody might want to use 
 complex with a library (fixed-point, arbitrary precision, 
 decimal, etc.) numeric type.
  Fractal generators, for example, are likely to use this a lot.

I agree that such any numeric type that effectively models a real number should be supported. In principle it ought to be sufficient to check that the required "floating-point-ish" operations (including sin and cos) are supported, plus maybe some tweaks to how internal temporary values are handled. However, I think relaxing the template constraints like this would best be done in the context of a library float-esque type (e.g. BigFloat) being implemented in Phobos, which could then be used to provide both proof-of-concept and the primary test case.
 Or even more exotically, use Complex!(Complex!real) to 
 implement hypercomplex numbers.

Perhaps best implemented as a type in its own right? :-)
Jan 01 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Tuesday, 19 November 2013 at 01:44:32 UTC, Craig Dillabaugh 
wrote:
 The complex type in std.complex restricts the real/imaginary
 parts of the number to be float/double/real.  I am curious to
 know if there is a reason why integral types are not permitted. 
 I
 checked the C++ equivalent and it does not have the same
 requirement.

Quoting the C++ standard, §26.4: "The effect of instantiating the template complex for any type other than float, double, or long double is unspecified." So even if some implementations support it, it is *not* standard C++. Lars
Jan 01 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Wednesday, 1 January 2014 at 19:55:58 UTC, Joseph Rushton 
Wakeling wrote:
 On Wednesday, 1 January 2014 at 12:29:35 UTC, Stewart Gordon 
 wrote:
 [...] why do you need an appropriate application in order not 
 to have this arbitrary restriction?  Or have I missed 
 something?

There are binary operations on complex numbers where the only sensible outcome seems to be non-integral real and imaginary parts. Addition, subtraction and multiplication are OK with integral types, but division really seems unpleasant to implement absent floating point, exponentiation even more so. I imagine there are ways to resolve this, but it certainly simplifies implementation to assume floating-point, and absent a compelling application there is not much reason to avoid that simplification.

I agree completely. This is the main reason why I chose to explicitly disallow integral types when I wrote std.complex.
 It isn't just integer types.  Somebody might want to use 
 complex with a library (fixed-point, arbitrary precision, 
 decimal, etc.) numeric type.
 Fractal generators, for example, are likely to use this a lot.

I agree that such any numeric type that effectively models a real number should be supported. In principle it ought to be sufficient to check that the required "floating-point-ish" operations (including sin and cos) are supported, plus maybe some tweaks to how internal temporary values are handled.

Agreed. Any "floating point-like" library type which is to be supported by std.complex must either have cos(), sin(), hypot(), etc. as member functions (since std.complex cannot import non-std modules), or the type needs to be supported by std.math as well. If so, std.math is the place to start, not std.complex.
Jan 01 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Wednesday, 1 January 2014 at 12:29:35 UTC, Stewart Gordon 
wrote:
  [...]

 Or even more exotically, use Complex!(Complex!real) to 
 implement hypercomplex numbers.

This is an extremely marginal use case. Currently Complex!(Complex!T) folds to Complex!T. I thought this was specified in the module documentation, but if it was, it's been removed. It is explained in a comment in the source code, though: /* Makes Complex!(Complex!T) fold to Complex!T. The rationale for this is that just like the real line is a subspace of the complex plane, the complex plane is a subspace of itself. Example of usage: --- Complex!T addI(T)(T x) { return x + Complex!T(0.0, 1.0); } --- The above will work if T is both real and complex. */ template Complex(T) if (is(T R == Complex!R)) { alias T Complex; }
Jan 01 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Friday, 20 December 2013 at 16:26:00 UTC, Joseph Rushton 
Wakeling wrote:
 On 26/11/13 17:28, Andrei Alexandrescu wrote:
 On 11/26/13 3:22 AM, Joseph Rushton Wakeling wrote:
 So, as other people have suggested, really the only thing we 
 can
 reasonably do is to define a separate Imaginary type

I agree. Let's move forward with this.

Provisional code: https://github.com/D-Programming-Language/phobos/pull/1797 I also realized there was no explicit issue report, so created one: https://d.puremagic.com/issues/show_bug.cgi?id=11787

I'm not 100% convinced that a pure imaginary type is the way to go. You may be interested in reading a previous discussion, where the subject of pure imaginary number support is brought up: http://forum.dlang.org/thread/flsf9u$jf$1 digitalmars.com Let me quote some of the replies in that thread: On Monday, 7 January 2008 at 08:05:48 UTC, Don Clugston wrote:
 I think the argument for pure imaginary types is extremely 
 weak. They are very annoying to work with, because they aren't 
 closed under multiplication. This is a nasty property to have 
 in a built-in type.

 Suppose you're writing a generic function product(T)(T[]) which 
 returns T[0]*T[1]*T[2]*... What's the return type?
 Suppose T is idouble. Then if the number of elements in T is 
 odd, the return type should be idouble, but if it's even, the 
 return type should be double!
 You could promote it to complex with a construction like 
 typeof(1?T*T:T), but that's inefficient, and negates most of 
 the benefits of having a imaginary type.

On Tuesday, 8 January 2008 at 03:13:34 UTC, Bill Baxter wrote:
 The bottom line is that there really is no useful mathematics 
 that can be performed using only imaginary numbers.  If to do 
 anything more than add, subtract, or multiply by a (real!) 
 scalar they become complex.  If you have any imaginary numbers 
 in a numerical code it means you're working in the complex 
 plane.

Jan 01 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Wednesday, 1 January 2014 at 23:54:05 UTC, Walter Bright wrote:
 On 11/26/2013 3:22 AM, Joseph Rushton Wakeling wrote:
 (It also leaves std.complex' behaviour incompatible with 
 std::complex and other
 library implementations, though I guess that's less of a 
 concern so long as we
 think what std.complex does is right.)

Making it behave differently than others would be a disaster.

Just as well that I was speaking hypothetically about things that _could_ be done, rather than _should_ be ... :-)
Jan 01 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Wednesday, 1 January 2014 at 23:32:58 UTC, Lars T. Kyllingstad 
wrote:
 I'm not 100% convinced that a pure imaginary type is the way to 
 go.  You may be interested in reading a previous discussion, 
 where the subject of pure imaginary number support is brought 
 up:

   http://forum.dlang.org/thread/flsf9u$jf$1 digitalmars.com

I have read through that thread already, but thank you for bringing it up -- it raises some issues that deserve an answer. My own thoughts on the matter go something like this: * There are calculations involving "pair of reals" implementations of complex numbers that generate spuriously incorrect results. They're corner cases, but they exist. * A purely imaginary type allows those calculations to be performed correctly. It also allows more precise calculations in general for cases where the real part of a complex number is known to be 0. * The idea of allowing imaginary literals but promoting them to complex when written to a variable sounds attractive but we can't discount the possibility that people will want to pass imaginary literals to functions or otherwise store them in variables, and later use them; if they are promoted to complex, calculations done with them may have lesser precision or even generate the aforementioned spurious results. I don't think simply writing to a variable should cause loss of precision like this. * You can do calculations involving purely real-valued numbers and complex numbers and not run into the same issues, because purely real values are supported. So you should be able to do the same with purely imaginary numbers. * The lack of closure under various operations is annoying but not insurmountable. I should add that as far as I'm concerned what I want is simply to find the best possible way to represent and use purely imaginary numbers in D. I don't mind if, having done that, the code winds up being rejected, as long as the exercise proves useful. Of course, an alternative would be to tweak the internals of Complex so that it can represent "purely real" and "purely imaginary" types precisely. I'm writing from a phone now so it's not really convenient to explain at length, but I'll try to expand on the above in the next days.
Jan 01 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Thursday, 2 January 2014 at 02:59:05 UTC, Andrei Alexandrescu
wrote:
 On 1/1/14 3:12 PM, Lars T. Kyllingstad wrote:
 Currently Complex!(Complex!T)
 folds to Complex!T.

What was the motivation?

To simplify generic code, where you have a function that takes both real and complex arguments and returns a complex number, e.g.: Complex!T fun(T)(T x) { return x * someComplexNumber; } In the above, T may be real or complex.
Jan 02 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Tuesday, 26 November 2013 at 12:52:30 UTC, John Colvin wrote:
 It seems to me that one really shouldn't special case for 
 absolute values when dealing with floating point.

Why not? std.math and std.mathspecial do it all the time.
Jan 02 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Tuesday, 26 November 2013 at 13:32:22 UTC, John Colvin wrote:
 Just had a chat with one of our local numerics experts: He 
 agreed that hidden special casing that breaks standards 
 compliance is bad and that whatever the solution was it should 
 be explicit.

I don't think we should worry too much about standards compliance. A library Complex type is quite different from a hardware floating-point type.
Jan 02 2014
prev sibling next sibling parent "Craig Dillabaugh" <cdillaba cg.scs.carleton.ca> writes:
On Wednesday, 1 January 2014 at 12:29:35 UTC, Stewart Gordon
wrote:
 On 19/11/2013 02:03, Andrei Alexandrescu wrote:> On 11/18/13 
 5:44 PM, Craig Dillabaugh wrote:
 <snip>
 Is there any reason why complex numbers in D's standard lib 
 must
 be of non-integral types?  I believe in C++ the type is 
 optimized
 for floating point values, but allows other types.

The simple reason is we couldn't find appropriate applications.

I don't understand. At the moment Complex appears to me to be type-agnostic - as long as a type supports the standard arithmetic operators and assignment of the value 0 to it, it will work. The only thing preventing it from working at the moment is this line struct Complex(T) if (isFloatingPoint!T) So why do you need an appropriate application in order not to have this arbitrary restriction? Or have I missed something? It isn't just integer types. Somebody might want to use complex with a library (fixed-point, arbitrary precision, decimal, etc.) numeric type. Fractal generators, for example, are likely to use this a lot. Or even more exotically, use Complex!(Complex!real) to implement hypercomplex numbers. Stewart.

Since I originally proposed this I should chime in again. My use-case for non floating point, complex numbers is radar image processing, where the radar intensity/phase are stored as complex numbers which are quantized as 16-bit integer values. Andrei suggested I come up with a proposal for non-FP complex, but reading this thread, and my experience working with integer valued complex values in C++ has now reversed my opinion. I think the current D requirement is good. Just a small example, the norm() method calculates the squared magnitude of the complex number (eg. you get a single real value). This is an operation that we use extensively in our work. This created problems with generating lots of NaN values due to integer overflows, that were somewhat tricky to find. I can imagine numerous issues like this popping up if using integral values with complex numbers.
Jan 02 2014
prev sibling next sibling parent =?UTF-8?B?U2ltZW4gS2rDpnLDpXM=?= <simen.kjaras gmail.com> writes:
On 2014-01-01 12:29, Stewart Gordon wrote:
 Or even more exotically, use Complex!(Complex!real) to implement
 hypercomplex numbers.

For a more generic solution to this, see Cayley-Dickson construction[1] and my implementation of such: https://github.com/Biotronic/Collectanea/blob/master/biotronic/CayleyDickson2.d It's nowhere near finished, has some bugs, and I haven't worked on it for at least a year, but it's an interesting way to create complex, hypercomplex, dual, and split-complex numbers (and combinations thereof). [1]: http://en.wikipedia.org/wiki/Cayley-Dickson_construction -- Simen
Jan 02 2014
prev sibling next sibling parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Thursday, 2 January 2014 at 11:37:22 UTC, Lars T. Kyllingstad 
wrote:
 I don't think we should worry too much about standards 
 compliance. A library Complex type is quite different from a 
 hardware floating-point type.

Are you sure? Sometimes you need to translate an algorithm, you don't understand the inner workings of, from a codebase/cookbook. If std.complex differs from the most used c++/fortran implementations people will be confused, and you also end up having (machine translated) algorithm libraries each supplying their own complex type. Use 3 different libraries and you have to deal with 3 different complex types. Floating point is rather sensitive to reordering of instructions, so I'd say you'll be better off mirroring one of the major existing implementations, otherwise accumulated discrepancies will be blamed on the language... A new tool that produce the same results as the old proven dinosaur tool look trustworthy. It makes you think that the conversion of your algorithms to the new tool was a success.
Jan 02 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Thursday, 2 January 2014 at 18:12:56 UTC, Stewart Gordon wrote:
 Then why not just disable division if it's a non-float type, 
 rather than preventing the whole complex template from being 
 used with that type? This is like cutting off somebody's arm 
 because they have a sore thumb.

Because that also seems to me to be an unpleasant option. A complex number implementation that fails to support ordinary arithmetic operations in all circumstances is pretty non-intuitive and will probably lead to unpleasant bugs in users' code.
 Moreover, we have no way in the general case of determining 
 whether T is an integral type, a library float-esque type, or 
 (for example) a Galois field type.  So disabling it _just in 
 case_ division doesn't work is crazy.  There must be better 
 ways to do it.

I don't follow your point here. We can constrain T however we see fit. The point isn't to have some perfect representation of every mathematical possibility, it's to have useful code that serves a good range of use-cases while being robust and easy to maintain. Restricting T to floating point types is a useful simplification here that has few costs in terms of expected use-cases.
 Exponentiation by a non-negative integer is straightforward.  
 So we should at least support this case for Gaussian integers.

Again, I don't see it as being useful to have a Complex!T whose support for binary operations may vary depending on the type of T.
 What do you mean by "in the context of", exactly?  Restricting 
 it to some float-esque type that is in Phobos would still be 
 overly restrictive.

No, I mean that until you have at least one actual float-esque type to test with, it is probably unwise to relax the template constraints that currently mandate built-in FP types.
 It would, but removing the restriction would simplify the 
 implementation of Hypercomplex(T) by enabling it to be a 
 wrapper for Complex!(Complex!T).

And complicate the implementation of Complex itself, for the sake of a likely very marginal special interest that could be supported quite well by an independent type.
Jan 02 2014
prev sibling next sibling parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Thursday, 2 January 2014 at 19:32:32 UTC, Joseph Rushton 
Wakeling wrote:
 No, I mean that until you have at least one actual float-esque 
 type to test with, it is probably unwise to relax the template 
 constraints that currently mandate built-in FP types.

Fixed-point support might not be important on x86, but it has been needed and used for complex numbers in the past on other CPUs.
Jan 02 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Thursday, 2 January 2014 at 19:50:46 UTC, Ola Fosheim Grøstad 
wrote:
 Fixed-point support might not be important on x86, but it has 
 been needed and used for complex numbers in the past on other 
 CPUs.

I agree that should be supported. I just think that template constraints should only be relaxed in the context of concrete test cases, otherwise you're making promises you don't know you can keep.
Jan 02 2014
prev sibling next sibling parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Thursday, 2 January 2014 at 19:58:27 UTC, Joseph Rushton 
Wakeling wrote:
 I agree that should be supported. I just think that template 
 constraints should only be relaxed in the context of concrete 
 test cases, otherwise you're making promises you don't know you 
 can keep.

Yes, that is very true. Fixed-point is very sensitive to precision-loss so it requires some careful thinking, so it is definitively not a drop-in replacement type. (It is actually useful on x86 too when you want reproducible results across binaries and hardware, such as in peer-to-peer syncing of computed state.)
Jan 02 2014
prev sibling next sibling parent "QAston" <qaston gmail.com> writes:
On Sunday, 24 November 2013 at 17:35:34 UTC, Shammah Chancellor 
wrote:
 On 2013-11-24 15:50:46 +0000, Joseph Rushton Wakeling said:

 On Saturday, 23 November 2013 at 15:13:22 UTC, Shammah 
 Chancellor wrote:
 I disagree.  I was using them for physics simulations.   They 
 are very useful for the computational physics community.   
 Just because most people are still using FORTRAN does not 
 mean they won't switch eventually.

Would it cause you any particular disadvantage to use the library std.complex rather than the built-in complex type?

It would if the they don't work correctly. There needs to be an Imaginary type and some proper operations between complex and imaginary types. That doesn't seem to be the case currently. I personally think having the built-in type is very helpful. However, I can understand from a language perspective that having "i" around is hard for the parser. Also, the argument "If complex/imaginary is built-in, why not have quaterions also" seems to imply that it should be a library type. -Shammah

You can have im!5.0 in a library type, to me it sounds good enough.
Jan 02 2014
prev sibling next sibling parent "Joseph Rushton Wakeling" <joseph.wakeling webdrake.net> writes:
On Thursday, 2 January 2014 at 20:23:40 UTC, Stewart Gordon wrote:
 The compiler rejecting the code is the most pleasant bug that's 
 possible IMO.

It's still unpleasant relative to the intuitive expectation that numerical types should support all basic numerical operations.
 Not being overly restrictive in what types you will allow it to 
 be used with is an important part of serving a good range of 
 use cases.

It's a question of balance: breadth of support vs. ease of implementation and maintenance. Some use cases may be better served in other ways than rolling them into one "catch-all" type.
 I'm sure we have a small handful already.  We just need to find 
 them. For instance, given the time I could probably dig up my 
 rational number implementation and update it to current D.

I'd really like to see that anyway, for its own sake; I've been working on a std.rational based on David Simcha's work, it would be good to compare.
 How would it complicate the implementation?  Removing the 
 undocumented rule whereby Complex!(Complex!T) folds to 
 Complex!T would be a slight simplification.

It would involve non-trivial revisions to the internals of Complex. If you want to submit a pull request supporting this I'm sure it'll be considered carefully, but it is not something to be done casually.
 Maybe the right course of action is to have a parameter to the 
 template that suppresses the type restriction and the folding 
 rule, so that people who want to use it on a type that might 
 not work properly can do so.  This would be a relatively small 
 change.

I do not know how familiar you are with the internals of std.complex.Complex but I think it would be a good idea to go through them in some detail before asserting that changes like this would be "relatively small".
Jan 02 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Thursday, 2 January 2014 at 01:17:54 UTC, Joseph Rushton 
Wakeling wrote:
    * There are calculations involving "pair of reals" 
 implementations of complex
      numbers that generate spuriously incorrect results. 
 They're corner cases,
      but they exist.

    * A purely imaginary type allows those calculations to be 
 performed
      correctly. It also allows more precise calculations in 
 general for cases
      where the real part of a complex number is known to be 0.

So we have three choices: 1. We add a dedicated Imaginary type, and leave Complex more or less the way it is. The pros are that it provides a way to represent imaginary numbers "correctly", while keeping Complex simple and performant. The cons are that Complex still gives wrong (or at least somewhat unexpected) results in the aforementioned corner cases, and that we add a new type with extremely limited usability. 2. Special-case Complex for imaginary numbers. The pros are that it solves the problems Imaginary was intended to solve, and we don't need a new type. The cons are that the Complex implementation becomes more complex (haha) and less performant, and that we actually change the semantics of IEEE floating-point "zero" somewhat. 3. Leave std.complex as-is, and make sure people know what the problematic cases are. They all suck. I don't know what is the lesser of three evils here, but I too am starting to lean towards 1. I'm probably going to continue playing devil's advocate, though. ;)
    * You can do calculations involving purely real-valued 
 numbers and complex
      numbers and not run into the same issues, because purely 
 real values are
      supported. So you should be able to do the same with 
 purely imaginary
      numbers.

That argument is fallacious. Imaginary numbers are quite different from real numbers.
 I should add that as far as I'm concerned what I want is simply 
 to find the best possible way to represent and use purely 
 imaginary numbers in D. I don't mind if, having done that, the 
 code winds up being rejected, as long as the exercise proves 
 useful.

I think it's great that you're doing this. :)
Jan 02 2014
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 02/01/14 22:08, QAston wrote:
 You can have im!5.0 in a library type, to me it sounds good enough.

As currently written, it'll be imaginary(5.0) or Imaginary!double(5.0) -- sorry for the verbosity, but "im" seems to me too likely to wind up being used as a variable name ... ;-)
Jan 02 2014
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 02/01/14 23:26, Lars T. Kyllingstad wrote:
 2. Special-case Complex for imaginary numbers. The pros are that it solves the
 problems Imaginary was intended to solve, and we don't need a new type. The
cons
 are that the Complex implementation becomes more complex (haha) and less
 performant, and that we actually change the semantics of IEEE floating-point
 "zero" somewhat.

Well, what I was mulling over was a Complex type that consists of two FP values (re and im) and two bools (let's call them hasRe, hasIm). The point of this is that if you assign to Complex from a regular real-valued number, then the value of re gets set, im gets set to zero, hasRe to true, and hasIm to false. So, you have a flag inside that says, "Hey, only worry about your real part." Then you could define an imaginary() helper function that creates a Complex number with re = 0.0, im set to whatever it needs to be, hasRe to false, and hasIm to true. So again, you can benefit from a Complex that _knows_ it's purely imaginary. Moreover, it should be possible to preserve that kind of knowledge under certain operations. e.g. if you multiply two Complex numbers which both have hasRe == false, you can _know_ that you're going to wind up with a number where hasRe == true and hasIm == false. OTOH when you perform certain other operations, you might end up with a zero re or im, but also with hasRe or hasIm to true, which basically reflects your _certainty_ that the real or imaginary part is zero. It might be interesting to try out, and it seems nicer to me than if ()'s based on whether re == 0 or im == 0, but on the other hand it enlarges the size of a Complex type and will result in a performance hit, and will break compatibility with C/C++ Complex types.
 They all suck. I don't know what is the lesser of three evils here, but I too
am
 starting to lean towards 1.  I'm probably going to continue playing devil's
 advocate, though. ;)

Please do. It's fun, and it also helps to improve things. :-)
    * You can do calculations involving purely real-valued numbers and complex
      numbers and not run into the same issues, because purely real values are
      supported. So you should be able to do the same with purely imaginary
      numbers.

That argument is fallacious. Imaginary numbers are quite different from real numbers.

Can you expand on that?
 I think it's great that you're doing this. :)

Thank you! :-)
Jan 02 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Thursday, 2 January 2014 at 18:37:36 UTC, Ola Fosheim Grøstad 
wrote:
 On Thursday, 2 January 2014 at 11:37:22 UTC, Lars T. 
 Kyllingstad wrote:
 I don't think we should worry too much about standards 
 compliance. A library Complex type is quite different from a 
 hardware floating-point type.

Are you sure?

Not at all. ;) I just think we should keep in mind why FP semantics are defined the way they are. Take 0.0*inf, for example. As has been mentioned, 0.0 may represent a positive real number arbitrarily close to zero, and inf may represent an arbitrarily large real number. The product of these is ill-defined, and hence represented by a NaN. 0.0+1.0i, on the other hand, represents a number which is arbitrarily close to i. Multiplying it with a very large real number gives you a number which has a very large imaginary part, but which is arbitrarily close to the imaginary axis, i.e. 0.0 + inf i. I think this is very consistent with FP semantics, and may be worth making a special case in std.complex.Complex.
 Sometimes you need to translate an algorithm, you don't 
 understand the inner workings of, from a codebase/cookbook. If 
 std.complex differs from the most used c++/fortran 
 implementations people will be confused, and you also end up 
 having (machine translated) algorithm libraries each supplying 
 their own complex type. Use 3 different libraries and you have 
 to deal with 3 different complex types.

I agree, but there is also a lot to be said for not repeating old mistakes, if we deem them to be such.
Jan 02 2014
prev sibling next sibling parent "Lars T. Kyllingstad" <public kyllingen.net> writes:
On Thursday, 2 January 2014 at 23:26:48 UTC, Joseph Rushton 
Wakeling wrote:
 On 02/01/14 23:26, Lars T. Kyllingstad wrote:
   * You can do calculations involving purely real-valued 
 numbers and complex
     numbers and not run into the same issues, because purely 
 real values are
     supported. So you should be able to do the same with 
 purely imaginary
     numbers.

That argument is fallacious. Imaginary numbers are quite different from real numbers.

Can you expand on that?

Mathematically, the real numbers are a ring, whereas the purely imaginary numbers are not. (I've been out of the abstract algebra game for a couple of years now, so please arrest me if I've remembered the terminology wrongly.) What it boils down to is that they are not closed under multiplication, which gives them radically different properties - or lack thereof.
Jan 02 2014
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 03/01/14 00:33, Stewart Gordon wrote:
 Please be specific.

You know, I'm starting to find your tone irritating. You are the one who's asking for functionality that goes beyond any Complex implementation that I'm aware of in any other language, and claiming that these things would be trivial to implement. I would expect a person who claims with confidence that something is trivial, to actually know the internals of the code well enough to understand what parts of it would need to be modified. On the basis of what you've written, I have no reason to believe that you do. There are numerous places inside current std.complex.Complex where temporary values are used mid-calculation. Those are all of type FPTemporary (which in practice means real). So, to handle library types (whether library floating-point types such as a BigFloat implementation, or a Complex!T so as to support hypercomplex numbers) you'd either have to special-case those functions or you'd have to provide an alternative Temporary template to handle the temporary internal values in the different cases. You'd also need to address questions of closure under operations (already an issue for the Imaginary type), and precision-related issues -- see e.g. the remarks by Craig Dillabaugh and Ola Fosheim Grøstad elsewhere in this thread. While it could be done, I don't think it would be simple to get right and I would regard it as a major revision of std.complex. I'd also be concerned that generalizing Complex in this way might lead to performance regressions for the standard Complex!T use-cases, due to the more complicated templates involved.
 Oh yes, some of the templated functions that take or return a complex would
need
 this extra parameter adding.

No. The internals would need significant rewriting, for reasons I've already elaborated.
 Another way it could be done is to have a Complex template that implements
these
 rules and which programmers would normally use, and have this aliasing
 ComplexImpl which actually provides the implementation and which programmers
can
 use directly if they want to bypass the restrictions.  The difficulty is
 documenting it in a way that would make sense to normal users....

If you want these things done in these ways, then I think the onus is on you to provide code which works and doesn't damage existing std.complex functionality. I don't personally see any rationale for implementing hypercomplex numbers as a specialization of Complex given that they can be just as well implemented as a type in their own right.
Jan 02 2014
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 03/01/14 01:04, Lars T. Kyllingstad wrote:
 Mathematically, the real numbers are a ring, whereas the purely imaginary
 numbers are not. (I've been out of the abstract algebra game for a couple of
 years now, so please arrest me if I've remembered the terminology wrongly.)
 What it boils down to is that they are not closed under multiplication, which
 gives them radically different properties - or lack thereof.

I've been out of the abstract algebra game for rather longer :-) but regardless of terminology, I understand what you mean. My point was meant to be somewhat simpler: if you have x * (a + bi) (i.e. real * complex in terms of the computer representation) then, this will come out with higher precision than (x + 0i) * (a + bi) ... because you can avoid unnecessary multiplications by zero and other such things. You can also avoid some nasty nan's that may arise in the latter case. You get to enjoy this extra-precision-and-avoid-nasty-errors because we already have built-in numerical types and std.complex.Complex defines binary operations relative to them as well as to other Complex types. I'm simply suggesting that the same opportunities to avoid those calculation errors should be available when you're dealing with purely-imaginary types. It's an implementation issue AFAICS, not a question of mathematical theory, although like you and Don I find the lack of closure for imaginary op imaginary to be very annoying. (I got round it simply by not defining opOpAssign for operations that could not be assigned back to an Imaginary type.)
Jan 02 2014
prev sibling next sibling parent "Ola Fosheim =?UTF-8?B?R3LDuHN0YWQi?= writes:
On Thursday, 2 January 2014 at 23:43:47 UTC, Lars T. Kyllingstad 
wrote:
 Not at all. ;)

I am never certain about anything related to FP. It is a very pragmatic hack, which is kind of viral. (but fun to talk about ;).
 I just think we should keep in mind why FP semantics are 
 defined the way they are.

Yes, unfortunately they are just kind-of-defined. 0.0 could represent anything from the minimum-denormal number to zero (Intel) to maximum-denormal number to zero (some other vendors). Then we have all the rounding-modes. And it gets worse with single than with double. I think the semantics of IEEE favours double over single, since detecting overflow is less important for double (it occurs seldom for doubles in practice so conflating overflow with 1.0/0.0 matters less for them than for single precision).
 Take 0.0*inf, for example.  As has been mentioned, 0.0 may 
 represent a positive real number arbitrarily close to zero, and 
 inf may represent an arbitrarily large real number. The product 
 of these is ill-defined, and hence represented by a NaN.

Yes, and it is consistent with having 0.0/0.0 evaluate to NaN. ( 0.0*(1.0/0.0) ought to give NaN as well )
 0.0+1.0i, on the other hand, represents a number which is 
 arbitrarily close to i. Multiplying it with a very large real 
 number gives you a number which has a very large imaginary 
 part, but which is arbitrarily close to the imaginary axis, 
 i.e. 0.0 + inf i. I think this is very consistent with FP 
 semantics, and may be worth making a special case in 
 std.complex.Complex.

I am too tired to figure out if you are staying within the max-min interval of potential values that can be represented (if you had perfect precision). I think that is the acid test. In order to reduce the unaccounted-for errors it is better to have a "wider" interval for each step to cover inaccuracies, and a bit dangerous if it gets more "narrow" than it should. I find it useful to try to think of floating point numbers as conceptual intervals of potential values (that get conceptually wider and wider the more you compute) and the actual FP value to be a "random" sample of that interval. For all I know, maybe some other implementations do what you suggest already, but my concern was more general than this specific issue. I think it would be a good idea to mirror a reference implementation that is widely used for scientific computation. Just to make sure that it is accepted. Imagine a team where the old boys cling to Fortran and the young guy wants D, if he can show the old boys that D produce the same results for what they do they are more likely to be impressed. Still, it is in the nature of FP that you should be able to configure and control expressions in order to overcome FP-related shortcomings. Like setting rounding-mode etc. So stuff like this ought to be handled the same way if it isn't standard practice. Not only for this stuff, but also for dealing with overflow/underflow and other "optional" aspects of FP computations.
 I agree, but there is also a lot to be said for not repeating 
 old mistakes, if we deem them to be such.

With templates you probably can find a clean way to throw in a compile-time switch for exception generation and other things that can be configured.
Jan 02 2014
prev sibling next sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 03/01/14 14:32, Stewart Gordon wrote:
 I wasn't asking for it to go beyond the existing complex implementation or any
 other.  I was proposing that the arbitrary restriction be removed so that the
 implementation we already have would work on them.

Yes, but it isn't an arbitrary restriction. Template constraints are fundamentally a promise to users about what can be expected to work. Integral types, or library types, won't work without significant modifications to the internals of the code. It would be a false promise to relax those constraints.
 FPTemporary is a template.  At the moment it's defined only for the built-in
 floating point types.  So what we really need is to define FPTemporary for
other
 types.  For int and smaller integral types, we can define it as real just as we
 do for float/double/real.  Whether it's adequate for long would be
 platform-dependent.  For other types, I suppose we can reference a property of
 the type, or just use the type itself if such a property isn't present.

 One possible idea (untested):
 -----
 template FPTemporary(F)
          if (is(typeof(-F.init + F.init * F.init - F.init))) {
      static if (isFloatingPoint!F || isIntegral!F) {
          alias real FPTemporary;
      } else static if (is(F.FPTemporary)) {
          alias F.FPTemporary FPTemporary;
      } else {
          alias F FPTemporary;
      }
 }
 -----

Yes, it ought to be possible to redefine FPTemporary (or define an alternative) to determine proper internal temporaries for any "float-esque" case. I was toying with something along the lines of, template FPTemporary(F) if (isNumeric!F || isFloatLike!F) { alias typeof(real.init * F.init) FPTemporary; } ... where isFloatLike would test for appropriate floating-point-like properties of F -- although this is probably far too simplistic. E.g. how do you handle the case of a float-like library type implemented as a class, not a struct? In any case, absent an appropriate test-case in Phobos, it would be premature to generalize the constraints for Complex or the design of FPTemporary.
 Oh yes, and it would be crazy to try and make it work for unsigned integer
 types.  But even if we were to resolve the issues with FPTemporary and that, it
 would still fall under my earlier suggestion of making it so that if people
want
 to use Complex on an unsupported type then they can explicitly suppress the
type
 restriction, but should understand that it might not work properly.

People who want to use Complex on an unsupported type can quite readily copy-paste the code and remove the constraints, if that's what they want to do. I think that's better than giving them an option which is essentially an invitation to shoot themselves in the foot, and which has very little chance of actually working. It doesn't matter if you document it as "This might not work", by providing the option you are still essentially saying, "This is an OK way to use the type." I think that's essentially an encouragement of bad code and a violation of the design principle that the easy thing to do should be the right thing to do.
 Really, I was just thinking that somebody who wants a quick-and-dirty
 hypercomplex number implementation for some app might try to do it that way.

I understand that, but quick-and-dirty solutions are often bad ones, and in this case, it just wouldn't work given the internals of Complex. If you would like to revise std.complex to support this approach, I'm sure your pull request will be considered carefully, but personally I don't see it as an effective way to pursue hypercomplex number support when there are other options on the table.
Jan 03 2014
prev sibling parent Joseph Rushton Wakeling <joseph.wakeling webdrake.net> writes:
On 03/01/14 21:21, Stewart Gordon wrote:
 How can isFloatLike be implemented?

I'm not sure. It's something that needs to be thought about and of course it also depends on whether you want it to test just for basic properties, or whether it is supported by mathematical operations (sin, cos, abs, exp, etc.). I think testing for basic properties should be enough, because if mathematical functions are needed but not supported, there will be a compilation error anyway.
 And how can we test for bigint or Galois field types?

For BigInt, it's conceivable to define an isIntegerLike template (David Simcha did this for his std.rational) that will handle library as well as built-in integral types. For Galois field types, I suggest that there needs to be an implementation before one talks about how to support them in std.complex.
 There are already violations of this principle in D, such as being able to cast
 away const.

Casting away const has valid applications and is a bit different from allowing the user to manually ignore template constraints on a library type, which will most likely almost always lead to problematic behaviour. I'm not aware of any library type in Phobos that does this, and if you _really_ want to override the template constraints, you can always copy and modify the code. Then, if it turns out to work well, you can submit patches to allow that case in Phobos too.
 Addition, subtraction and multiplication would work.

They wouldn't, because when you relaxed the template constraints to allow Complex!(Complex!T), the code would fail to compile. And I don't think you should correct that by stripping out basic arithmetic operations.
 So the programmer could just copy the code and reimplement division and
exponentiation so that they work
 (or just get rid of them if they aren't needed).

As I keep saying, you're asking for extra complication to be added to a type that supports its intended (probably the vast majority of) use cases well, for the sake of a niche use case _that can be implemented without any problem as an entirely independent type_. What you perceive as conceptual/notational elegance -- Complex!(Complex!T) -- isn't worth it.
Jan 03 2014