|
111
|
1 ------------------------------------------------------------------------------
|
|
|
2 -- --
|
|
|
3 -- GNAT RUN-TIME COMPONENTS --
|
|
|
4 -- --
|
|
|
5 -- A D A . N U M E R I C S . G E N E R I C _ C O M P L E X _ T Y P E S --
|
|
|
6 -- --
|
|
|
7 -- B o d y --
|
|
|
8 -- --
|
|
131
|
9 -- Copyright (C) 1992-2018, Free Software Foundation, Inc. --
|
|
111
|
10 -- --
|
|
|
11 -- GNAT is free software; you can redistribute it and/or modify it under --
|
|
|
12 -- terms of the GNU General Public License as published by the Free Soft- --
|
|
|
13 -- ware Foundation; either version 3, or (at your option) any later ver- --
|
|
|
14 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
|
|
|
15 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
|
|
|
16 -- or FITNESS FOR A PARTICULAR PURPOSE. --
|
|
|
17 -- --
|
|
|
18 -- As a special exception under Section 7 of GPL version 3, you are granted --
|
|
|
19 -- additional permissions described in the GCC Runtime Library Exception, --
|
|
|
20 -- version 3.1, as published by the Free Software Foundation. --
|
|
|
21 -- --
|
|
|
22 -- You should have received a copy of the GNU General Public License and --
|
|
|
23 -- a copy of the GCC Runtime Library Exception along with this program; --
|
|
|
24 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
|
|
|
25 -- <http://www.gnu.org/licenses/>. --
|
|
|
26 -- --
|
|
|
27 -- GNAT was originally developed by the GNAT team at New York University. --
|
|
|
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
|
|
|
29 -- --
|
|
|
30 ------------------------------------------------------------------------------
|
|
|
31
|
|
|
32 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
|
|
|
33
|
|
|
34 package body Ada.Numerics.Generic_Complex_Types is
|
|
|
35
|
|
|
36 subtype R is Real'Base;
|
|
|
37
|
|
|
38 Two_Pi : constant R := R (2.0) * Pi;
|
|
|
39 Half_Pi : constant R := Pi / R (2.0);
|
|
|
40
|
|
|
41 ---------
|
|
|
42 -- "*" --
|
|
|
43 ---------
|
|
|
44
|
|
|
45 function "*" (Left, Right : Complex) return Complex is
|
|
|
46
|
|
|
47 Scale : constant R := R (R'Machine_Radix) ** ((R'Machine_Emax - 1) / 2);
|
|
|
48 -- In case of overflow, scale the operands by the largest power of the
|
|
|
49 -- radix (to avoid rounding error), so that the square of the scale does
|
|
|
50 -- not overflow itself.
|
|
|
51
|
|
|
52 X : R;
|
|
|
53 Y : R;
|
|
|
54
|
|
|
55 begin
|
|
|
56 X := Left.Re * Right.Re - Left.Im * Right.Im;
|
|
|
57 Y := Left.Re * Right.Im + Left.Im * Right.Re;
|
|
|
58
|
|
|
59 -- If either component overflows, try to scale (skip in fast math mode)
|
|
|
60
|
|
|
61 if not Standard'Fast_Math then
|
|
|
62
|
|
|
63 -- Note that the test below is written as a negation. This is to
|
|
|
64 -- account for the fact that X and Y may be NaNs, because both of
|
|
|
65 -- their operands could overflow. Given that all operations on NaNs
|
|
|
66 -- return false, the test can only be written thus.
|
|
|
67
|
|
|
68 if not (abs (X) <= R'Last) then
|
|
|
69 X := Scale**2 * ((Left.Re / Scale) * (Right.Re / Scale) -
|
|
|
70 (Left.Im / Scale) * (Right.Im / Scale));
|
|
|
71 end if;
|
|
|
72
|
|
|
73 if not (abs (Y) <= R'Last) then
|
|
|
74 Y := Scale**2 * ((Left.Re / Scale) * (Right.Im / Scale)
|
|
|
75 + (Left.Im / Scale) * (Right.Re / Scale));
|
|
|
76 end if;
|
|
|
77 end if;
|
|
|
78
|
|
|
79 return (X, Y);
|
|
|
80 end "*";
|
|
|
81
|
|
|
82 function "*" (Left, Right : Imaginary) return Real'Base is
|
|
|
83 begin
|
|
|
84 return -(R (Left) * R (Right));
|
|
|
85 end "*";
|
|
|
86
|
|
|
87 function "*" (Left : Complex; Right : Real'Base) return Complex is
|
|
|
88 begin
|
|
|
89 return Complex'(Left.Re * Right, Left.Im * Right);
|
|
|
90 end "*";
|
|
|
91
|
|
|
92 function "*" (Left : Real'Base; Right : Complex) return Complex is
|
|
|
93 begin
|
|
|
94 return (Left * Right.Re, Left * Right.Im);
|
|
|
95 end "*";
|
|
|
96
|
|
|
97 function "*" (Left : Complex; Right : Imaginary) return Complex is
|
|
|
98 begin
|
|
|
99 return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
|
|
|
100 end "*";
|
|
|
101
|
|
|
102 function "*" (Left : Imaginary; Right : Complex) return Complex is
|
|
|
103 begin
|
|
|
104 return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
|
|
|
105 end "*";
|
|
|
106
|
|
|
107 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
|
|
|
108 begin
|
|
|
109 return Left * Imaginary (Right);
|
|
|
110 end "*";
|
|
|
111
|
|
|
112 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
|
|
|
113 begin
|
|
|
114 return Imaginary (Left * R (Right));
|
|
|
115 end "*";
|
|
|
116
|
|
|
117 ----------
|
|
|
118 -- "**" --
|
|
|
119 ----------
|
|
|
120
|
|
|
121 function "**" (Left : Complex; Right : Integer) return Complex is
|
|
|
122 Result : Complex := (1.0, 0.0);
|
|
|
123 Factor : Complex := Left;
|
|
|
124 Exp : Integer := Right;
|
|
|
125
|
|
|
126 begin
|
|
|
127 -- We use the standard logarithmic approach, Exp gets shifted right
|
|
|
128 -- testing successive low order bits and Factor is the value of the
|
|
|
129 -- base raised to the next power of 2. For positive exponents we
|
|
|
130 -- multiply the result by this factor, for negative exponents, we
|
|
|
131 -- divide by this factor.
|
|
|
132
|
|
|
133 if Exp >= 0 then
|
|
|
134
|
|
|
135 -- For a positive exponent, if we get a constraint error during
|
|
|
136 -- this loop, it is an overflow, and the constraint error will
|
|
|
137 -- simply be passed on to the caller.
|
|
|
138
|
|
|
139 while Exp /= 0 loop
|
|
|
140 if Exp rem 2 /= 0 then
|
|
|
141 Result := Result * Factor;
|
|
|
142 end if;
|
|
|
143
|
|
|
144 Factor := Factor * Factor;
|
|
|
145 Exp := Exp / 2;
|
|
|
146 end loop;
|
|
|
147
|
|
|
148 return Result;
|
|
|
149
|
|
|
150 else -- Exp < 0 then
|
|
|
151
|
|
|
152 -- For the negative exponent case, a constraint error during this
|
|
|
153 -- calculation happens if Factor gets too large, and the proper
|
|
|
154 -- response is to return 0.0, since what we essentially have is
|
|
|
155 -- 1.0 / infinity, and the closest model number will be zero.
|
|
|
156
|
|
|
157 begin
|
|
|
158 while Exp /= 0 loop
|
|
|
159 if Exp rem 2 /= 0 then
|
|
|
160 Result := Result * Factor;
|
|
|
161 end if;
|
|
|
162
|
|
|
163 Factor := Factor * Factor;
|
|
|
164 Exp := Exp / 2;
|
|
|
165 end loop;
|
|
|
166
|
|
|
167 return R'(1.0) / Result;
|
|
|
168
|
|
|
169 exception
|
|
|
170 when Constraint_Error =>
|
|
|
171 return (0.0, 0.0);
|
|
|
172 end;
|
|
|
173 end if;
|
|
|
174 end "**";
|
|
|
175
|
|
|
176 function "**" (Left : Imaginary; Right : Integer) return Complex is
|
|
|
177 M : constant R := R (Left) ** Right;
|
|
|
178 begin
|
|
|
179 case Right mod 4 is
|
|
|
180 when 0 => return (M, 0.0);
|
|
|
181 when 1 => return (0.0, M);
|
|
|
182 when 2 => return (-M, 0.0);
|
|
|
183 when 3 => return (0.0, -M);
|
|
|
184 when others => raise Program_Error;
|
|
|
185 end case;
|
|
|
186 end "**";
|
|
|
187
|
|
|
188 ---------
|
|
|
189 -- "+" --
|
|
|
190 ---------
|
|
|
191
|
|
|
192 function "+" (Right : Complex) return Complex is
|
|
|
193 begin
|
|
|
194 return Right;
|
|
|
195 end "+";
|
|
|
196
|
|
|
197 function "+" (Left, Right : Complex) return Complex is
|
|
|
198 begin
|
|
|
199 return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
|
|
|
200 end "+";
|
|
|
201
|
|
|
202 function "+" (Right : Imaginary) return Imaginary is
|
|
|
203 begin
|
|
|
204 return Right;
|
|
|
205 end "+";
|
|
|
206
|
|
|
207 function "+" (Left, Right : Imaginary) return Imaginary is
|
|
|
208 begin
|
|
|
209 return Imaginary (R (Left) + R (Right));
|
|
|
210 end "+";
|
|
|
211
|
|
|
212 function "+" (Left : Complex; Right : Real'Base) return Complex is
|
|
|
213 begin
|
|
|
214 return Complex'(Left.Re + Right, Left.Im);
|
|
|
215 end "+";
|
|
|
216
|
|
|
217 function "+" (Left : Real'Base; Right : Complex) return Complex is
|
|
|
218 begin
|
|
|
219 return Complex'(Left + Right.Re, Right.Im);
|
|
|
220 end "+";
|
|
|
221
|
|
|
222 function "+" (Left : Complex; Right : Imaginary) return Complex is
|
|
|
223 begin
|
|
|
224 return Complex'(Left.Re, Left.Im + R (Right));
|
|
|
225 end "+";
|
|
|
226
|
|
|
227 function "+" (Left : Imaginary; Right : Complex) return Complex is
|
|
|
228 begin
|
|
|
229 return Complex'(Right.Re, R (Left) + Right.Im);
|
|
|
230 end "+";
|
|
|
231
|
|
|
232 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
|
|
|
233 begin
|
|
|
234 return Complex'(Right, R (Left));
|
|
|
235 end "+";
|
|
|
236
|
|
|
237 function "+" (Left : Real'Base; Right : Imaginary) return Complex is
|
|
|
238 begin
|
|
|
239 return Complex'(Left, R (Right));
|
|
|
240 end "+";
|
|
|
241
|
|
|
242 ---------
|
|
|
243 -- "-" --
|
|
|
244 ---------
|
|
|
245
|
|
|
246 function "-" (Right : Complex) return Complex is
|
|
|
247 begin
|
|
|
248 return (-Right.Re, -Right.Im);
|
|
|
249 end "-";
|
|
|
250
|
|
|
251 function "-" (Left, Right : Complex) return Complex is
|
|
|
252 begin
|
|
|
253 return (Left.Re - Right.Re, Left.Im - Right.Im);
|
|
|
254 end "-";
|
|
|
255
|
|
|
256 function "-" (Right : Imaginary) return Imaginary is
|
|
|
257 begin
|
|
|
258 return Imaginary (-R (Right));
|
|
|
259 end "-";
|
|
|
260
|
|
|
261 function "-" (Left, Right : Imaginary) return Imaginary is
|
|
|
262 begin
|
|
|
263 return Imaginary (R (Left) - R (Right));
|
|
|
264 end "-";
|
|
|
265
|
|
|
266 function "-" (Left : Complex; Right : Real'Base) return Complex is
|
|
|
267 begin
|
|
|
268 return Complex'(Left.Re - Right, Left.Im);
|
|
|
269 end "-";
|
|
|
270
|
|
|
271 function "-" (Left : Real'Base; Right : Complex) return Complex is
|
|
|
272 begin
|
|
|
273 return Complex'(Left - Right.Re, -Right.Im);
|
|
|
274 end "-";
|
|
|
275
|
|
|
276 function "-" (Left : Complex; Right : Imaginary) return Complex is
|
|
|
277 begin
|
|
|
278 return Complex'(Left.Re, Left.Im - R (Right));
|
|
|
279 end "-";
|
|
|
280
|
|
|
281 function "-" (Left : Imaginary; Right : Complex) return Complex is
|
|
|
282 begin
|
|
|
283 return Complex'(-Right.Re, R (Left) - Right.Im);
|
|
|
284 end "-";
|
|
|
285
|
|
|
286 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
|
|
|
287 begin
|
|
|
288 return Complex'(-Right, R (Left));
|
|
|
289 end "-";
|
|
|
290
|
|
|
291 function "-" (Left : Real'Base; Right : Imaginary) return Complex is
|
|
|
292 begin
|
|
|
293 return Complex'(Left, -R (Right));
|
|
|
294 end "-";
|
|
|
295
|
|
|
296 ---------
|
|
|
297 -- "/" --
|
|
|
298 ---------
|
|
|
299
|
|
|
300 function "/" (Left, Right : Complex) return Complex is
|
|
|
301 a : constant R := Left.Re;
|
|
|
302 b : constant R := Left.Im;
|
|
|
303 c : constant R := Right.Re;
|
|
|
304 d : constant R := Right.Im;
|
|
|
305
|
|
|
306 begin
|
|
|
307 if c = 0.0 and then d = 0.0 then
|
|
|
308 raise Constraint_Error;
|
|
|
309 else
|
|
|
310 return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
|
|
|
311 Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
|
|
|
312 end if;
|
|
|
313 end "/";
|
|
|
314
|
|
|
315 function "/" (Left, Right : Imaginary) return Real'Base is
|
|
|
316 begin
|
|
|
317 return R (Left) / R (Right);
|
|
|
318 end "/";
|
|
|
319
|
|
|
320 function "/" (Left : Complex; Right : Real'Base) return Complex is
|
|
|
321 begin
|
|
|
322 return Complex'(Left.Re / Right, Left.Im / Right);
|
|
|
323 end "/";
|
|
|
324
|
|
|
325 function "/" (Left : Real'Base; Right : Complex) return Complex is
|
|
|
326 a : constant R := Left;
|
|
|
327 c : constant R := Right.Re;
|
|
|
328 d : constant R := Right.Im;
|
|
|
329 begin
|
|
|
330 return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
|
|
|
331 Im => -((a * d) / (c ** 2 + d ** 2)));
|
|
|
332 end "/";
|
|
|
333
|
|
|
334 function "/" (Left : Complex; Right : Imaginary) return Complex is
|
|
|
335 a : constant R := Left.Re;
|
|
|
336 b : constant R := Left.Im;
|
|
|
337 d : constant R := R (Right);
|
|
|
338
|
|
|
339 begin
|
|
|
340 return (b / d, -(a / d));
|
|
|
341 end "/";
|
|
|
342
|
|
|
343 function "/" (Left : Imaginary; Right : Complex) return Complex is
|
|
|
344 b : constant R := R (Left);
|
|
|
345 c : constant R := Right.Re;
|
|
|
346 d : constant R := Right.Im;
|
|
|
347
|
|
|
348 begin
|
|
|
349 return (Re => b * d / (c ** 2 + d ** 2),
|
|
|
350 Im => b * c / (c ** 2 + d ** 2));
|
|
|
351 end "/";
|
|
|
352
|
|
|
353 function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
|
|
|
354 begin
|
|
|
355 return Imaginary (R (Left) / Right);
|
|
|
356 end "/";
|
|
|
357
|
|
|
358 function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
|
|
|
359 begin
|
|
|
360 return Imaginary (-(Left / R (Right)));
|
|
|
361 end "/";
|
|
|
362
|
|
|
363 ---------
|
|
|
364 -- "<" --
|
|
|
365 ---------
|
|
|
366
|
|
|
367 function "<" (Left, Right : Imaginary) return Boolean is
|
|
|
368 begin
|
|
|
369 return R (Left) < R (Right);
|
|
|
370 end "<";
|
|
|
371
|
|
|
372 ----------
|
|
|
373 -- "<=" --
|
|
|
374 ----------
|
|
|
375
|
|
|
376 function "<=" (Left, Right : Imaginary) return Boolean is
|
|
|
377 begin
|
|
|
378 return R (Left) <= R (Right);
|
|
|
379 end "<=";
|
|
|
380
|
|
|
381 ---------
|
|
|
382 -- ">" --
|
|
|
383 ---------
|
|
|
384
|
|
|
385 function ">" (Left, Right : Imaginary) return Boolean is
|
|
|
386 begin
|
|
|
387 return R (Left) > R (Right);
|
|
|
388 end ">";
|
|
|
389
|
|
|
390 ----------
|
|
|
391 -- ">=" --
|
|
|
392 ----------
|
|
|
393
|
|
|
394 function ">=" (Left, Right : Imaginary) return Boolean is
|
|
|
395 begin
|
|
|
396 return R (Left) >= R (Right);
|
|
|
397 end ">=";
|
|
|
398
|
|
|
399 -----------
|
|
|
400 -- "abs" --
|
|
|
401 -----------
|
|
|
402
|
|
|
403 function "abs" (Right : Imaginary) return Real'Base is
|
|
|
404 begin
|
|
|
405 return abs R (Right);
|
|
|
406 end "abs";
|
|
|
407
|
|
|
408 --------------
|
|
|
409 -- Argument --
|
|
|
410 --------------
|
|
|
411
|
|
|
412 function Argument (X : Complex) return Real'Base is
|
|
|
413 a : constant R := X.Re;
|
|
|
414 b : constant R := X.Im;
|
|
|
415 arg : R;
|
|
|
416
|
|
|
417 begin
|
|
|
418 if b = 0.0 then
|
|
|
419
|
|
|
420 if a >= 0.0 then
|
|
|
421 return 0.0;
|
|
|
422 else
|
|
|
423 return R'Copy_Sign (Pi, b);
|
|
|
424 end if;
|
|
|
425
|
|
|
426 elsif a = 0.0 then
|
|
|
427
|
|
|
428 if b >= 0.0 then
|
|
|
429 return Half_Pi;
|
|
|
430 else
|
|
|
431 return -Half_Pi;
|
|
|
432 end if;
|
|
|
433
|
|
|
434 else
|
|
|
435 arg := R (Atan (Double (abs (b / a))));
|
|
|
436
|
|
|
437 if a > 0.0 then
|
|
|
438 if b > 0.0 then
|
|
|
439 return arg;
|
|
|
440 else -- b < 0.0
|
|
|
441 return -arg;
|
|
|
442 end if;
|
|
|
443
|
|
|
444 else -- a < 0.0
|
|
|
445 if b >= 0.0 then
|
|
|
446 return Pi - arg;
|
|
|
447 else -- b < 0.0
|
|
|
448 return -(Pi - arg);
|
|
|
449 end if;
|
|
|
450 end if;
|
|
|
451 end if;
|
|
|
452
|
|
|
453 exception
|
|
|
454 when Constraint_Error =>
|
|
|
455 if b > 0.0 then
|
|
|
456 return Half_Pi;
|
|
|
457 else
|
|
|
458 return -Half_Pi;
|
|
|
459 end if;
|
|
|
460 end Argument;
|
|
|
461
|
|
|
462 function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
|
|
|
463 begin
|
|
|
464 if Cycle > 0.0 then
|
|
|
465 return Argument (X) * Cycle / Two_Pi;
|
|
|
466 else
|
|
|
467 raise Argument_Error;
|
|
|
468 end if;
|
|
|
469 end Argument;
|
|
|
470
|
|
|
471 ----------------------------
|
|
|
472 -- Compose_From_Cartesian --
|
|
|
473 ----------------------------
|
|
|
474
|
|
|
475 function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
|
|
|
476 begin
|
|
|
477 return (Re, Im);
|
|
|
478 end Compose_From_Cartesian;
|
|
|
479
|
|
|
480 function Compose_From_Cartesian (Re : Real'Base) return Complex is
|
|
|
481 begin
|
|
|
482 return (Re, 0.0);
|
|
|
483 end Compose_From_Cartesian;
|
|
|
484
|
|
|
485 function Compose_From_Cartesian (Im : Imaginary) return Complex is
|
|
|
486 begin
|
|
|
487 return (0.0, R (Im));
|
|
|
488 end Compose_From_Cartesian;
|
|
|
489
|
|
|
490 ------------------------
|
|
|
491 -- Compose_From_Polar --
|
|
|
492 ------------------------
|
|
|
493
|
|
|
494 function Compose_From_Polar (
|
|
|
495 Modulus, Argument : Real'Base)
|
|
|
496 return Complex
|
|
|
497 is
|
|
|
498 begin
|
|
|
499 if Modulus = 0.0 then
|
|
|
500 return (0.0, 0.0);
|
|
|
501 else
|
|
|
502 return (Modulus * R (Cos (Double (Argument))),
|
|
|
503 Modulus * R (Sin (Double (Argument))));
|
|
|
504 end if;
|
|
|
505 end Compose_From_Polar;
|
|
|
506
|
|
|
507 function Compose_From_Polar (
|
|
|
508 Modulus, Argument, Cycle : Real'Base)
|
|
|
509 return Complex
|
|
|
510 is
|
|
|
511 Arg : Real'Base;
|
|
|
512
|
|
|
513 begin
|
|
|
514 if Modulus = 0.0 then
|
|
|
515 return (0.0, 0.0);
|
|
|
516
|
|
|
517 elsif Cycle > 0.0 then
|
|
|
518 if Argument = 0.0 then
|
|
|
519 return (Modulus, 0.0);
|
|
|
520
|
|
|
521 elsif Argument = Cycle / 4.0 then
|
|
|
522 return (0.0, Modulus);
|
|
|
523
|
|
|
524 elsif Argument = Cycle / 2.0 then
|
|
|
525 return (-Modulus, 0.0);
|
|
|
526
|
|
|
527 elsif Argument = 3.0 * Cycle / R (4.0) then
|
|
|
528 return (0.0, -Modulus);
|
|
|
529 else
|
|
|
530 Arg := Two_Pi * Argument / Cycle;
|
|
|
531 return (Modulus * R (Cos (Double (Arg))),
|
|
|
532 Modulus * R (Sin (Double (Arg))));
|
|
|
533 end if;
|
|
|
534 else
|
|
|
535 raise Argument_Error;
|
|
|
536 end if;
|
|
|
537 end Compose_From_Polar;
|
|
|
538
|
|
|
539 ---------------
|
|
|
540 -- Conjugate --
|
|
|
541 ---------------
|
|
|
542
|
|
|
543 function Conjugate (X : Complex) return Complex is
|
|
|
544 begin
|
|
|
545 return Complex'(X.Re, -X.Im);
|
|
|
546 end Conjugate;
|
|
|
547
|
|
|
548 --------
|
|
|
549 -- Im --
|
|
|
550 --------
|
|
|
551
|
|
|
552 function Im (X : Complex) return Real'Base is
|
|
|
553 begin
|
|
|
554 return X.Im;
|
|
|
555 end Im;
|
|
|
556
|
|
|
557 function Im (X : Imaginary) return Real'Base is
|
|
|
558 begin
|
|
|
559 return R (X);
|
|
|
560 end Im;
|
|
|
561
|
|
|
562 -------------
|
|
|
563 -- Modulus --
|
|
|
564 -------------
|
|
|
565
|
|
|
566 function Modulus (X : Complex) return Real'Base is
|
|
|
567 Re2, Im2 : R;
|
|
|
568
|
|
|
569 begin
|
|
|
570
|
|
|
571 begin
|
|
|
572 Re2 := X.Re ** 2;
|
|
|
573
|
|
|
574 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
|
|
|
575 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
|
|
|
576 -- squaring does not raise constraint_error but generates infinity,
|
|
|
577 -- we can use an explicit comparison to determine whether to use
|
|
|
578 -- the scaling expression.
|
|
|
579
|
|
|
580 -- The scaling expression is computed in double format throughout
|
|
|
581 -- in order to prevent inaccuracies on machines where not all
|
|
|
582 -- immediate expressions are rounded, such as PowerPC.
|
|
|
583
|
|
|
584 -- ??? same weird test, why not Re2 > R'Last ???
|
|
|
585 if not (Re2 <= R'Last) then
|
|
|
586 raise Constraint_Error;
|
|
|
587 end if;
|
|
|
588
|
|
|
589 exception
|
|
|
590 when Constraint_Error =>
|
|
|
591 return R (Double (abs (X.Re))
|
|
|
592 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
|
|
|
593 end;
|
|
|
594
|
|
|
595 begin
|
|
|
596 Im2 := X.Im ** 2;
|
|
|
597
|
|
|
598 -- ??? same weird test
|
|
|
599 if not (Im2 <= R'Last) then
|
|
|
600 raise Constraint_Error;
|
|
|
601 end if;
|
|
|
602
|
|
|
603 exception
|
|
|
604 when Constraint_Error =>
|
|
|
605 return R (Double (abs (X.Im))
|
|
|
606 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
|
|
|
607 end;
|
|
|
608
|
|
|
609 -- Now deal with cases of underflow. If only one of the squares
|
|
|
610 -- underflows, return the modulus of the other component. If both
|
|
|
611 -- squares underflow, use scaling as above.
|
|
|
612
|
|
|
613 if Re2 = 0.0 then
|
|
|
614
|
|
|
615 if X.Re = 0.0 then
|
|
|
616 return abs (X.Im);
|
|
|
617
|
|
|
618 elsif Im2 = 0.0 then
|
|
|
619
|
|
|
620 if X.Im = 0.0 then
|
|
|
621 return abs (X.Re);
|
|
|
622
|
|
|
623 else
|
|
|
624 if abs (X.Re) > abs (X.Im) then
|
|
|
625 return
|
|
|
626 R (Double (abs (X.Re))
|
|
|
627 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
|
|
|
628 else
|
|
|
629 return
|
|
|
630 R (Double (abs (X.Im))
|
|
|
631 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
|
|
|
632 end if;
|
|
|
633 end if;
|
|
|
634
|
|
|
635 else
|
|
|
636 return abs (X.Im);
|
|
|
637 end if;
|
|
|
638
|
|
|
639 elsif Im2 = 0.0 then
|
|
|
640 return abs (X.Re);
|
|
|
641
|
|
|
642 -- In all other cases, the naive computation will do
|
|
|
643
|
|
|
644 else
|
|
|
645 return R (Sqrt (Double (Re2 + Im2)));
|
|
|
646 end if;
|
|
|
647 end Modulus;
|
|
|
648
|
|
|
649 --------
|
|
|
650 -- Re --
|
|
|
651 --------
|
|
|
652
|
|
|
653 function Re (X : Complex) return Real'Base is
|
|
|
654 begin
|
|
|
655 return X.Re;
|
|
|
656 end Re;
|
|
|
657
|
|
|
658 ------------
|
|
|
659 -- Set_Im --
|
|
|
660 ------------
|
|
|
661
|
|
|
662 procedure Set_Im (X : in out Complex; Im : Real'Base) is
|
|
|
663 begin
|
|
|
664 X.Im := Im;
|
|
|
665 end Set_Im;
|
|
|
666
|
|
|
667 procedure Set_Im (X : out Imaginary; Im : Real'Base) is
|
|
|
668 begin
|
|
|
669 X := Imaginary (Im);
|
|
|
670 end Set_Im;
|
|
|
671
|
|
|
672 ------------
|
|
|
673 -- Set_Re --
|
|
|
674 ------------
|
|
|
675
|
|
|
676 procedure Set_Re (X : in out Complex; Re : Real'Base) is
|
|
|
677 begin
|
|
|
678 X.Re := Re;
|
|
|
679 end Set_Re;
|
|
|
680
|
|
|
681 end Ada.Numerics.Generic_Complex_Types;
|
