55

I've been finding a strange behaviour of log functions in C++ and numpy about the behaviour of log function handling complex infinite numbers. Specifically, log(inf + inf * 1j) equals (inf + 0.785398j) when I expect it to be (inf + nan * 1j).

When taking the log of a complex number, the real part is the log of the absolute value of the input and the imaginary part is the phase of the input. Returning 0.785398 as the imaginary part of log(inf + inf * 1j) means it assumes the infs in the real and the imaginary part have the same length. This assumption does not seem to be consistent with other calculation, for example, inf - inf == nan, inf / inf == nan which assumes 2 infs do not necessarily have the same values.

Why is the assumption for log(inf + inf * 1j) different?

Reproducing C++ code:

#include <complex>
#include <limits>
#include <iostream>
int main() {
    double inf = std::numeric_limits<double>::infinity();
    std::complex<double> b(inf, inf);
    std::complex<double> c = std::log(b);
    std::cout << c << "\n";
}

Reproducing Python code (numpy):

import numpy as np

a = complex(float('inf'), float('inf'))
print(np.log(a))

EDIT: Thank you for everyone who's involved in the discussion about the historical reason and the mathematical reason. All of you turn this naive question into a really interesting discussion. The provided answers are all of high quality and I wish I can accept more than 1 answers. However, I've decided to accept @simon's answer as it explains in more detail the mathematical reason and provided a link to the document explaining the logic (although I can't fully understand it).

Firman
  • 928
  • 7
  • 15
  • 2
    [python](https://ideone.com/IePnI7) and [cpp](https://godbolt.org/z/KnM83r5MM) – Marek R Dec 14 '22 at 12:58
  • 8
    [C++] _"...Error handling and special values..."_ _"...If z is (+∞,+∞), the result is (+∞,π/4).."_ https://en.cppreference.com/w/cpp/numeric/complex/log – Richard Critten Dec 14 '22 at 13:02
  • 2
    @RichardCritten that should be posted as an answer. The reason is probably "_The semantics of this function are intended to be consistent with the C function `clog`_". – Andras Deak -- Слава Україні Dec 14 '22 at 13:03
  • 1
    It would be interesting how far back in history this behaviour goes. [At least gfortran](https://tio.run/##S8svKilKzNNNT4Mw/v8vKMrMK1HXUPJIzcnJ11EIzy/KSVFU0lTnSs1L@f8fAA) also does this, though this may be more because it shares code with the C-focused GCC than because of Fortran itself. – leftaroundabout Dec 15 '22 at 10:09

3 Answers3

51

The free final draft of the C99 specification says on page 491

clog(+∞, +i∞) returns +∞ + iπ/4.

This is still the case currently. The C++ specification explains the same rules with the note

The semantics of this function are intended to be consistent with the C function clog.

I agree the behaviour is confusing from a math point of view, and arguably inconsistent with other inf semantics, as you pointed out. But pragmatically, it's part of the C standard, which makes it part of the C++ standard, and since NumPy normally relies on C behaviour (even in confusing cases), this is inherited in the Python example.

The standard-library cmath.log() function has the same behaviour (if you test it right...):

>>> import cmath

>>> cmath.log(complex(float('inf'), float('inf')))
(inf+0.7853981633974483j)

I have no means to investigate the rationale that went into the C standard. I assume there were pragmatic choices being made here, potentially when considering how these complex functions interact with each other.

Andras Deak -- Слава Україні
  • 33,737
  • 11
  • 83
  • 111
  • 1
    It's not that numpy relies on C behaviour, it's because most of Python libraries are developed in C (for performances purposes). – Fareanor Dec 14 '22 at 13:24
  • 2
    @Fareanor there's quite some baggage from C directly in NumPy, and I would expect huge pushback if someone proposed to "fix" this. The most egregious example is [how `np.int_` is platform-dependent](https://github.com/numpy/numpy/issues/9464#issuecomment-318185121) because the C `long` type is platform-dependent. Of course it's hard to disentangle "C does this" from "we've been doing it this way for 15 years". – Andras Deak -- Слава Україні Dec 14 '22 at 13:27
  • 1
    Yes but my point is that it's not a _"C does this ! So do we !"_. The real reason why Python libs are so dependent of C is that because they are developed in C under the hood. At the beginning, they weren't, and Python was too slow a language to be worthy to be used as it is today, so they decided to rewrite most of the Python libs in C (hence why the python implementation have that many shared objects embedded). My knowledge of Python is relatively poor, but I'm pretty sure it's the reason why there are such bounds between Python and C. – Fareanor Dec 14 '22 at 13:38
  • 2
    @Fareanor Of course, the premise is that all these extensions are in C. But this is an implementation detail, and libraries could hide weird choices made by C even though they are implemented in it. For example by special-casing `log` for complex numbers when the arguments are infinite. But point taken, and this is probably mostly a philosophical question. Pragmatically if you write your library in C you naturally do what it does. – Andras Deak -- Слава Україні Dec 14 '22 at 13:41
  • 4
    You're right. My guess is that they probably just forwarded the call to the c log function without doing anything more (that would explain that the results are the same). – Fareanor Dec 14 '22 at 13:43
  • 11
    @Fareanor: "it's not a "C does this ! So do we !"." As the person who implemented this, I can tell you that it is _exactly_ that. :-) The CPython cmath module doesn't wrap the C complex math operations; it implements everything directly. Special-case handling was deliberately chosen (by me) to match what's in the C99 standard. – Mark Dickinson Dec 15 '22 at 11:36
  • @MarkDickinson Oh right, my guess was wrong then. Thank you for the clarification. – Fareanor Dec 15 '22 at 11:46
  • 1
    @Fareanor: It wasn't a bad guess - if we could rely on some level of consistency from the various platform C libraries then it would indeed be better to just be wrapping those. But unfortunately reality intrudes, and it turns out that we can't rely on that consistency. – Mark Dickinson Dec 15 '22 at 11:51
  • @MarkDickinson Yes, that makes sense. Thank you for having taken the time to explain how and why things are done this way. – Fareanor Dec 15 '22 at 11:54
  • 2
    @MarkDickinson that's awesome insight, thanks. You should consider posting an answer, because even though it only addresses the Python part of the question, it's the least guesswork/dunno answer (and a very interesting part of relevant history). – Andras Deak -- Слава Україні Dec 15 '22 at 19:50
39

The value of 0.785398 (actually pi/4) is consistent with at least some other functions: as you said, the imaginary part of the logarithm of a complex number is identical with the phase angle of the number. This can be reformulated to a question of its own: what is the phase angle of inf + j * inf?

We can calculate the phase angle of a complex number z by atan2(Im(z), Re(z)). With the given number, this boils down to calculating atan2(inf, inf), which is also 0.785398 (or pi/4), both for Numpy and C/C++. So now a similar question could be asked: why is atan2(inf, inf) == 0.785398?

I do not have an answer to the latter (except for "the C/C++ specifications say so", as others already answered), I only have a guess: as atan2(y, x) == atan(y / x) for x > 0, probably someone made the decision in this context to not interpret inf / inf as "undefined" but instead as "a very large number divided by the same very large number". The result of this ratio would be 1, and atan(1) == pi/4 by the mathematical definition of atan.

Probably this is not a satisfying answer, but at least I could hopefully show that the log definition in the given edge case is not completely inconsistent with similar edge cases of related function definitions.

Edit: As I said, consistent with some other functions: it is also consistent with np.angle(complex(np.inf, np.inf)) == 0.785398, for example.

Edit 2: Looking at the source code of an actual atan2 implementation brought up the following code comment:

note that the non obvious cases are y and x both infinite or both zero. for more information, see Branch Cuts for Complex Elementary Functions, or Much Ado About Nothing's Sign Bit, by W. Kahan

I dug up the referenced document, you can find a copy here. In Chapter 8 of this reference, called "Complex zeros and infinities", William Kahan (who is both mathematician and computer scientist and, according to Wikipedia, the "Father of Floating Point") covers the zero and infinity edge cases of complex numbers and arrives at pi/4 for feeding inf + j * inf into the arg function (arg being the function that calculates the phase angle of a complex number, just like np.angle above). You will find this result on page 17 in the linked PDF. I am not mathematician enough for being able to summarize Kahan's rationale (which is to say: I don't really understand it), but maybe someone else can.

simon
  • 1,503
  • 8
  • 16
  • 3
    Thanks for this. It seems like if it's angle related, `inf` in real and in imaginary are assumed to be the equal length. – Firman Dec 14 '22 at 14:59
  • 1
    @Firman That's exactly what I tried to say with my overly convolved answer :D – simon Dec 14 '22 at 15:00
  • 10
    For a reason why `atan2(inf, inf) = pi/4` might be defined this way: it gives the correct quadrant, which is known even though the exact angle is unknown. Returning `nan` would give less information. – jpa Dec 15 '22 at 11:54
  • 1
    @Firman I added a link to a document that provides more insight. See my 2nd edit. – simon Dec 15 '22 at 14:34
  • 2
    @jpa you could make the same argument for `inf/inf`. That's still _clearly non-negative_, so defining it as +1 would give more information than NaN. But evidently they didn't find this convincing in that occasion. – leftaroundabout Dec 15 '22 at 18:34
  • @simon: From what I get from Kahan's book, there are 9 complex infinities: (inf ± j*k), (inf + j*inf), (inf - j*inf), (-inf ± j*k), (-inf + j*inf), (-inf - j*inf), (±k + j*inf), (±k - j*inf) and (nan plus inf). Each of those, when applying a math function, such as `log` deserves its own asymptotic value, and that extends also to `atan2` because `atan2` and complex `log` are linked. So for example in python `atan2(k, inf) = 0`, `atan2(inf, k) = pi/2`, `atan2(inf, inf) = pi/4`, and so on... [He](https://en.wikipedia.org/wiki/William_Kahan) seems to know what he is talking about. – rodrigo Dec 15 '22 at 20:44
  • Good find on that book chapter, simon. @rodrigo what isn't clear to me is why it's not an option to choose nan-valued complexes. If we want to assign a "real" number, it might as well be pi/4. The question is why `log` _and_ `atan2` don't return nan for doubly infinite cases. But admittedly I didn't read too much context around the final decision in the book. – Andras Deak -- Слава Україні Dec 15 '22 at 21:41
  • @AndrasDeak--СлаваУкраїні: I don't pretend to understand it fully, but I think it has to do with asymptotic behavior. For example `0 / inf = 0`, that is no surprise. Similarly, it is quite natural that `atan2(0, inf) = 0` instead of NaN as it is an horizontal asymptote; similarly `atan2(inf, 0) = pi/2` as it is a vertical asymptote. Thus it is not so unexpected that `atan2(inf, inf) = pi/4`, just in the middle of the other two, a 45º asymptote. – rodrigo Dec 15 '22 at 22:33
  • 1
    No, I think that's the crux of the issue, @rodrigo. When you have _two_ infinities you can't normally assume that they are "the same size". Infinities don't exist as such among real and complex numbers, they are just a limit. But then when you have `atan2(inf, inf)` that's actually `lim_{t->inf} atan2(f(t), g(t))` or something similar, and depending on `f` and `g` this could be anything. This is the exact reason why `inf/inf` and `0*inf` are `nan`. I wonder if there's a good explanation based on projective planes (where infinity exists, but what I'm unfamiliar with)... – Andras Deak -- Слава Україні Dec 15 '22 at 22:37
  • 1
    @AndrasDeak--СлаваУкраїні: Yes, I admit that the `inf/inf` vs `atan2(inf, inf)` looks weird. I think Kahan's book reasons in the Riemann sphere. Let me quote it: _Thus, any coherent scheme for computing complex products, quotients and logarithms at zero and infinity can be regarded as a scheme that rounds `arg(z)` into one of the ten values above [fractions of pi] when `z` is zero or infinite._ But note that this in general affects only the `phase` of a complex number, the magnitude will still be infinity for every of the 8 non-NaN infinity complex numbers... – rodrigo Dec 15 '22 at 23:22
  • 1
    ... The `log` complex function is funny because the imaginary part of the result is the phase of the argument. Thus the imaginary part of the `log` of any of those 8 non-NaN complex infinities should be a multiple of `pi/4`. Similarly `atan2(y, x)` can be defined as `phase(x + y*j)`, so if `x` or `y` are `inf` (but not NaN), you will have a multiple of `pi/4`, too. – rodrigo Dec 15 '22 at 23:25
  • @AndrasDeak--СлаваУкраїні and rodrigo: Thank you both for your feedback. Yes, indeed, Kahan being both a mathematician and computer scientist seems to know what he's talking about – only I did not understand him :D It seems like, you rodrigo, did better there. Maybe the mathematical reasoning and reference could be included into your answer, AndrasDeak? Then it would be both a technical one (referring to the C spec) and a mathematical one (referring to Kahan's reasoning). What do you think? – simon Dec 16 '22 at 09:06
  • @rodrigo I added more detail on Kahan to the answer. Thanks a lot for the hint! – simon Dec 16 '22 at 09:37
  • Thanks simon, but my answer doesn't try to explain the original reason, just the heritage. Your investigation's results should be in your answer :) – Andras Deak -- Слава Україні Dec 16 '22 at 09:45
  • @AndrasDeak--СлаваУкраїні Thank you as well, and I understand you, but in a way I was hoping for a "complete" answer to the "why" of the question. I think it can be understood in multiple ways – of which your answer gives the heritage aspect and mine (to some extent) points to the mathematical rationale. I feel an answer that is ready for acceptance would cover both, and I feel yours would be the better basis for that. But I don't have strong opinions there – in any case, the readers already can find both aspects covered as-is, if they take the time :) – simon Dec 16 '22 at 09:51
  • Yeah, having multiple relevant answers is by design of the platform. And I think your answer is a lot better answer to the "why" the question actually asks about. – Andras Deak -- Слава Україні Dec 16 '22 at 10:48
5

If we consider this from a pure maths point of view then we can look at operations in terms of limits e.g. as x goes to infinity, 1/x tends to 0 (denoted lim(x => inf) 1/x = 0), which is what we observe with floating points.

For operations on 2 infinities, we consider each infinity separately. Thus:

lim(x => inf) x/1 = inf
lim(x => inf) 1/x = 0

and in general we say inf/x = inf, and x/inf = 0. Thus:

lim(x => inf) inf/x = inf
lim(x => inf) x/inf = 0

Which of these 2 should we prefer? The spec for floats sidesteps by declaring it a NaN.

For complex logs however we observe:

lim(x=>inf) log(x + 0j) = inf + 0j
lim(x=>inf) log(0 + xj) = inf + pi/2j
lim(x=>inf) log(inf + xj) = inf + 0j
lim(x=>inf) log(x + infj) = inf + pi/2j

There is still a contradiction, but rather than being between 0 and inf, it is between 0 and pi/2, so the spec authors chose to split the difference. Why they made this chooice I couldn't say, but floating point infinities are not mathematical infinities, but instead represent 'this number is too large for representation'. Given that uses of log(complex) are likely to be more pure mathematical than those for subtraction and division, the authors may have felt that retaining the identity im(log(x+xj)) == pi/4 was useful.

rhellen
  • 159
  • 1
  • 2
    Worse, if `x` and `y` are separate variables then as `(x, y) -> (inf, inf)`, you can get any angle you like for `x + i*y` in the range `[0, pi/2]`. – kaya3 Dec 16 '22 at 16:53