How odd...

On to the complex number case. From the standard for Complex: # log z = log (magnitude z) :+ phase z # phase (0 :+ 0) = 0 This is a special case for the phase of zero. # (x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x')

...

log (0 :+ 0) (-Infinity) :+ 0.0

...

log (0 :+ 0) * 2 (-Infinity) :+ NaN

Which is the source of the problem. The imaginary part involves multiplying (-Infinity) by the imaginary part of 2 (i.e. 0), which is NaN. So no, its not a bug, its according to the standard. Whether the standard ought to be modified for multiplication by numbers of the form (x :+ 0) is another question. Arguably the correct value in the last case should have been (-Infinity) :+ 0.0 on the grounds that (x :+ y) * 2 should be equal to (x * 2 :+ y * 2). In most cases this holds, but if one of the complex multiplicands contains an infinity then you get a NaN instead. While working through this I also came across the following case which technically is a bug:

...

0 ** 0 1.0

...

exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case. Paul.

Andrew Coppin

11:38 a.m.

Paul Johnson wrote:

...

Andrew Coppin wrote:

...
...
0**2 0

...
(0 :+ 0)**2 NaN :+ NaN

(Is this a bug?) According to the Standard Prelude, # x ** y = exp (log x * y)

I had a feeling this would be the cause.

...

...
log 0 -Infinity

Oh. So... since when does Haskell know about infinity? BTW, I recently had some code like this: foo x | x < 0 = ... | x == 0 = ... | x > 0 = ... I was most perplexed when I got a "non-exhaustive patterns" exception... It turns out there was a NaN in there. I forget about that.

...

...
exp (log 0 * 2) 0.0

Well that's interesting. I did wonder why it *doesn't* break in the real case...

...

On to the complex number case. From the standard for Complex:

# log z = log (magnitude z) :+ phase z

# phase (0 :+ 0) = 0 This is a special case for the phase of zero.

# (x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x')

...
log (0 :+ 0) (-Infinity) :+ 0.0

...
log (0 :+ 0) * 2 (-Infinity) :+ NaN

Which is the source of the problem. The imaginary part involves multiplying (-Infinity) by the imaginary part of 2 (i.e. 0), which is NaN.

Um... why would infinity * 0 be NaN? That doesn't make sense...

...

So no, its not a bug, its according to the standard.

So I'm the only person who was expecting zero squared to be zero? (IMHO the standard should try to implement mathematical operations in a mathematically sensible way...)

...

While working through this I also came across the following case which technically is a bug:

...
0 ** 0 1.0

...
exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case.

Now that really *is* odd...

Lennart Augustsson

11:47 a.m.

Haskell doesn't know much about infinity, but Haskell implementations are allowed to use IEEE floating point which has infinity. And to get things right, there needs to be a few changes to the library to do the right thing for certain numbers, this is not news. In fact I filed a bug report a while back about it. -- Lennart On 8/4/07, Andrew Coppin wrote:

...

Paul Johnson wrote:

...
Andrew Coppin wrote:

...
...
0**2 0

...
(0 :+ 0)**2 NaN :+ NaN

(Is this a bug?) According to the Standard Prelude, # x ** y = exp (log x * y)

I had a feeling this would be the cause.

...
...
log 0 -Infinity

Oh. So... since when does Haskell know about infinity?

BTW, I recently had some code like this:

foo x | x < 0 = ... | x == 0 = ... | x > 0 = ...

I was most perplexed when I got a "non-exhaustive patterns" exception... It turns out there was a NaN in there. I forget about that.

...
...
exp (log 0 * 2) 0.0

Well that's interesting. I did wonder why it *doesn't* break in the real case...

...
On to the complex number case. From the standard for Complex:

# log z = log (magnitude z) :+ phase z

# phase (0 :+ 0) = 0 This is a special case for the phase of zero.

# (x:+y) * (x':+y') = (x*x'-y*y') :+ (x*y'+y*x')

...
log (0 :+ 0) (-Infinity) :+ 0.0

...
log (0 :+ 0) * 2 (-Infinity) :+ NaN

Which is the source of the problem. The imaginary part involves multiplying (-Infinity) by the imaginary part of 2 (i.e. 0), which is NaN.

Um... why would infinity * 0 be NaN? That doesn't make sense...

...
So no, its not a bug, its according to the standard.

So I'm the only person who was expecting zero squared to be zero? (IMHO the standard should try to implement mathematical operations in a mathematically sensible way...)

...
While working through this I also came across the following case which technically is a bug:

...
0 ** 0 1.0

...
exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case.

Now that really *is* odd...

_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

Andrew Coppin

12:16 p.m.

Lennart Augustsson wrote:

...

Haskell doesn't know much about infinity, but Haskell implementations are allowed to use IEEE floating point which has infinity. And to get things right, there needs to be a few changes to the library to do the right thing for certain numbers, this is not news. In fact I filed a bug report a while back about it.

It's news to me that the IEEE floating point standard defines infinity. (NaN I knew about...)

Paul Johnson

1:06 p.m.

Andrew Coppin wrote:

...

It's news to me that the IEEE floating point standard defines infinity. (NaN I knew about...) See http://en.wikipedia.org/wiki/IEEE_754

Paul.

Paul Johnson

1:16 p.m.

Also see http://hal.archives-ouvertes.fr/hal-00128124 before you start on any serious numerical software. Paul.

Simon Peyton-Jones

9 Aug 9 Aug

5:53 a.m.

Indeed! Getting the library numerics to do the Right Thing is something that can only be done by people who know about numerics. People who build compilers aren't, alas. It's quite a specialised subject, and very easy to screw up. And there's performance to worry about too in the common case when the corner cases don't appear. So if there are folk out there who care about getting numerics correct, we would welcome your direct involvement. Look at the code, make it right, send patches to libraries@haskell.orgmailto:libraries@haskell.org. It's all open source, and you'll be benefiting lots of people. Simon From: haskell-cafe-bounces@haskell.org [mailto:haskell-cafe-bounces@haskell.org] On Behalf Of Lennart Augustsson Sent: 04 August 2007 16:47 To: Andrew Coppin Cc: haskell-cafe@haskell.org Subject: Re: [Haskell-cafe] How odd... Haskell doesn't know much about infinity, but Haskell implementations are allowed to use IEEE floating point which has infinity. And to get things right, there needs to be a few changes to the library to do the right thing for certain numbers, this is not news. In fact I filed a bug report a while back about it. -- Lennart On 8/4/07, Andrew Coppin mailto:andrewcoppin@btinternet.com> wrote: Paul Johnson wrote:

...

Andrew Coppin wrote:

...
...
0**2 0

...
(0 :+ 0)**2 NaN :+ NaN

(Is this a bug?) According to the Standard Prelude, # x ** y = exp (log x * y)

Paul Johnson

4 Aug 4 Aug

1:07 p.m.

Andrew Coppin wrote:

...

Paul Johnson wrote:

...
...
log 0 -Infinity Oh. So... since when does Haskell know about infinity? I should have mentioned that the underlying platform in my case is an Intel P4. Haskell does not specify a floating point implementation; the assumption is that it uses whatever the platform provides because anything else would be horribly inefficient. The P4 implements IEEE floating point, which as Andrew pointed out includes Infinity, -Infinity and NaN as special cases of floating point values. It also seems to distinguish between 0.0 and (-0.0), although they are still equal. For instance:

...

(-0.0) -0.0

...

1.0 / 0.0 Infinity

...

1.0 / (-0.0) -Infinity

(Aside: should Infinity etc. be included in the Haskell standard? Or possibly in a Data.Numeric.IEEE extension? They look like constructors, and it would be nice to be able to pattern match on them.) So if you tried these experiments on a non-IEEE platform then you would get different results. You might get exceptions or just wierd big numbers. ISTR someone posted results from a Sun Sparc along these lines some time ago.

...

foo x | x < 0 = ... | x == 0 = ... | x > 0 = ... I was most perplexed when I got a "non-exhaustive patterns" exception... It turns out there was a NaN in there. I forget about that. Nasty. I'll have to bear that one in mind myself.

You can detect infinities by equality, but not NaN.

...

let inf = (1.0 / 0.0) let nan = inf * 0.0 inf == inf True nan == nan False

...

...
...
exp (log 0 * 2) 0.0 Well that's interesting. I did wonder why it *doesn't* break in the real case... I haven't perused the IEEE standard, but I expect it defines something like this:

Infinity * 0 = NaN Infinity * _ = Infinity exp Infinity = Infinity exp (-Infinity) = 0

...

Um... why would infinity * 0 be NaN? That doesn't make sense... Infinity times anything is Infinity. Zero times anything is zero. So what should Infinity * zero be? There isn't one right answer. In this case the "morally correct" answer is zero, but in other contexts it might be Infinity or even some finite number other than zero.

Consider 0.0 / 0.0, which also evaluates to NaN. What should its value be? If you take the limit of (x / x) as x -> 0 then the right answer is 0. On the other hand if you take the limit of (0 / x) as x -> 0 then the right answer is infinity. Worse yet, if you take the limit of (2x / x) as x-> 0 then the right answer is 2. You can pick any finite or infinite value you like. So the only possible answer is NaN.

...

...
So no, its not a bug, its according to the standard.

So I'm the only person who was expecting zero squared to be zero? (IMHO the standard should try to implement mathematical operations in a mathematically sensible way...)

...

...
While working through this I also came across the following case which technically is a bug:

...
0 ** 0 1.0

...
exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case.

Now that really *is* odd... When I said "built in" I meant "built in to the hardware". This is

It does *try*. I'm not sure if IEEE arithmetic was actually defined back in 98. It certainly wasn't widely implemented. There might well be a case for revisiting the standard to allow for IEEE values, but you can't mandate them because not all platforms support IEEE arithmetic. probably another special case defined in the IEEE standard, which is not the same as the Haskell 98 definition. The reason why the IEEE standard was defined in the first place was that floating point software portability was being seriously limited by these corner cases. You would get some numerical code working on a Sun, and then find it broke on a PC because one of them defined (0.0 / 0.0) as 1 and the other defined it as 0. Worse yet, you might find it worked on Intel chips but not IBM or AMD. Programmers also had to code explicit checks for division by zero and implement their own versions of Infinity and NaN in cases where they might appear, which cluttered up the code and slowed down execution. One way out of this morass would be to define Haskell floating point arithmetic as IEEE standard, and document non-conformance for certain platforms. In the long term that is probably the way forwards (do any currently sold chips *not* implement IEEE?). It would also let us include Infinity and NaN as constructors. However it would lead to significant problems when compiling code that used these values on non-IEEE platforms. What do you do then? Paul.

Lennart Augustsson

1:21 p.m.

The Haskell type class RealFloat has a reasonable number of operations to test for NaN, denormalized numbers, etc. You can also ask if the implementation uses IEEE. -- Lennart On 8/4/07, Paul Johnson wrote:

...

Andrew Coppin wrote:

...
Paul Johnson wrote:

...
...
log 0 -Infinity Oh. So... since when does Haskell know about infinity? I should have mentioned that the underlying platform in my case is an Intel P4. Haskell does not specify a floating point implementation; the assumption is that it uses whatever the platform provides because anything else would be horribly inefficient. The P4 implements IEEE floating point, which as Andrew pointed out includes Infinity, -Infinity and NaN as special cases of floating point values. It also seems to distinguish between 0.0 and (-0.0), although they are still equal. For instance:

...
(-0.0) -0.0

...
1.0 / 0.0 Infinity

...
1.0 / (-0.0) -Infinity

(Aside: should Infinity etc. be included in the Haskell standard? Or possibly in a Data.Numeric.IEEE extension? They look like constructors, and it would be nice to be able to pattern match on them.)

So if you tried these experiments on a non-IEEE platform then you would get different results. You might get exceptions or just wierd big numbers. ISTR someone posted results from a Sun Sparc along these lines some time ago.

...
foo x | x < 0 = ... | x == 0 = ... | x > 0 = ... I was most perplexed when I got a "non-exhaustive patterns" exception... It turns out there was a NaN in there. I forget about that. Nasty. I'll have to bear that one in mind myself.

You can detect infinities by equality, but not NaN.

...
let inf = (1.0 / 0.0) let nan = inf * 0.0 inf == inf True nan == nan False

...
...
...
exp (log 0 * 2) 0.0 Well that's interesting. I did wonder why it *doesn't* break in the real case... I haven't perused the IEEE standard, but I expect it defines something like this:

Infinity * 0 = NaN Infinity * _ = Infinity exp Infinity = Infinity exp (-Infinity) = 0

...
Um... why would infinity * 0 be NaN? That doesn't make sense... Infinity times anything is Infinity. Zero times anything is zero. So what should Infinity * zero be? There isn't one right answer. In this case the "morally correct" answer is zero, but in other contexts it might be Infinity or even some finite number other than zero.

Consider 0.0 / 0.0, which also evaluates to NaN. What should its value be? If you take the limit of (x / x) as x -> 0 then the right answer is 0. On the other hand if you take the limit of (0 / x) as x -> 0 then the right answer is infinity. Worse yet, if you take the limit of (2x / x) as x-> 0 then the right answer is 2. You can pick any finite or infinite value you like. So the only possible answer is NaN.

...
...
So no, its not a bug, its according to the standard.

So I'm the only person who was expecting zero squared to be zero? (IMHO the standard should try to implement mathematical operations in a mathematically sensible way...)

...
...
While working through this I also came across the following case which technically is a bug:

...
0 ** 0 1.0

...
exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case.

Now that really *is* odd... When I said "built in" I meant "built in to the hardware". This is

It does *try*. I'm not sure if IEEE arithmetic was actually defined back in 98. It certainly wasn't widely implemented. There might well be a case for revisiting the standard to allow for IEEE values, but you can't mandate them because not all platforms support IEEE arithmetic. probably another special case defined in the IEEE standard, which is not the same as the Haskell 98 definition.

The reason why the IEEE standard was defined in the first place was that floating point software portability was being seriously limited by these corner cases. You would get some numerical code working on a Sun, and then find it broke on a PC because one of them defined (0.0 / 0.0) as 1 and the other defined it as 0. Worse yet, you might find it worked on Intel chips but not IBM or AMD. Programmers also had to code explicit checks for division by zero and implement their own versions of Infinity and NaN in cases where they might appear, which cluttered up the code and slowed down execution.

One way out of this morass would be to define Haskell floating point arithmetic as IEEE standard, and document non-conformance for certain platforms. In the long term that is probably the way forwards (do any currently sold chips *not* implement IEEE?). It would also let us include Infinity and NaN as constructors. However it would lead to significant problems when compiling code that used these values on non-IEEE platforms. What do you do then?

Paul.

_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

Ian Lynagh

2:41 p.m.

On Sat, Aug 04, 2007 at 04:38:08PM +0100, Andrew Coppin wrote:

...

BTW, I recently had some code like this:

foo x | x < 0 = ... | x == 0 = ... | x > 0 = ...

I was most perplexed when I got a "non-exhaustive patterns" exception... It turns out there was a NaN in there. I forget about that.

If you use guards then you can get non-exhaustive pattern warnings anyway, e.g.: myNot x | x == True = False | x == False = True gives Warning: Pattern match(es) are non-exhaustive In the definition of `myNot': Patterns not matched: _ as it is (in general) undecidable whether all patterns are covered.

...

...
While working through this I also came across the following case which technically is a bug:

...
0 ** 0 1.0

...
exp (log 0 * 0) NaN

I suspect that GHCi is using a built-in exponentiation operator that doesn't quite conform to the standard in this case.

Hmm, according to the report the Floating class has a default method x ** y = exp (log x * y) but the Double instance is just instance Floating Double where ... so could define something different. 6.4.3 says 0**y is undefined. Thanks Ian

Neil Mitchell

2:53 p.m.

Hi If you just use Catch (http://www-users.cs.york.ac.uk/~ndm/catch/):

...

...
foo x | x < 0 = ... | x == 0 = ... | x > 0 = ...

This gives an error. Something identical to this code is in Data.FiniteMap, and indeed, when using floats and NaN's (or just silly Ord classes) you can cause Data.FiniteMap to pattern match error. See section 6.3 of the draft paper on the Catch website for details.

...

myNot x | x == True = False | x == False = True

This is determined to be exhaustive.

...

as it is (in general) undecidable whether all patterns are covered.

In general, yes. But its possible to do a lot better than GHC's current warnings (I'm not saying its worth changing, though) Thanks Neil

Albert Y. C. Lai

2:59 p.m.

Andrew Coppin wrote:

...

...
0^2 0

...
(0 :+ 0)^2 0 :+ 0

...
0**2 0

...
(0 :+ 0)**2 NaN :+ NaN

There is nothing wrong AFAIK. (How much do I know? BSc in math, went through classes on real analysis and complex analysis.) There is no reason to expect complex ** to agree with real **. Real x**y is first defined for natural y (agreeing with x^y), then extend to integer y (agreeing with x^^y), then extend to rational y (taking nth root when y = m/n), then extend to real y by continuity wherever possible. You can still expect real 0**2 = 0^2 = 0. Complex x**y involves picking a phase angle of x. "Phase angle" is an ill notion for 0. Complex 0**y is better left undefined. You said

...

So I'm the only person who was expecting zero squared to be zero?

Are you trying to be sensational, or are you going hyperbole? If you want zero squared, you're welcome to use 0^2. Complex 0**2 does not have to be related to zero squared. You said

...

(IMHO the standard should try to implement mathematical operations in a mathematically sensible way...)

But AFAIK it is already a mathematically sensible way - it is what I learned from my math classes. Perhaps you mean highschoolly sensible. I understand that highschool math is a bit simple-minded, e.g., you can greatly confuse someone by 1 = ((-1)**2)**0.5 = ((-1)**0.5)**2) = i**2 = -1 "So I'm the only one expecting x square-root squared to be x?" Thank God complex ** is supposed to be different from real **.

Dan Piponi

3:47 p.m.

On 8/4/07, Albert Y. C. Lai wrote:

...

There is no reason to expect complex ** to agree with real **.

There's every reason. It is standard mathematical practice to embed the integers in the rationals in the reals in the complex numbers and it is nice to have as many functions as possible respect that embedding. In particular, it's nice to get a kind of commutativity when we start casting between different numeric types. There's an old paper on this subject by Reynolds (http://lambda-the-ultimate.org/node/2078) though I'm more familiar with the treatment in Pierce's "Basic Category for Computer Scientists". If you look through a variety of complex analysis books you'll find a number of different approaches to handling 0^x for complex x. Many agree that this is defined for Re(x)>0. Some say it is defined for integer x>0. Some say it's not defined at all, and then in the following paragraphs talk about x^2 where x takes values in the entire complex plane. So it seems to me that it is entirely reasonable to hope that (0:+0)**2=0:+0.

...

Complex x**y involves picking a phase angle of x. "Phase angle" is an ill notion for 0. Complex 0**y is better left undefined.

There is no right answer here, it's a matter of taste. Among mathematicians there doesn't seem to be a uniform consensus. But I think that in computing a commutativity principle is nice to have, so I think we should have x**fromInteger y = x^y, where defined, so that, in some sense, (^) is embedded in (**). -- Dan

Stephen Forrest

4:35 p.m.

On 8/4/07, Dan Piponi wrote:

...

On 8/4/07, Albert Y. C. Lai wrote:

...
There is no reason to expect complex ** to agree with real **.

There's every reason. It is standard mathematical practice to embed the integers in the rationals in the reals in the complex numbers and it is nice to have as many functions as possible respect that embedding.

A example I have seen before that illustrates some the difficulties with preserving such behaviour is (-1)^(1/3). Of course, taking the nth root is multi-valued, so if you're to return a single value, you must choose a convention. Many implementations I have seen choose the solution with lowest argument (i.e. the first solution encounted by a counterclockwise sweep through the plane starting at (1,0).) With this interpretation, (-1)^(1/3) = 0.5 + sqrt(3)/2 * i. If you go with the real solution (-1) you might need to do so carefully in order to preserve other useful properties of ^, like continuity. Steve

Dan Piponi

5:06 p.m.

On 8/4/07, Stephen Forrest wrote:

...

Of course, taking the nth root is multi-valued, so if you're to return a single value, you must choose a convention. Many implementations I have seen choose the solution with lowest argument (i.e. the first solution encounted by a counterclockwise sweep through the plane starting at (1,0).)

This is one approach but not the one used by some languages. Another way of looking at it is this: a function like log is multivalued as you pass around the origin. So what you can do is 'cut' the complex plane so that you're not allowed to cross the cuts. Now you have a region on which a single-valued function can be unambiguously defined. It's fairly standard practice, when documenting functions of a complex variable, to specify precisely which 'branch cuts' are being used. Here's a quote from the Mathematica documentation describing their Log function: "Log[z] has a branch cut discontinuity in the complex z plane running from -infinity to 0".

...

With this interpretation, (-1)^(1/3) = 0.5 + sqrt(3)/2 * i. If you go with the real solution (-1) you might need to do so carefully in order to preserve other useful properties of ^, like continuity.

You can guarantee this by making sure you make the right 'cuts'. -- Dan

Jacques Carette

11 Aug 11 Aug

10:21 a.m.

[Sorry for the long quote, but context is important] Dan Piponi wrote:

...

It's fairly standard practice, when documenting functions of a complex variable, to specify precisely which 'branch cuts' are being used. Here's a quote from the Mathematica documentation describing their Log function: "Log[z] has a branch cut discontinuity in the complex z plane running from -infinity to 0".

...
With this interpretation, (-1)^(1/3) = 0.5 + sqrt(3)/2 * i. If you go with the real solution (-1) you might need to do so carefully in order to preserve other useful properties of ^, like continuity.

You can guarantee this by making sure you make the right 'cuts'.

If only it were that simple. Up until rather recently, it was not even known if there was a set of "right cuts" for ln and all the arc-trig functions that would work "together" properly. See `"According to Abramowitz and Stegun" or arccoth needn't be uncouth', by Corless, Jeffrey, Watt and Davenport, available from http://citeseer.ist.psu.edu/corless00according.html or http://www.apmaths.uwo.ca/~djeffrey/Offprints/couth.pdf This is by no means a solved problem. For example, there is still no consensus as to where the branch cuts for the various versions of the LegendreQ function should go. Jacques http://www.apmaths.uwo.ca/%7Edjeffrey/Offprints/couth.pdf

6495

Age (days ago)

6502

Last active (days ago)

List overview

Download

16 comments

10 participants

participants (10)

Albert Y. C. Lai
Andrew Coppin
Dan Piponi
Ian Lynagh
Jacques Carette
Lennart Augustsson
Neil Mitchell
Paul Johnson
Simon Peyton-Jones
Stephen Forrest

How odd...

Andrew Coppin

Andrew Coppin

Andrew Coppin

Dan Piponi

Dan Piponi

tags

participants (10)