[CLN-list] [GiNaC-devel] evalf() on cube roots

[Bruno Haible]

Hi Jan,

> maybe I should explain why I am considering this. My background is in
> engineering (mechanical) and here complex numbers are considered as
> extremely exotic. For me as an engineer the idea that the cubic root of
> -1 should be complex sounds absurd. Most engineers probably have the
> opinion that -1 has only one cubic root anyway!

You are surprising me. Even in mechanical engineering, I would expect waves
to appear; and with waves you have to consider wave propagation, interference,
and resonance. Some of these topics are best modeled with Laplacians, which
imply complex numbers [1].

But anyway, yes I understand there are applications where you have real
numbers and want/expect real numbers as a result. This is why CLN has
the type cl_R and many functions that, given a cl_R as input, are guaranteed
to produce a cl_R result.

> > I would like to submit a patch that allows choosing the convention for
> > negative cube roots between the solution on the principal branch and the
> > real-valued solution.
> > 
> > That is, in the cln library, file complex/transcendental/cl_C_expt_C.cc
> > I would like to handle the following case separately:
> > 
> > y rational and y ratio m/n
> > x rational and x < 0
> > n odd

The choice which branch to use for a transcendental function is not a
decision that can/should be done for a particular x value without considering
the other x values.

Rather, one needs to consider the function as a whole and - where applicable -
related transcendental functions as well. [2]

The branch cuts implemented in CLN follow the specification in CLtL2 [3].
While this specification is arbitrary (e.g. why is the branch cut along the
negative real axis continuous with the upper half plane, not with the lower
half plane?) and in one place outright mistaken, it is nevertheless worth
implementing and following.

In particular, [3] and [4] specify that expt(x,y), for a fixed y, as a function
of x, has its branch cut along the negative real axis, continuous with
quadrant II (i.e. the upper half plane):

(expt #C(-8.0 -0.001) 1/3)   => #C(1.0000721 -1.7320092)
(expt #C(-8.0 0.0) 1/3)      => #C(1.0       1.7320508)
(expt #C(-8.0 0.001) 1/3)    => #C(1.0000721 1.7320092)

> > What about introducing a flag called (e.g.)
> > cln::cl_odd_roots_of_rationals_real

Aside from Richy's point that a flag is not a terribly good concept for
a library:
This approach would be changing expt as a cl_N -> cl_N function.
But what you actually want is a cl_R -> cl_R function.

So, I would suggest to instead introduce a new function
  root: (cl_R x, cl_I n) -> cl_R
that is defined as follows:
  - for x >= 0, it is the same as expt(x,1/n).
  - for x < 0 and n odd, it is the same as - root(-x, n).
  - for x < 0 and n even, it is undefined.

The latter case will require adding an invalid_argument_exception to the
hierarchy:
  runtime_exception
   +-- invalid_argument_exception
        +-- division_by_0_exception

This would be similar to the cbrt() function that many numerical libraries
have, and generalize it to n=5, n=7, etc. and negative exponents as well.

Bruno

[1] https://en.wikipedia.org/wiki/Laplace_transform#Examples_and_applications
[2] http://www.ai.mit.edu/projects/iiip/doc/CommonLISP/HyperSpec/Body/sec_12-1-5-4.html
[3] https://www.cs.cmu.edu/Groups/AI/html/cltl/clm/node129.html
[4] https://www.cs.cmu.edu/Groups/AI/html/cltl/clm/node127.html