1
2
3
4
5
6
7
8
9
10 csqrt_: proc (number) returns (complex float bin (27));
11
12 dcl (number, a, b) complex float bin (27);
13
14 a = number;
15 real (b) = abs (real (a));
16
17 if imag (a) = 0.0e0
18 then do;
19 real (b) = sqrt (real (b));
20 imag (b) = 0.0e0;
21 end;
22 else do;
23 real (b) = sqrt ((abs (a) + real (b))*0.5e0);
24
25 if real (a)<0.0e0
26 & imag (a)<0.0e0
27 then real (b) = -real (b);
28
29 imag (b) = imag (a) * 0.5e0/real (b);
30 end;
31
32 if real (a) >= 0.0e0
33 then return (b);
34
35 real (a) = imag (b);
36 imag (a) = real (b);
37
38 return (a);
39
40 end csqrt_;