Hmm, I forgot to add explicit test with zero for the mp context... Also note using mpf_ge() to correctly handle unsigned zeros.

But we have in Mathematica:

In[7]:= ArcSin[Infinity + I]
Out[7]= -I Infinity

In[8]:= ArcSin[Infinity]
Out[8]= -I Infinity

In[9]:= ArcSin[10.^10 + I]
Out[9]= 1.5708 + 23.719 I

Ditto Diofant/SymPy:

In [11]: asin(oo)
Out[11]: -∞⋅ⅈ

In [12]: asin(1e+10)
Out[12]: 1.5707963267949 - 23.7189981105004⋅ⅈ

So, it's a compatibility break:
https://github.com/sympy/sympy/blob/a9a188ab40f260d83df5790a96cfac9d6d34d8ec/sympy/functions/elementary/trigonometric.py#L2178-L2181

In this case, I think, it's worth to do. I don't see good reasons why real component was zeroed c.f. case of finite big argument: Mathematica just discards finite components of a complex number with presence of infinity, so it looks consistent for this CAS, but not for diofant/sympy. @oscarbenjamin ?

Well, we have a disgrepancy for diofant/sympy evaluation rules and mpmath already in v1.3.0:

In [1]: asin(oo, evaluate=False).evalf()  # in sympy 1.13.dev
Out[1]: 1.5707963267949 - ∞⋅ⅈ

In [2]: import mpmath;mpmath.__version__
Out[2]: '1.3.0'

In [3]: asin(oo)
Out[3]: -∞⋅ⅈ

Probably this affects a lot of functions (but in many cases it might be a mpmath fault, see e.g. #776). I believe, this example is a bug on CAS's side. Following invariant should be preserved: f(*args, evaluate=False).evalf() == f(*args, evaluate=True).evalf(). What do you think?

According to https://en.cppreference.com/w/cpp/numeric/complex/asin special cases,

asin(complex(+oo, +x)) = -I * asinh(complex(-x, +oo))
  = -I * conj(asinh(complex(-x, -oo)))
  = -I * conj(-asinh(complex(+x, +oo)))
  = -I * conj(-complex(+oo, pi/2))
  = I * complex(+oo, -pi/2) = complex(pi/2, +oo)

asin(complex(+oo, -x)) = -I * asinh(complex(+x, +oo))
  = -I * complex(+oo, pi/2) = complex(pi/2, -oo)

for any positive finite x.

While the special cases do not specify the case where x is 0 but assuming identity 0 == +0, asin(complex(+oo, +0)) = complex(pi/2, +oo).

Notice that the sign difference on the branch cut that has been discussed in #786 which also demonstrates the systematic inconsistency between FP and symbolic software.

Sympy asin(oo) -> -oo * I is likely a bug or a consequence of some internal manipulation where in expressions containing infinities, all (presumably) finite symbols/terms are discarded while ignoring the fact that terms may have different unit-vectors (of real and imaginary axes).

According to https://en.cppreference.com/w/cpp/numeric/complex/asin special cases

I think, these are taken into account correctly.

While the special cases do not specify the case where x is 0

In fact, they do, same rule "If z is x+∞i (for any positive finite x), the result is +∞+π/2" + conjugation for asinh // https://en.cppreference.com/w/c/numeric/complex/casinh

You are correct, for zero you just have to account it's sign.

Notice that the sign difference on the branch cut that has been discussed in #786 which also demonstrates the systematic inconsistency between FP and symbolic software.

It's not an inconsistency. FP arithmetic with signed zero just can't represent the value directly on the branch cut. I'm trying to document that. CCC rule does work for elementary functions, but it's tricky to specify on which side we choose to be continuous in general case...

Sympy asin(oo) -> -oo * I is likely a bug or a consequence of some internal manipulation where in expressions containing infinities

I would guess it's a blindly copy from Mathematica, where e.g.

In[1]:= Infinity + I
Out[1]= Infinity

In[2]:= Pi + I Infinity
Out[2]= I Infinity

Here above trick does make sense; Mathematica has only concept of DirectedInfinity/ComplexInfinity, but in the Diofant/SymPy:

In [7]: oo + 1
Out[7]: ∞

In [8]: I*oo + 1
Out[8]: 1 + ∞⋅ⅈ

Draft PR for the Diofant to fix evaluation rules: diofant/diofant#1396

I think, that the SymPy should implement something like that. Current evaluation rules could be broken also by rewrite rules, for example:

In [1]: asin(oo)
Out[1]: -∞⋅ⅈ

In [2]: asin(oo, evaluate=False).rewrite(acos)
Out[2]: 
π      
─ - ∞⋅ⅈ
2      

Mathematica has only concept of DirectedInfinity/ComplexInfinity, but in the Diofant/SymPy:

In [7]: oo + 1
Out[7]: ∞

In [8]: I*oo + 1
Out[8]: 1 + ∞⋅ⅈ

I have contemplated whether this should be changed so that there is only directed infinity and complex infinity. I'm not sure that 1 + oo*I is really meaningful or useful.

I have contemplated whether this should be changed so that there is only directed infinity and complex infinity.

Yes, maybe that should be changed that way.

Pros: this will match current limit output. For example, in the mentioned Diofant's pr:

In [1]: asin(oo)
Out[1]: 
π      
─ - ∞⋅ⅈ
2      

In [2]: asin(z).limit(z, oo)
Out[2]: -∞⋅ⅈ

This could be fixed, however. Or documented as a feature of the limit algorithm:)

I'm not sure that 1 + oo*I is really meaningful

Why not? It's just as meaningful as "oo*I" == "positive infinity on the imaginary line". This one - on the shifted line.

Of course, there is possible another meaning: "directed infinity with argument pi/2". But above interpretation seems to be useful and widespread notation, e.g. for integrals:
https://docs.sympy.org/latest/modules/integrals/g-functions.html#the-inverse-laplace-transform-of-a-g-function

or useful.

Sort of. Clearly, you lose something e.g. in the following example (current SymPy):

In [6]: asin(1.e+10)
Out[6]: 1.5707963267949 - 23.7189981105004⋅ⅈ

In [7]: asin(1.e+100)
Out[7]: 1.5707963267949 - 230.951656479965⋅ⅈ

In [8]: asin(oo)
Out[8]: -∞⋅ⅈ

And if the mp context will adopt that convention - there will be yet another difference wrt the fp context.

1 + oo*I in Mathematica is DirectedInfinity[1 + I] - that is wrong, sorry for the noise.

Within the rectangular form of complex numbers that mpmath.mpc and Python/Numpy/C++... complex data types are all about, the complex infinity (the top point of Riemann sphere) has no representation. While Sympy/Mathematica/any CAS may invent the corresponding symbol, in the context of mpmath it is basically OT as it uses the rectangular form.

Within the rectangular form of complex numbers that mpmath.mpc

I think that's an implementation detail. E.g. we could trivially change back and forth to the polar form in finite case.

the context of mpmath it is basically OT as it uses the rectangular form.

Infinite complex numbers per the c11 have 8 variants as directed infinities, e.g. complex('(inf+infj)') is a directed infinity with argument pi/4, complex('(inf+1j)') is a directed infinity with argument 0, etc... Perhaps, we could introduce special form for more generic notion of directed infinity in the mpmath. Not sure if it worth.

Infinite complex numbers per the c11 have 8 variants as directed infinities, ...

Absolutely. The OT remark was for the complex infinity as the top point of Riemann sphere that Mathematica ComplexInfinity basically represents.

Perhaps, we could introduce special form for more generic notion of directed infinity in the mpmath.

I think that is a topic for libraries that use mpmath, such as sympy, but for mpmath which is about multiple-precision arithmetic, it would be off-topic, IMHO.

But above interpretation seems to be useful and widespread notation, e.g. for integrals:

It is widespread notation to use c + I*oo in integral limits but I consider it an abuse of notation. That is fine but I don't think that it means we should try to make c + I*oo have a well defined meaning as an expression on its own. A well defined expression could be e.g. BromwichIntegral(expr, z) and could display using c + I*oo without any semantic ambiguity.

It is widespread notation to use c + I*oo in integral limits

As well as oo + I*c.

I consider it an abuse of notation

I see. Lets then put #793 on hold, I'll factor out other commits to a separate pr. I would appreciate if you could discuss this issue in the SymPy community. (Maybe you could mention also #786: the current decision clearly will not make sense if the mp context got signed zeros, see #167)

I doubt you convinced me so far. The notion of complex numbers with infinite components seems to be well defined, used in textbooks and in numerical libraries. Etc, see #795 (comment) I think that CAS's could support this notion together with directed infinity, which is already possible in the Diofant/SymPy:

In [6]: oo*exp(I*pi/6)
Out[6]: 
   ⅈ⋅π
   ───
    6 
∞⋅ℯ   

Regarding limits issue, you will have to fix something in the SymPy anyway, because:

In [2]: limit(asin(x), x, oo)
Out[2]: -∞⋅ⅈ

In [3]: limit(re(asin(x)), x, oo)
Out[3]: 
π
─
2

In [4]: limit(im(asin(x)), x, oo)
Out[4]: -∞

A well defined expression could be e.g. BromwichIntegral(expr, z)

What about integrals over (-oo + I*c, oo + I*c)?

could display using c + I*oo without any semantic ambiguity

I.e. "lets just use that for pretty-printing". Ok.

BTW, Mathematica's convention doesn't look consistent (i.e. there is more than just DirectedInfinity and zoo): addition/subtraction of DirectedInfinity'es left intact:

In[1]:= DirectedInfinity[Pi/2] + DirectedInfinity[Pi/3]

Out[1]= Infinity

In[2]:= DirectedInfinity[Exp[I Pi/2]] + DirectedInfinity[Exp[I Pi/3]]

                               I/3 Pi
Out[2]= I Infinity + Infinity E

I think that CAS's could support this notion together with directed infinity, which is already possible in the Diofant/SymPy:

I think that we might be misunderstanding each other when talking about directed infinities. To me a directed infinity is exp(I*theta)*oo for some real theta. I think that this notion of directed infinity and also zoo are the only infinities that should be needed (in Diofant/SymPy).

Maybe Mathematica's DirectedInfinity is something different and maybe also this (my) notion of directed infinity does not really work in mpmath...

Maybe Mathematica's DirectedInfinity is something different

No, it's exactly that one. And how you can see, Mathematica uses more infinity types. Not sure why pi+ I*oo is worse than I*oo + I*exp(I*pi/3)...

BTW:

>>> gmpy2.asin(cmath.inf)
mpc('1.5707963267948966+infj')
>>> flint.acb(cmath.inf).asin()  # I hope flint.ctx has no option like gmpy2's allow_complex?
nan + nanj

In fact, it seems that the algorithm is correct in the imaginary part (if the imaginary component of argument is not a nan). So, here is a simplified patch, that fix the real part at the end.

Maybe Mathematica's DirectedInfinity is something different

No, it's exactly that one. And how you can see, Mathematica uses more infinity types. Not sure why pi+ I*oo is worse than I*oo + I*exp(I*pi/3)...

It is just not clear what model of "infinity" is being used when we have something like pi + I*oo. In some sense oo is a shorthand for some sort of unbounded limit but I'm not sure that complex numbers with extended reals for real and imaginary parts is well defined.

In my notion of directed infinity I*oo + I*exp(I*pi/3) should be undefined or nan or something like oo + (-oo) is. I suppose leaving it unevaluated allows some operations to work like abs or possibly re/im.

With SymPy:

In [1]: oo + (-oo)
Out[1]: nan

In [2]: I*oo + oo
Out[2]: ∞ + ∞⋅ⅈ

In [3]: re(I*oo + oo)
Out[3]: ∞

In [4]: arg(I*oo + oo)
Out[4]: 
π
─
4

I don't like the arg(I*oo + oo) = pi/4 result. It could be nan or an error. Potentially AccumBounds(0, pi/2) is acceptable to make it somewhat well defined. It should not be a definite answer like pi/4 though.

It is just not clear what model of "infinity" is being used when we have something like pi + I*oo. In some sense oo is a shorthand for some sort of unbounded limit but I'm not sure that complex numbers with extended reals for real and imaginary parts is well defined.

The re(asin(z)) is just a real-valued function, isn't? Just as im(asin(z)). It's perfectly meaningful to consider component-wise limits, even if they are unbounded. In that sense pi + I*oo seems to be a useful notion and I don't see here something mathematically wrong. If the I*oo will eat finite real component per default evaluation rules - we lose some information.

I suppose leaving it unevaluated allows some operations to work like abs or possibly re/im.

Indeed. In somewhat a funny way...

In[1]:= I Infinity + I Exp[I Pi/3] Infinity

                                 I/3 Pi
Out[1]= I Infinity + I Infinity E

In[2]:= Abs[%1]

                                     I/3 Pi
Out[2]= Abs[I Infinity + I Infinity E      ]

In[3]:= Re[%1]

Out[3]= -Infinity

In[4]:= Im[%1]

Out[4]= Infinity

I would prefer to keep last two cases unevaluated too. On contrary, Abs may return Infinity.

I don't like the arg(I*oo + oo) = pi/4 result.

That's consistent with https://en.cppreference.com/w/c/numeric/complex/carg:

>>> math.atan2(math.inf, math.inf)  # as in the mpmath's master; v1.3.0 shows ``nan``.
0.7853981633974483 

(Diofant/SymPy says it's 0...) See also Kahan's "“Branch cuts for complex elementary functions" paper.

This looks as a special case for a directed infinity (in this case expansion/collection is legal, as well as for other similar 3 cases):

(1 + I)/sqrt(2)*oo == (1 + I)*oo = oo + I*oo

This looks as a special case for a directed infinity (in this case expansion/collection is legal, as well as for other similar 3 cases):

(1 + I)/sqrt(2)*oo == (1 + I)*oo = oo + I*oo

But we also have:

In [1]: expand((2+I)*oo)
Out[1]: ∞ + ∞⋅ⅈ

In [2]: arg(_)
Out[2]: 
π
─
4

This is inconsistent with my idea of directed infinity for which the argument is the defining feature i.e. arg((2+I)*oo) = arg(2+I).

The re(asin(z)) is just a real-valued function, isn't?

Yes, and the limit of the composite function (re(asin(z)).limit(z, oo)) is well defined. That does not necessarily mean that we should be able to evaluate the function at the limit in composite steps like a = asin(oo); r = re(a) though.

Evaluation with infinities should only work if all interpretations of the evaluation as a limit are equivalent which is why oo - oo needs to be undefined. I think that it is reasonable for many possible evaluation results to just be undefined for infinite arguments if that means that we can get to a place where all results are either well defined or clearly undefined.

But we also have

I didn't say such expansions are correct in all cases. This example is a bug. Seems to be an instance of sympy/sympy#4596, which is really solvable.

That does not necessarily mean that we should be able to evaluate the function at the limit in composite steps like a = asin(oo); r = re(a) though.

I doubt we have many options here. Either we could consider C as a metric space (Direction->Complexes in the Mathematica) and use that do define limits, as usual. Or we can consider limits of a function in a given direction on the complex plane (Direction->Exp[I theta] in the Mathematica), from the given origin (not just 0!), i.e. essentially as complex-valued functions of a real argument.

In the first case, we have only ComplexInfinity.

But in the second, we have all kind of DirectedInfinity's and also things like pi + I*oo. Just as for ±oo on the extended real line, it's easy to deduce right limit direction from such special infinity end point. E.g.:

• limit(f(z), z, I*oo) - on the imaginary axis, in direction pi/2 (this is something, that Diofant already has)
• limit(f(z), z, 1 + I*oo) - along the line, parallel to the imaginary axes, shifted by 1

But we also have

I didn't say such expansions are correct in all cases. This example is a bug. Seems to be an instance of sympy/sympy#4596, which is really solvable.

No it doesn't have anything to do with automatic distribution:

In [1]: e = (2+I)*oo

In [2]: e
Out[2]: ∞⋅(2 + ⅈ)

In [3]: arg(e)
Out[3]: atan(1/2)

In [4]: e.expand()
Out[4]: ∞ + ∞⋅ⅈ

In [5]: arg(e.expand())
Out[5]: 
π
─
4

If Out [4] is correct then Out [5] is not because we cannot conclude any specific argument for oo + I*oo if we have no idea what limit generated the two separate infinities e.g.:

In [7]: limit(arg(x+I*x**2), x, oo)
Out[7]: 
π
─
2

In [8]: limit(arg(x**2+I*x), x, oo)
Out[8]: 0

No it doesn't have anything to do with automatic distribution:

Oh, I see, that's fixed in the SymPy. I would definitely want to backport that for the Diofant.

If Out [4] is correct then Out [5] is not

IIUC, "mul" helper for expand just blindly do multiplication term-by-term. Again, this is a bug: in this case, SymPy shouldn't do it (there is kwarg force to enforce things like that).

if we have no idea what limit generated the two separate infinities e.g.

Hmm, why not?

limit(arg(x+I*x**2), x, oo) = limit(arg(x**2*(I + x**-1)), x, oo)
                            = limit(arg(I + x), x, 0-) = arg(I) = pi/2
limit(arg(x**2+I*x), x, oo) = limit(arg(x**2*(1 + I*x**-1)), x, oo)
                            = limit(arg(1 + I*x), x, 0-) = arg(1) = 0

In [1]: Limit(x + I*x**2, x, oo)
Out[1]: 
    ⎛   2    ⎞
lim ⎝ⅈ⋅x  + x⎠
x─→∞

In [2]: _.doit(heuristics=False)
Out[2]: ∞⋅ⅈ

In [3]: Limit(x**2 + I*x, x, oo)
Out[3]: 
    ⎛ 2      ⎞
lim ⎝x  + ⅈ⋅x⎠
x─→∞

In [4]: _.doit(heuristics=False)
Out[4]: ∞

No it doesn't have anything to do with automatic distribution:

Oh, I see, that's fixed in the SymPy.

No, it is not fixed. It just doesn't apply in this particular case because neither factor is a Rational.

if we have no idea what limit generated the two separate infinities e.g.

Hmm, why not?

The limits are all well defined but not the same. The problem is that if we don't have a limit but rather an expression with infinities like oo + I*oo then there are many possible limits limit(f(x)+I*g(x), x, oo) that the expression might represent. Those different possibilities have different limiting values for arg so arg(oo + I*oo) is not well defined. The different limits can give anything in [0, pi/2] so evaluating arg(oo + I*oo) = pi/4 is incorrect.

As I said above:

Evaluation with infinities should only work if all interpretations of the evaluation as a limit are equivalent

No, it is not fixed. It just doesn't apply in this particular case

Sure, I meant that. I realize, that sympy/sympy#4596 is not fixed.

(Though, it's interesting how many test failures are if distribution will be disabled for non-integer Rational's.)

there are many possible limits limit(f(x)+I*g(x), x, oo) that the expression might represent.

Just as with oo.

Those different possibilities have different limiting values

No. All such limits are equal to oo*exp(I*pi/4). In this particular case the exponent could be expanded over multiplication wrt oo.

Unfortunately, I was wrong in another point: such limits can't represent values like 1+I*oo. To get that, you have to understand limits component-wise, i.e.:

limit(f(z), z, z0) = limit(re(f(z)), z, z0) + I*limit(im(f(z)), z, z0)

So, returning to our special case, both

In [1]: limit(asin(z), z, oo)
Out[1]: -∞⋅ⅈ

In [2]: limit(re(asin(z)), z, oo)
Out[2]: 
π
─
2

have sense. But which one we will use to define asin(oo)?

Both kinds of infinities are useful: oo*exp(I*phi) have defined argument, pi + I*oo - provide information about real and imaginary component. But in general we can't have defined both arg(z) and re(z)/im(z) for infinities.

I think useful definition of asin(oo) should preserve maximum details about behaviour of asin for z->oo.

Ok, I'm going to merge this tomorrow: #793 (comment)

In short, it seems there is no compatibility break: we just extend the present (e.g. in v1.3) convention to other arguments with infinite components, besides asin(±oo); and we do this consistently, following existing standards for numerical libraries.