Thread: [pure-lang-users] complex 0

SourceForge Headquarters 1320 Columbia Street Suite 310 San Diego, CA 92101 +1 (858) 422-6466

While piddling with some math problems, I noticed the following
behavior. 

 > (2.0+:1.0)/0.0;
 inf+:inf
 > (2.0+:1.0)/(0.0+:0.0);
 nan+:nan

Um. Should this be consistent?

e.r.

Eddie Rucker wrote:
>  > (2.0+:1.0)/0.0;
>  inf+:inf
>  > (2.0+:1.0)/(0.0+:0.0);
>  nan+:nan
> 
> Um. Should this be consistent?

Maybe, not sure. AFAICS, there are three possible perspectives:

- (2.0+:1.0)/0.0 has an exact 0 in the imaginary part of the divisior, 
so it's actually a vector multiplied with an infinite scalar, whereas in 
the case of (2.0+:1.0)/(0.0+:0.0) it's an inexact zero so the general 
formula must be used which yields nan+:nan because a zero meets a pole 
in both components. (OTOH, if we take this view then probably 
(2.0+:1.0)/(0.0+:0) should yield the same result as the former, whereas 
right now it yields the same result as the latter.)

- The divisor should be promoted to a complex 0.0+:0.0 in the first 
case, so you get nan+:nan in either case.

- The divisior should be promoted to a real 0.0 in the second case, 
yielding inf+:inf in either case. That's what I get with mzscheme, but 
it doesn't make any real sense (pun intended) to me because 0.0+:0.0 is 
*not* the same as 0.0+:0 (or just 0.0) in the IEEE 754 sense (the 
imaginary zero might as well be just an underflow).

I haven't checked the Goldberg paper, though. Maybe there's a "more 
defined" formula in the complex divisor case?

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

On Tue, 2008-09-09 at 00:31 +0200, Albert Graef wrote:
> Eddie Rucker wrote:
> >  > (2.0+:1.0)/0.0;
> >  inf+:inf
> >  > (2.0+:1.0)/(0.0+:0.0);
> >  nan+:nan
> > 
> > Um. Should this be consistent?
> 
> Maybe, not sure. AFAICS, there are three possible perspectives:
> 
> - (2.0+:1.0)/0.0 has an exact 0 in the imaginary part of the divisior, 
> so it's actually a vector multiplied with an infinite scalar, whereas in 
> the case of (2.0+:1.0)/(0.0+:0.0) it's an inexact zero so the general 
> formula must be used which yields nan+:nan because a zero meets a pole 
> in both components. (OTOH, if we take this view then probably 
> (2.0+:1.0)/(0.0+:0) should yield the same result as the former, whereas 
> right now it yields the same result as the latter.)

In Pure, ln 0 => -inf and ln 0.0 -> -inf.  The former is an exact zero
so in real math ln 0 should be undefined. The second is a number close
to zero which yields -inf. So I'm thinking the exact 0 for the imaginary
part of 0.0 should also be promoted to double just like in ln 0. If I'm
wrong then I'm Ok with that.

> - The divisor should be promoted to a complex 0.0+:0.0 in the first 
> case, so you get nan+:nan in either case.
> 
> - The divisior should be promoted to a real 0.0 in the second case, 
> yielding inf+:inf in either case. That's what I get with mzscheme, but 
> it doesn't make any real sense (pun intended) to me because 0.0+:0.0 is 
> *not* the same as 0.0+:0 (or just 0.0) in the IEEE 754 sense (the 
> imaginary zero might as well be just an underflow).

Hm. if we start worrying about 0.0+:0 vs 0.0+:0.0, then should we also
worry about ln 0 vs ln 0.0 too? I think I've already mentioned that
mzscheme says that (log 0) is undefined and (log 0.0) => -inf.0.
Am I'm flogging a dead horse here? Oh woe woe woe is me!

e.r.

Albert Graef scripsit:

> - The divisor should be promoted to a complex 0.0+:0.0 in the first 
> case, so you get nan+:nan in either case.

Division with slash effectively promotes its arguments to floats, so
0.0 should be promoted to 0.0+:0.0.

-- 
As you read this, I don't want you to feel      John Cowan
sorry for me, because, I believe everyone       co...@cc...
will die someday.                               http://www.ccil.org/~cowan
        --From a Nigerian-type scam spam

On Mon, 2008-09-08 at 20:14 -0400, John Cowan wrote:
> Albert Graef scripsit:
> 
> > - The divisor should be promoted to a complex 0.0+:0.0 in the first 
> > case, so you get nan+:nan in either case.
> 
> Division with slash effectively promotes its arguments to floats, so
> 0.0 should be promoted to 0.0+:0.0.

Allrighty then! I know Albert already has enough flies to swat without
me fanning the turds :( Sorry.

e.r.

John Cowan wrote:
> Division with slash effectively promotes its arguments to floats, so
> 0.0 should be promoted to 0.0+:0.0.

Good point! In fact gcc seems to do that, too (I get nan+:nan for /0.0 
just as well as for /(0.0+I*0.0).

However, that's *not* what ISO/IEC9899, Annex G [1] says (see G.5.1, 
especially Example 2). It requires an inf+I*inf result in the /0.0 case, 
if I read it right. That has the "informative" sticker (i.e., 
"recommended" but not "normative"), so implementations may differ. 
AFAICT, my implementation in math.pure conforms with Annex G (at least 
the multiplicative operators, still have to check the others). So I 
guess I'll just leave it at that for now.

[1] I'm referring to the 2005 draft, couldn't find anything newer. 
Here's the URL: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

Albert Graef wrote:
> However, that's *not* what ISO/IEC9899, Annex G [1] says (see G.5.1, 
> especially Example 2). It requires an inf+I*inf result in the /0.0 case, 
> if I read it right.

Uhmm, of course it also requires a inf+:inf result in the /(0+:0) case, 
I'll have to fix that in math.pure...

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

Albert Graef wrote:
> [1] I'm referring to the 2005 draft, couldn't find anything newer. 
> Here's the URL: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf

Here's a newer version, linked to from the Wikipedia C99 page:
http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

I know I'm the one that brought this up, but how about we revisit the
issue later when we have time to think about this. IIRC you needed to
finish Pd for you classes. Later, we can compare Pure results to GSL's
results. I think consistency with GSL should be a must.

I'll note math things that look funny and we can hash them later. I'll
just report other types of errors for now.

e.r.

On Thu, 2008-09-11 at 12:34 +0200, Albert Graef wrote:
> I'm beginning to wonder whether it's really worth the hassle to follow 
> the ISOC99 recommendations concerning the treatment of infinities and 
> NaNs for complex * and /. The algorithms in Annex G do look like a kind 
> of kludge to me, and they employ C99 functions for extracting and 
> scaling the exponents of floating point numbers which might not be 
> available on some systems (such as Windows), at least not directly. 
> Hence special support in the runtime would be needed to implement these 
> in a portable way across all platforms.
> 
> That seems like an awful lot of kludges to just make those nan+:nan's go 
> away, considering that Pure also offers the polar representation where 
> you'll get the correct infinite results for division by zero and similar 
> cases, as long as the phase angles are finite.
> 
> So actually John's suggestion seems to be the most reasonable: just 
> promote real operands of * and / to complex, then the results will at 
> least be consistent. Or does anyone here really need the ISOC 
> recommended behaviour?
> 
> Albert
> 

Eddie Rucker wrote:
> I know I'm the one that brought this up, but how about we revisit the
> issue later when we have time to think about this.

Well, in any case it was bad to have these two cases return inconsistent 
results, so I fixed that now in the way I described. Making it conform 
with the ISOC99 recommendations is now just a matter of editing the two 
rules for complex * and /, if we lateron decide that we need that.

> Later, we can compare Pure results to GSL's
> results. I think consistency with GSL should be a must.

The only way to achieve 100% compatibility with GSL there would be to 
actually use the complex multiplication and division routines employed 
by GSL itself. But presumably GSL just uses whatever complex arithmetic 
the C compiler provides (or Fortran, in the case of BLAS). In the case 
of the GNU compilers (which are the ones we use on all major platforms 
except the *BSDs) the results should be compatible with what math.pure 
provides now.

An alternative would be to actually use the complex double routines 
provided by the C runtime and library, and add marshalling of complex 
results to the C interface. That would be useful to have anyway, for 
interfacing to GSL. I'll have a look at that in the next release.

> I'll note math things that look funny and we can hash them later. I'll
> just report other types of errors for now.

Ok, thanks.

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

On Tue, 2008-09-09 at 13:00 +0200, Albert Graef wrote:

> Yes, you are, I'm afraid. :) The Scheme way doesn't make much sense 
> there because log is always an inexact operation.

What about ln 1? Scheme has a different philosophy from Pure: all
operations *try* to return an exact value if the operation's argument(s)
are exact. 

(log 1) => 0 ; an exact 0 not an approximated 0.0
(sin 0) => 0
(cos 0) => 1 ; an exact 1
(sqrt 4) => 2 ; an exact 2
(sqrt -4) => 0+2i ; no approximations on the img part.
(/ 1 3) => 1/3 
(expt 16 1/2) => 4 
(expt 8 1/3) => 2.0 ; uh oh? My Casio gives the exact value of 2 here.
What does you HP give?

I'm not disagreeing with your reasoning that Pure should promote its
argument(s) to inexact values like C does though. I think I've totally
chased the rabbit away from my original question which was about
consistency for which you found the answer in a later post. My arm is
getting tired of flogging. I'll shut up and go back to work now ;-) 

e.r.

On Tue, 2008-09-09 at 10:30 -0400, John Cowan wrote:
> R5RS only requires this of a particular list of standard procedures,
> which does not include / or the transcendental functions.  R6RS, which
> requires support for arbitrary ratios as R5RS does not, does require /
> to produce exact results on exact arguments, and expt to produce exact
> results if the first argument is an exact real value and the second
> argument an exact integer.
> 
> R6RS also explicitly says that sqrt, exp, log, sin, cos, tan, asin, acos,
> atan, expt, make-polar, magnitude, and angle may (but are not required
> to) return inexact results even when given exact arguments.

I got this wrong then. Apologies to Albert. I assumed all Scheme
implementations implemented arithmetic like MzScheme. You'd think by now
I would not assume anything ;-)

e.r.

Eddie Rucker scripsit:

> I got this wrong then. Apologies to Albert. I assumed all Scheme
> implementations implemented arithmetic like MzScheme. You'd think by now
> I would not assume anything ;-)

R5RS Scheme implementations are really separate languages with a common
small core, not like C or Ruby or Python implementations.

-- 
Mark Twain on Cecil Rhodes:                    John Cowan
I admire him, I freely admit it,               http://www.ccil.org/~cowan
and when his time comes I shall                co...@cc...
buy a piece of the rope for a keepsake.

Eddie Rucker wrote:
> I got this wrong then. Apologies to Albert.

There's no reason to apologize, and there's no easy right or wrong in 
these questions, that's what makes them so difficult. Of course we don't 
want to violate mathematical truths (and even that's not guaranteed with 
floating point numbers!), but apart from that there's still a lot of 
leeway in which formulas and algebraic identities we're going to use.

The acclaimed but essentially vacuous "principle of the least surprise" 
may yield different results here, depending on whether you look at it 
from the pov of functional analysis, numerics, or just plain software 
engineering.

That's why I try to stick to established standards, but even ISO/IEC9899 
says that they're not sure about their "informational" design decisions 
in that area either, that's why they tagged them that way. I can 
understand that, I think that I'm a reasonably well-educated 
mathematician, but I'm not sure about the right way to implement these 
things either. :)

For the time being, I'll try to follow the ISO/IEC9899 recommendations, 
so if you notice anything which is at odds with these, please let me know!

Cheers,
Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

Eddie Rucker wrote:
> In Pure, ln 0 => -inf and ln 0.0 -> -inf.

Yes, that's because int args of those floating point operations are 
generally promoted to double. It's the most sensible way to do it. The 
same is true for / (which is *always* an inexact floating point division 
in Pure, just like in C).

> Hm. if we start worrying about 0.0+:0 vs 0.0+:0.0, then should we also
> worry about ln 0 vs ln 0.0 too? I think I've already mentioned that
> mzscheme says that (log 0) is undefined and (log 0.0) => -inf.0.
> Am I'm flogging a dead horse here?

Yes, you are, I'm afraid. :) The Scheme way doesn't make much sense 
there because log is always an inexact operation.

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

Albert Graef scripsit:

> Yes, that's because int args of those floating point operations are 
> generally promoted to double. It's the most sensible way to do it. The 
> same is true for / (which is *always* an inexact floating point division 
> in Pure, just like in C).

If only.  In fact, int/int in C is truncating integer division, but that's
okay, because C is obviously broken there.  I know Pure is supposed to
"do it the C way", but that's carrying things too far.

-- 
How they ever reached any  conclusion at all    <co...@cc...>
is starkly unknowable to the human mind.        http://www.ccil.org/~cowan
        --"Backstage Lensman", Randall Garrett

John Cowan wrote:
> If only.  In fact, int/int in C is truncating integer division, but that's
> okay, because C is obviously broken there.  I know Pure is supposed to
> "do it the C way", but that's carrying things too far.

Right, thanks for correcting me. s/C/Pascal/ (which does it in a much 
more sensible way than C IMHO.)

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

Eddie Rucker scripsit:

> What about ln 1? Scheme has a different philosophy from Pure: all
> operations *try* to return an exact value if the operation's argument(s)
> are exact. 

R5RS only requires this of a particular list of standard procedures,
which does not include / or the transcendental functions.  R6RS, which
requires support for arbitrary ratios as R5RS does not, does require /
to produce exact results on exact arguments, and expt to produce exact
results if the first argument is an exact real value and the second
argument an exact integer.

R6RS also explicitly says that sqrt, exp, log, sin, cos, tan, asin, acos,
atan, expt, make-polar, magnitude, and angle may (but are not required
to) return inexact results even when given exact arguments.

> (log 1) => 0 ; an exact 0 not an approximated 0.0

In fact, this produces 0.0 in all of Scsh, Guile, and Chicken (both with
and without the numbers egg).

> (sin 0) => 0
> (cos 0) => 1 ; an exact 1
> (sqrt 4) => 2 ; an exact 2
> (sqrt -4) => 0+2i ; no approximations on the img part.
> (expt 16 1/2) => 4 
> (expt 8 1/3) => 2.0 ; uh oh? My Casio gives the exact value of 2 here.

And likewise with all of these.

> (/ 1 3) => 1/3 

This is indeed exact on Chicken with the numbers egg, Guile, and Scsh;
Chicken without the numbers egg prints a lexical-syntax warning, treats
1/2 as 0.5, and produces an inexact result.

-- 
John Cowan        http://ccil.org/~cowan   co...@cc...
Lope de Vega: "It wonders me I can speak at all.  Some caitiff rogue did
rudely yerk me on the knob, wherefrom my wits still wander."
An Englishman: "Ay, a filchman to the nab betimes 'll leave a man
crank for a spell." --Harry Turtledove, Ruled Britannia

Eddie Rucker wrote:
> What about ln 1? Scheme has a different philosophy from Pure: all
> operations *try* to return an exact value if the operation's argument(s)
> are exact.

That's true, of course. But in the case ln 0 there's no exact result. In 
fact, Scheme's approach to this raises many difficult questions, because 
it's not clear immediately what, in a precise mathematical sense, should 
be an exact result. That borders on computer algebra which is a *very* 
difficult area which goes way beyond just rational numbers and the like.

> I think I've totally
> chased the rabbit away from my original question which was about
> consistency for which you found the answer in a later post.

I agree, that should be the goal. And in fact you've been right, the 
definitions in math.pure have some bugs there.

There are many tradeoffs in numerical mathematics, even more so in 
computer algebra, design decisions are difficult, and sometimes I'm not 
even aware of them. So please continue to raise those questions, that 
always gives me an opportunity to rethink those things!

Thanks,
Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag

On Wed, 2008-09-10 at 00:43 +0200, Albert Graef wrote:
> There are many tradeoffs in numerical mathematics, even more so in 
> computer algebra, design decisions are difficult, and sometimes I'm not 

So true!

> even aware of them. So please continue to raise those questions, that 
> always gives me an opportunity to rethink those things!

OK, I'll keep fanning the poo pile. Keep your trusty flyswatter at
hand ;-)

e.r.

Albert Graef wrote:
> Good point! In fact gcc seems to do that, too (I get nan+:nan for /0.0 
> just as well as for /(0.0+I*0.0).

As a matter of fact, Haskell gives the same results, too (at least Hugs 
does), and the Haskell report uses the same formulas as math.pure in the 
complex/complex case.

I'm beginning to wonder whether it's really worth the hassle to follow 
the ISOC99 recommendations concerning the treatment of infinities and 
NaNs for complex * and /. The algorithms in Annex G do look like a kind 
of kludge to me, and they employ C99 functions for extracting and 
scaling the exponents of floating point numbers which might not be 
available on some systems (such as Windows), at least not directly. 
Hence special support in the runtime would be needed to implement these 
in a portable way across all platforms.

That seems like an awful lot of kludges to just make those nan+:nan's go 
away, considering that Pure also offers the polar representation where 
you'll get the correct infinite results for division by zero and similar 
cases, as long as the phase angles are finite.

So actually John's suggestion seems to be the most reasonable: just 
promote real operands of * and / to complex, then the results will at 
least be consistent. Or does anyone here really need the ISOC 
recommended behaviour?

Albert

-- 
Dr. Albert Gr"af
Dept. of Music-Informatics, University of Mainz, Germany
Email:  Dr....@t-..., ag...@mu...
WWW:    http://www.musikinformatik.uni-mainz.de/ag