1
2
3
4
5
6
7
8
9
10 dcsqrt_: proc (number) returns (complex float bin (63));
11
12 dcl (number, a, b) complex float bin (63);
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 then return (b);
33
34 real (a) = imag (b);
35 imag (a) = real (b);
36
37 return (a);
38
39 end dcsqrt_;