0
$\begingroup$

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))
$\endgroup$
  • $\begingroup$ Yes, Matlab uses the principal square root. Why do you say it doesn't work for the second case? $\endgroup$ – Robert Israel Sep 11 '15 at 16:58
  • $\begingroup$ @RobertIsrael: Please see my Edit. I am not sure if i can catch all cases by this construct... $\endgroup$ – Roland Sep 11 '15 at 17:04
  • $\begingroup$ You want the real part $\ge 0$. The sign of the imaginary part only matters when the real part is $0$. $\endgroup$ – Robert Israel Sep 11 '15 at 17:06
  • $\begingroup$ 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] $\endgroup$ – Roland Sep 11 '15 at 17:16
  • $\begingroup$ If you're doing this in C, don't write it yourself. Use complex.h and csqrt. $\endgroup$ – horchler Sep 11 '15 at 17:45
1
$\begingroup$

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.

enter image description here enter image description here

$\endgroup$
  • $\begingroup$ 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 :-( $\endgroup$ – Roland Sep 11 '15 at 16:58
  • $\begingroup$ You seem to be trying to avoid having a positive imaginary part. Why would you do that? $\endgroup$ – Robert Israel Sep 11 '15 at 17:03
0
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.