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

Martin Maechler m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Fri Sep 6 10:37:36 CEST 2024


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