[R] BUG: atan(1i) / 5 = NaN+Infi ?

[Martin Maechler]

>>>>> Richard O'Keefe 
>>>>>     on Sat, 7 Sep 2024 02:40:29 +1200 writes:

    > G.5.1 para 2 can be found in the C17 standard -- I
    > actually have the final draft not the published standard.

Ok.  Thank you.

A direct hopefully stable link to that final draft's  Appendix G
seems to be    
  https://www.open-std.org/JTC1/SC22/WG14/www/docs/n2310.pdf#chapter.14

which is good to have available.

    > 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.

indeed, currently we only require C11 for R.

A longer term solution that we (the R core team) will probably
look into is to start making use of current C standards for
complex arithmetic.

As mentioned, all this was not yet available when R started and
already came with complex numbers as base type  ....
It may (or may not, I'm not the expert) be a bit challenging
trying to remain back compatible (e.g. with save complex number
R objects)  and still use C standard complex headers ...

But mid to long term I guess that would be the way to go.

Martin

    > 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).


Thank you.
Indeed, I'd  like the idea to consider  LIA  standards as much
as (sensibly) possible.

Martin


    > 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