How calculates Matlab the sqrt of a complex number?

Question

i am just analysing a C Implementation of my Algorithm vs the Matlab-Algorithm. It works quite fine, exceptionally when it Comes to calculate the square root of a complex number.

I tried already 3 different implementations on how to calculate a complex square root in C, but None of this implementation Matches the matlab result.

The Problem is, that the sign of the imaginary parts differs from my implementation. Here are two examples:

s = -1.1721375 - 0.0000000i
sqrt(s) = 0.0000000 - 1.0826530i

s = -1.7648506 - 0.0478944i
sqrt(s) = 0.0180244 - 1.3285991i

I tried to pick the principal square root, which is the one with the positive real part. This works for the first example, but not for the second case.

I Need to have the same result as calculated by matlab.

Thanks in Advance!

EDIT: I am trying to get all cases correct by the following Code

boolean sHavePositiveImaginaryPart = (s.imag > 0);
s = sqrt(s);
if (sHavePositiveImaginaryPart && s.imag > 0)
s = -1 * s;

Is this correct?

EDIT2: I have collected some test cases. It seems the sign of the imaginary part always remain unchanged:

s =
 -1.1721375 - 0.0000000i
squareRootOfs =
  0.0000000 - 1.0826530i

s =
 -1.7648506 - 0.0478944i
squareRootOfs =
  0.0180244 - 1.3285991i

s =
 -0.0193461 + 0.3215795i
squareRootOfs =
  0.3891110 + 0.4132233i

s =
  1.9100717e+02 - 7.1220589e+01i
squareRootOfs =
 14.0509853 - 2.5343628i

s =
 -1.7817555e+08 - 1.0207493e+08i
squareRootOfs =
  3.6856221e+03 - 1.3847721e+04i

s =
 -6.5224486e+19 + 1.5745077e+20i
squareRootOfs =
  7.2526336e+09 + 1.0854731e+10i

s =
  9.1788721e+23 + 6.5826579e+23i
squareRootOfs =
  1.0117840e+12 + 3.2529957e+11i

s =
  1.9785707e+26 + 1.0730929e+26i
squareRootOfs =
  1.4542022e+13 + 3.6896273e+12i

s =
  1.2809769e+28 - 7.7697397e+27i
squareRootOfs =
  1.1788071e+14 - 3.2955940e+13i

s =
 -5.8533863e+04 + 1.9665835e+03i
squareRootOfs =
  4.0636616e+00 + 2.4197185e+02i

s =
  5.5913438e+12 + 1.7300604e+13i
squareRootOfs =
  3.4476833e+06 + 2.5090188e+06i

s =
  1.3380493e+27 + 1.4976028e+27i
squareRootOfs =
  4.0904337e+13 + 1.8306161e+13i

s =
 -1.5978770 - 2.3017066i
squareRootOfs =
  0.7759182 - 1.4832147i


s =
 -1.5974181e+14 - 2.1267702e+14i
squareRootOfs =
  7.2885170e+06 - 1.4589869e+07i

EDIT3: Found this in the sqrt-doc of matlab:

%SQRT   Square root of fi object, computed using a bisection algorithm
%   C = SQRT(A) returns the square root of fi object A. Intermediate
%   quantities are calculated using the fimath associated with A.
%   The numerictype object of C is determined automatically for you using
%   an internal rule (see below).

[...]

%   Internal Rule:
%   For syntaxes where the numerictype object of the output is not 
%   specified as an input to the sqrt function, it is automatically 
%   calculated according to the following internal rule:
%     sign(C) = sign(A)
%     word-length(C) = ceil((word-length(A)/2))
%     fraction-length(C) = word-length(C) - 
%                           ceil(((word-length(A) - fraction-length(A))/2))

Yes, Matlab uses the principal square root. Why do you say it doesn't work for the second case? — Robert Israel, Sep 11 '15 at 16:58
@RobertIsrael: Please see my Edit. I am not sure if i can catch all cases by this construct... — Roland, Sep 11 '15 at 17:04
You want the real part $\ge 0$. The sign of the imaginary part only matters when the real part is $0$. — Robert Israel, Sep 11 '15 at 17:06
Analysing the real part does not seem to work for me... This does NOT work: [Code] boolean sHavePositiveRealPart = (s.real >= 0); s = sqrt(s); if (!sHavePositiveRealPart && s.real >= 0) s = -1 * s; [/Code] — Roland, Sep 11 '15 at 17:16
If you're doing this in C, don't write it yourself. Use complex.h and csqrt. — horchler, Sep 11 '15 at 17:45
It is supposed to be used on an embedded Hardware. I am not allowed to use complex.h, so i have to write the function on my own. — Roland, Sep 11 '15 at 17:47
I don't see any reason that complex.h wouldn't work on an embedded platform. In any case, you can implement a form of csqrt yourself. Note that this won't be exactly like Matlab's implementation, which is doing several things to be numerically careful and handling special cases. — horchler, Sep 11 '15 at 17:56
The Implementation in your link is one of the alternatives i have already tried. Is there any way to see the Code of the matlab sqrt function? — Roland, Sep 11 '15 at 17:59
You want to change the sign if $s.real < 0$ or ($s.real = 0$ and $s.imag < 0$). — Robert Israel, Sep 11 '15 at 18:17
@RobertIsrael: Sorry, that does not work either. [Code]boolean changeSign = (s.real < 0 || (s.imag < 0 && (abs(s.real) <= precision)))[/Code] does not filter all cases... :-( — Roland, Sep 11 '15 at 18:26
@RobertIsrael: The 2nd Case does not work. Please have a look at my EDIT3, too. — Roland, Sep 11 '15 at 19:10
You mean s = 1.2809769e+28 - 7.7697397e+27i ? The two square roots are 1.1788e+14 - 3.2956e+13i and -1.1788e+14 + 3.2956e+13i. The one you want is the one with positive real part, namely 1.1788e+14 - 3.2956e+13i. If you got the other one, multiply by $-1$. — Robert Israel, Sep 11 '15 at 19:17
No, i meant the case s=-1.7648506 - 0.0478944i. SquareRootOfSInMatlab = 0.0180244 - 1.3285991i; SquareRootOfSInC = -0.0182260163+1.33396506i! Your Condition (real < 0 || (real == 0 && imag < 0)) caused multiplication by -1 which was obviously wrong — Roland, Sep 11 '15 at 19:28
Your Edit 3: "Square root of fi object" – that documentation is for the fixed-point toolbox. You appear to be using floating point. Is there any way to see the Code of the matlab sqrt function? Sure, use coder/codegen to compile a function that take the square root of an arbitrary complex value. Then look at and decipher the generated code. — horchler, Sep 11 '15 at 22:14
How did you get -0.018226013+1.33396506i? If it was originally that, the real part is negative so you would multiply by $-1$ and get the correct result. If you started with the correct result, the real part is positive so you would not multiply by $-1$. — Robert Israel, Sep 13 '15 at 03:18
Oh! I think you're taking the sign of the real part before the square root. First take the square root, then look at the real part of that square root. — Robert Israel, Sep 13 '15 at 03:22

score 1 · Answer 1 · edited Jun 12 '20 at 10:38

1

I run octave right now which should be Matlab compatible and vast majority of basic algorithms are well tested to be Matlab compatible.

Using the code:

[xs,ys] = meshgrid(-16:15,-16:15);

figure(1); imagesc(angle((xs-i*ys))); colormap gray; colorbar;

title('Argument of complex numbers.');

figure(2); imagesc(angle(sqrt(xs-i*ys))); colormap gray; colorbar;

title('Argument of square root.');

The plots below indicate that the branch cut is along the negative real line.

edited Jun 12 '20 at 10:38

Community

1

answered Sep 11 '15 at 16:43

mathreadler

25,824

So, i can do it like that? [code] boolean sHavePositiveImaginaryPart = (s.imag > 0); s = sqrt(s); if (sHavePositiveImaginaryPart && s.imag > 0) s = -1*s; [/code]
Sorry, cannot Format the snippet properly :-(
– Roland Sep 11 '15 at 16:58
You seem to be trying to avoid having a positive imaginary part. Why would you do that? – Robert Israel Sep 11 '15 at 17:03

score 0 · Answer 2 · answered Sep 11 '15 at 16:37

0

Maybe it uses the polar form: $z^{1/2}=(re^{i\operatorname{Arg}z})^{1/2}=r^{1/2}e^{i\operatorname{Arg}z/2}$.

To do this, you would need to compute a real square root, the argument of $z$, and $e^x$ (modulo some more FLOPs for multiplies).

answered Sep 11 '15 at 16:37

parsiad

25,154

How calculates Matlab the sqrt of a complex number?

2 Answers2