[R] BUG: atan(1i) / 5 = NaN+Infi ?
Richard O'Keefe
r@oknz @end|ng |rom gm@||@com
Fri Sep 6 16:40:29 CEST 2024
G.5.1 para 2 can be found in the C17 standard -- I actually have the
final draft not the published standard. It's in earlier standards, I
just didn't check earlier standards. Complex arithmetic was not in
the first C standard (C89) but was in C99.
The complex numbers do indeed form a field, and Z*W invokes an
operation in that field when Z and W are both complex numbers. Z*R
and R*Z, where R is real-but-not-complex, is not that field operation;
it's the scalar multiplication from the vector spec view.
One way to characterise the C and Ada view is that real numbers x can
be viewed as (x,ZERO!!!) and imaginary numbers y*i can be viewed as
(ZERO!!!, y) where ZERO!!! is a real serious HARD zero, not an IEEE
+epsilon or -epsilon, In fact this is important for getting "sign of
zero" right.
x + (u,v) = (x+u, v) *NOT* (x+u, v+0) and this does matter in IEEE arithmetic.
R is of course based on S, and S was not only designed before C got
complex numbers, but before there was an IEEE standard with things
like NaN, Inf, and signed zeros But there *is* an IEEE standard now,
and even IBM mainframes offer IEEE-compliant arithmetic, so worrying
about the sign of zero etc is not something we can really overlook
these days.
You are of course correct that the one-point compactification of the
complex numbers involves adjoining just one infinity and that whacking
IEEE infinities into complex numbers does not give you anything
mathematically interesting (unless you count grief and vexation as
"interesting" (:-)). Since R distinguishes between 0+Infi and
NaN+Infi, it's not clear that the one-point compactification has any
relevance to what R does. And it's not just an unexpected NaN; with
all numbers finite you can get zeros with the wrong sign. (S having
been designed before the sign of zero was something you needed to be
aware of.)
For what it's worth, the ISO "Language Independent Arithmetic"
standard, part 3, defines separate real, imaginary, and complex types,
and defines x*(z,w) to be (x*z, x*w) directly, just like C and Ada.
So it is quite clear that R does not currently conform to LIA-3.
LIA-3 (ISO/IEC 10967-3:2006) is the nearest we have to a definition of
what "right" answers are for floating-point complex arithmetic, and
what R does cannot be called "right" by that definition. But of
course R doesn't claim conformance to any part of the LIA standard.
Whether the R community *want* R and C to give the same answers is not
for me to say. I can only say that *I* found it reassuring that C
gave the expected answers when R did not, or, to put it another way,
disconcerting that R did not agree with C (or LIA-3).
What really annoys me is that I wrote an entire technical report on
(some of the) problems with complex arithmetic, and this whole "just
treat x as (x, +0.0)" thing completely failed to occur to me as
something anyone might do.
On Fri, 6 Sept 2024 at 20:37, Martin Maechler
<maechler using stat.math.ethz.ch> wrote:
>
> >>>>> Richard O'Keefe
> >>>>> on Fri, 6 Sep 2024 17:24:07 +1200 writes:
>
> > The thing is that real*complex, complex*real, and complex/real are not
> > "complex arithmetic" in the requisite sense.
>
> > The complex numbers are a vector space over the reals,
>
> Yes, but they _also_ are field (and as others have argued mathematically only
> have one infinity point),
> and I think here we are fighting with which definition should
> take precedence here.
> The English Wikipedia page is even more extensive and precise,
> https://en.wikipedia.org/wiki/Complex_number (line breaking by me):
>
> " The complex numbers form a rich structure that is simultaneously
> - an algebraically closed field,
> - a commutative algebra over the reals, and
> - a Euclidean vector space of dimension two."
>
> our problem "of course" is that we additionally add +/- Inf for
> the reals and for storage etc treating them as a 2D vector space
> over the reals is "obvious".
>
> > and complex*real and real*complex are vector*scalar and scalar*vector.
> > For example, in the Ada programming language, we have
> > function "*" (Left, Right : Complex) return Complex;
> > function "*" (Left : Complex; Right : Real'Base) return Complex;
> > function "*" (Left : Real'Base; Right : Complex) return Complex;
> > showing that Z*R and Z*W involve *different* functions.
>
> > It's worth noting that complex*real and real*complex just require two
> > real multiplications,
> > no other arithmetic operations, while complex*complex requires four
> > real multiplications,
> > an addition, and a subtraction. So implementing complex*real by
> > conventing the real
> > to complex is inefficient (as well as getting the finer points of IEEE
> > arithmetic wrong).
>
> I see your point.
>
> > As for complex division, getting that *right* in floating-point is
> > fiendishly difficult (there are
> > lots of algorithms out there and the majority of them have serious flaws)
> > and woefully costly.
>
> > It's not unfair to characterise implementing complex/real
> > by conversion to complex and doing complex/complex as a
> > beginner's bungle.
>
> ouch! ... but still I tend to acknowledge your point, incl the "not unfair" ..
>
> > There are good reasons why "double", "_Imaginary double", and "_Complex double"
> > are distinct types in standard C (as they are in Ada),
>
> interesting. OTOH, I think standard C did not have strict
> standards about complex number storage etc in the mid 1990s
> when R was created.
>
> > and the definition of multiplication
> > in G.5.1 para 2 is *direct* (not via complex*complex).
>
> I see (did not know about) -- where can we find 'G.5.1 para 2'
>
> > Now R has its own way of doing things, and if the judgement of the R
> > maintainers is
> > that keeping the "convert to a common type and then operate" model is
> > more important
> > than getting good answers, well, it's THEIR language, not mine.
>
> Well, it should also be the R community's language,
> where we, the R core team, do most of the "base" work and also
> emphasize guaranteeing long term stability.
>
> Personally, I think that
> "convert to a common type and then operate"
> is a good rule and principle in many, even most places and cases,
> but I hate it if humans should not be allowed to break good
> rules for even better reasons (but should rather behave like algorithms ..).
>
> This may well be a very good example of re-considering.
> As mentioned above, e.g., I was not aware of the C language standard
> being so specific here and different than what we've been doing
> in R.
>
>
> > But let's not pretend
> > that the answers are *right* in any other sense.
>
> I think that's too strong -- Jeff's computation (just here below)
> is showing one well defined sense of "right" I'd say.
> (Still I know and agree the Inf * 0 |--> NaN
> rule *is* sometimes undesirable)
>
> > On Fri, 6 Sept 2024 at 11:07, Jeff Newmiller via R-help
> > <r-help using r-project.org> wrote:
> >>
> >> atan(1i) -> 0 + Inf i
> >> complex(1/5) -> 0.2 + 0i
> >> atan(1i) -> (0 + Inf i) * (0.2 + 0i)
> -> 0*0.2 + 0*0i + Inf i * 0.2 + Inf i * 0i
> >> infinity times zero is undefined
> -> 0 + 0i + Inf i + NaN * i^2
> -> 0 + 0i + Inf i - NaN
> -> NaN + Inf i
> >>
> >> I am not sure how complex arithmetic could arrive at another answer.
> >>
> >> I advise against messing with infinities... use atan2() if you don't actually need complex arithmetic.
> >>
> >> On September 5, 2024 3:38:33 PM PDT, Bert Gunter <bgunter.4567 using gmail.com> wrote:
> >> >> complex(real = 0, imaginary = Inf)
> >> >[1] 0+Infi
> >> >
> >> >> Inf*1i
> >> >[1] NaN+Infi
> >> >
> >> >>> complex(real = 0, imaginary = Inf)/5
> >> >[1] NaN+Infi
> >> >
> >> >See the Note in ?complex for the explanation, I think. Duncan can correct
> >> >if I'm wrong.
> >> >
> >> >-- Bert
>
> [...................]
>
> Martin
>
> --
> Martin Maechler
> ETH Zurich and R Core team
More information about the R-help
mailing list