Re: [CLN-list] [GiNaC-devel] evalf() on cube roots
Hi, This is a CLN topic. Let's take this thread to cln-list. Your proposal seems to suggest that CLN pick another branch cut for roots if a certain global flag is set. I'm sceptical about introducing global flags which control the computation in a way that results in very different numerical results - it looks like calling for problems. Admitted, there are global flags (cl_inhibit_floating_point_underflow and those dealing with I/O), and rounding control is well known from IEEE754. But returning a complex number which is far away from the conventional result in the complex plane is a different thing. I'm not aware of any such flag in another system; is there one? I'ld like to hear what others think. -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/> On 31.10.18 18:02, Jan Rheinländer wrote:
Picking up this old thread...
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
What about introducing a flag called (e.g.) cln::cl_odd_roots_of_rationals_real (similar to cln::cl_inhibit_floating_point_underflow)?
If the flag is false (default), the above case evaluates to exp(log(x) * y) as before
If the flag is true, it evaluates to (-1)^m * (-x)^(m/n)
Would such a patch be acceptable?
Greetings,
Jan Rheinländer
Am 02.08.2015 um 21:55 schrieb Richard B. Kreckel:
On 08/02/2015 08:59 PM, Jan Rheinländer wrote:
There isn't. GiNaC, like almost all other systems, returns the solution on the principal branch, compatible with exp(log(-1)/3). I suppose there is some excellent mathematical reason for this... Actually, there isn't. Branch cuts are a matter of convention. Nothing more but also nothing less.
but for me it means that I can't use GiNaC to verify Cardano's formula for a cubic function:
x = (-1 - sqrt(2))^(1/3) + (-1 + sqrt(2))^(1/3); (-1-sqrt(2))^(1/3)+(-1+sqrt(2))^(1/3)
evalf(x^3 + 3 * x + 2); 3.3544445609126528848+8.9073474964875349776*I
Surely there must be a way to verify such a formula in GiNaC??? It should be straightforward to replace such rational expressions by the form you want them to be in before calling evalf(), using techniques such as these: <http://www.ginac.de/FAQ.html#advanced>
Cheers -richy.
_______________________________________________ GiNaC-devel mailing list GiNaC-devel@ginac.de https://www.cebix.net/mailman/listinfo/ginac-devel
Hi, 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! I am using GiNaC for programming a tool that I use in engineering. Am 21.09.19 um 18:00 schrieb Richard B. Kreckel:
Hi,
This is a CLN topic. Let's take this thread to cln-list.
Your proposal seems to suggest that CLN pick another branch cut for roots if a certain global flag is set.
I'm sceptical about introducing global flags which control the computation in a way that results in very different numerical results - it looks like calling for problems. Admitted, there are global flags (cl_inhibit_floating_point_underflow and those dealing with I/O), and rounding control is well known from IEEE754. But returning a complex number which is far away from the conventional result in the complex plane is a different thing. I'm not aware of any such flag in another system; is there one?
I'ld like to hear what others think.
-richy.
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... [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
participants (3)
-
Bruno Haible
-
Jan Rheinländer
-
Richard B. Kreckel