1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
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 --
9 -- Copyright (C) 1992-2020, Free Software Foundation, Inc. --
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. --
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. --
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/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
32 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
34 package body Ada.Numerics.Generic_Complex_Types is
36 subtype R is Real'Base;
38 Two_Pi : constant R := R (2.0) * Pi;
39 Half_Pi : constant R := Pi / R (2.0);
45 function "*" (Left, Right : Complex) return Complex is
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.
56 X := Left.Re * Right.Re - Left.Im * Right.Im;
57 Y := Left.Re * Right.Im + Left.Im * Right.Re;
59 -- If either component overflows, try to scale (skip in fast math mode)
61 if not Standard'Fast_Math then
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.
68 if not (abs (X) <= R'Last) then
70 (CodePeer, Intentional,
71 "test always false", "test for infinity");
73 X := Scale**2 * ((Left.Re / Scale) * (Right.Re / Scale) -
74 (Left.Im / Scale) * (Right.Im / Scale));
77 if not (abs (Y) <= R'Last) then
79 (CodePeer, Intentional,
80 "test always false", "test for infinity");
82 Y := Scale**2 * ((Left.Re / Scale) * (Right.Im / Scale)
83 + (Left.Im / Scale) * (Right.Re / Scale));
90 function "*" (Left, Right : Imaginary) return Real'Base is
92 return -(R (Left) * R (Right));
95 function "*" (Left : Complex; Right : Real'Base) return Complex is
97 return Complex'(Left.Re * Right, Left.Im * Right);
100 function "*" (Left : Real'Base; Right : Complex) return Complex is
102 return (Left * Right.Re, Left * Right.Im);
105 function "*" (Left : Complex; Right : Imaginary) return Complex is
107 return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
110 function "*" (Left : Imaginary; Right : Complex) return Complex is
112 return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
115 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
117 return Left * Imaginary (Right);
120 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
122 return Imaginary (Left * R (Right));
129 function "**" (Left : Complex; Right : Integer) return Complex is
130 Result : Complex := (1.0, 0.0);
131 Factor : Complex := Left;
132 Exp : Integer := Right;
135 -- We use the standard logarithmic approach, Exp gets shifted right
136 -- testing successive low order bits and Factor is the value of the
137 -- base raised to the next power of 2. For positive exponents we
138 -- multiply the result by this factor, for negative exponents, we
139 -- divide by this factor.
143 -- For a positive exponent, if we get a constraint error during
144 -- this loop, it is an overflow, and the constraint error will
145 -- simply be passed on to the caller.
148 if Exp rem 2 /= 0 then
149 Result := Result * Factor;
152 Factor := Factor * Factor;
160 -- For the negative exponent case, a constraint error during this
161 -- calculation happens if Factor gets too large, and the proper
162 -- response is to return 0.0, since what we essentially have is
163 -- 1.0 / infinity, and the closest model number will be zero.
167 if Exp rem 2 /= 0 then
168 Result := Result * Factor;
171 Factor := Factor * Factor;
175 return R'(1.0) / Result;
178 when Constraint_Error =>
184 function "**" (Left : Imaginary; Right : Integer) return Complex is
185 M : constant R := R (Left) ** Right;
188 when 0 => return (M, 0.0);
189 when 1 => return (0.0, M);
190 when 2 => return (-M, 0.0);
191 when 3 => return (0.0, -M);
192 when others => raise Program_Error;
200 function "+" (Right : Complex) return Complex is
205 function "+" (Left, Right : Complex) return Complex is
207 return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
210 function "+" (Right : Imaginary) return Imaginary is
215 function "+" (Left, Right : Imaginary) return Imaginary is
217 return Imaginary (R (Left) + R (Right));
220 function "+" (Left : Complex; Right : Real'Base) return Complex is
222 return Complex'(Left.Re + Right, Left.Im);
225 function "+" (Left : Real'Base; Right : Complex) return Complex is
227 return Complex'(Left + Right.Re, Right.Im);
230 function "+" (Left : Complex; Right : Imaginary) return Complex is
232 return Complex'(Left.Re, Left.Im + R (Right));
235 function "+" (Left : Imaginary; Right : Complex) return Complex is
237 return Complex'(Right.Re, R (Left) + Right.Im);
240 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
242 return Complex'(Right, R (Left));
245 function "+" (Left : Real'Base; Right : Imaginary) return Complex is
247 return Complex'(Left, R (Right));
254 function "-" (Right : Complex) return Complex is
256 return (-Right.Re, -Right.Im);
259 function "-" (Left, Right : Complex) return Complex is
261 return (Left.Re - Right.Re, Left.Im - Right.Im);
264 function "-" (Right : Imaginary) return Imaginary is
266 return Imaginary (-R (Right));
269 function "-" (Left, Right : Imaginary) return Imaginary is
271 return Imaginary (R (Left) - R (Right));
274 function "-" (Left : Complex; Right : Real'Base) return Complex is
276 return Complex'(Left.Re - Right, Left.Im);
279 function "-" (Left : Real'Base; Right : Complex) return Complex is
281 return Complex'(Left - Right.Re, -Right.Im);
284 function "-" (Left : Complex; Right : Imaginary) return Complex is
286 return Complex'(Left.Re, Left.Im - R (Right));
289 function "-" (Left : Imaginary; Right : Complex) return Complex is
291 return Complex'(-Right.Re, R (Left) - Right.Im);
294 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
296 return Complex'(-Right, R (Left));
299 function "-" (Left : Real'Base; Right : Imaginary) return Complex is
301 return Complex'(Left, -R (Right));
308 function "/" (Left, Right : Complex) return Complex is
309 a : constant R := Left.Re;
310 b : constant R := Left.Im;
311 c : constant R := Right.Re;
312 d : constant R := Right.Im;
315 if c = 0.0 and then d = 0.0 then
316 raise Constraint_Error;
318 return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
319 Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
323 function "/" (Left, Right : Imaginary) return Real'Base is
325 return R (Left) / R (Right);
328 function "/" (Left : Complex; Right : Real'Base) return Complex is
330 return Complex'(Left.Re / Right, Left.Im / Right);
333 function "/" (Left : Real'Base; Right : Complex) return Complex is
334 a : constant R := Left;
335 c : constant R := Right.Re;
336 d : constant R := Right.Im;
338 return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
339 Im => -((a * d) / (c ** 2 + d ** 2)));
342 function "/" (Left : Complex; Right : Imaginary) return Complex is
343 a : constant R := Left.Re;
344 b : constant R := Left.Im;
345 d : constant R := R (Right);
348 return (b / d, -(a / d));
351 function "/" (Left : Imaginary; Right : Complex) return Complex is
352 b : constant R := R (Left);
353 c : constant R := Right.Re;
354 d : constant R := Right.Im;
357 return (Re => b * d / (c ** 2 + d ** 2),
358 Im => b * c / (c ** 2 + d ** 2));
361 function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
363 return Imaginary (R (Left) / Right);
366 function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
368 return Imaginary (-(Left / R (Right)));
375 function "<" (Left, Right : Imaginary) return Boolean is
377 return R (Left) < R (Right);
384 function "<=" (Left, Right : Imaginary) return Boolean is
386 return R (Left) <= R (Right);
393 function ">" (Left, Right : Imaginary) return Boolean is
395 return R (Left) > R (Right);
402 function ">=" (Left, Right : Imaginary) return Boolean is
404 return R (Left) >= R (Right);
411 function "abs" (Right : Imaginary) return Real'Base is
413 return abs R (Right);
420 function Argument (X : Complex) return Real'Base is
421 a : constant R := X.Re;
422 b : constant R := X.Im;
431 return R'Copy_Sign (Pi, b);
443 arg := R (Atan (Double (abs (b / a))));
462 when Constraint_Error =>
470 function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
473 return Argument (X) * Cycle / Two_Pi;
475 raise Argument_Error;
479 ----------------------------
480 -- Compose_From_Cartesian --
481 ----------------------------
483 function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
486 end Compose_From_Cartesian;
488 function Compose_From_Cartesian (Re : Real'Base) return Complex is
491 end Compose_From_Cartesian;
493 function Compose_From_Cartesian (Im : Imaginary) return Complex is
495 return (0.0, R (Im));
496 end Compose_From_Cartesian;
498 ------------------------
499 -- Compose_From_Polar --
500 ------------------------
502 function Compose_From_Polar (
503 Modulus, Argument : Real'Base)
507 if Modulus = 0.0 then
510 return (Modulus * R (Cos (Double (Argument))),
511 Modulus * R (Sin (Double (Argument))));
513 end Compose_From_Polar;
515 function Compose_From_Polar (
516 Modulus, Argument, Cycle : Real'Base)
522 if Modulus = 0.0 then
525 elsif Cycle > 0.0 then
526 if Argument = 0.0 then
527 return (Modulus, 0.0);
529 elsif Argument = Cycle / 4.0 then
530 return (0.0, Modulus);
532 elsif Argument = Cycle / 2.0 then
533 return (-Modulus, 0.0);
535 elsif Argument = 3.0 * Cycle / R (4.0) then
536 return (0.0, -Modulus);
538 Arg := Two_Pi * Argument / Cycle;
539 return (Modulus * R (Cos (Double (Arg))),
540 Modulus * R (Sin (Double (Arg))));
543 raise Argument_Error;
545 end Compose_From_Polar;
551 function Conjugate (X : Complex) return Complex is
553 return Complex'(X.Re, -X.Im);
560 function Im (X : Complex) return Real'Base is
565 function Im (X : Imaginary) return Real'Base is
574 function Modulus (X : Complex) return Real'Base is
582 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
583 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
584 -- squaring does not raise constraint_error but generates infinity,
585 -- we can use an explicit comparison to determine whether to use
586 -- the scaling expression.
588 -- The scaling expression is computed in double format throughout
589 -- in order to prevent inaccuracies on machines where not all
590 -- immediate expressions are rounded, such as PowerPC.
592 -- ??? same weird test, why not Re2 > R'Last ???
593 if not (Re2 <= R'Last) then
594 raise Constraint_Error;
598 when Constraint_Error =>
599 pragma Assert (X.Re /= 0.0);
600 return R (Double (abs (X.Re))
601 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
607 -- ??? same weird test
608 if not (Im2 <= R'Last) then
609 raise Constraint_Error;
613 when Constraint_Error =>
614 pragma Assert (X.Im /= 0.0);
615 return R (Double (abs (X.Im))
616 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
619 -- Now deal with cases of underflow. If only one of the squares
620 -- underflows, return the modulus of the other component. If both
621 -- squares underflow, use scaling as above.
634 if abs (X.Re) > abs (X.Im) then
636 R (Double (abs (X.Re))
637 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
640 R (Double (abs (X.Im))
641 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
652 -- In all other cases, the naive computation will do
655 return R (Sqrt (Double (Re2 + Im2)));
663 function Re (X : Complex) return Real'Base is
672 procedure Set_Im (X : in out Complex; Im : Real'Base) is
677 procedure Set_Im (X : out Imaginary; Im : Real'Base) is
686 procedure Set_Re (X : in out Complex; Re : Real'Base) is
691 end Ada.Numerics.Generic_Complex_Types;