1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME 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 --
11 -- Copyright (C) 1992-2001 Free Software Foundation, Inc. --
13 -- GNAT is free software; you can redistribute it and/or modify it under --
14 -- terms of the GNU General Public License as published by the Free Soft- --
15 -- ware Foundation; either version 2, or (at your option) any later ver- --
16 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
17 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
18 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
19 -- for more details. You should have received a copy of the GNU General --
20 -- Public License distributed with GNAT; see file COPYING. If not, write --
21 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
22 -- MA 02111-1307, USA. --
24 -- As a special exception, if other files instantiate generics from this --
25 -- unit, or you link this unit with other files to produce an executable, --
26 -- this unit does not by itself cause the resulting executable to be --
27 -- covered by the GNU General Public License. This exception does not --
28 -- however invalidate any other reasons why the executable file might be --
29 -- covered by the GNU Public License. --
31 -- GNAT was originally developed by the GNAT team at New York University. --
32 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
34 ------------------------------------------------------------------------------
36 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
37 package body Ada.Numerics.Generic_Complex_Types is
39 subtype R is Real'Base;
41 Two_Pi : constant R := R (2.0) * Pi;
42 Half_Pi : constant R := Pi / R (2.0);
48 function "*" (Left, Right : Complex) return Complex is
53 X := Left.Re * Right.Re - Left.Im * Right.Im;
54 Y := Left.Re * Right.Im + Left.Im * Right.Re;
56 -- If either component overflows, try to scale.
58 if abs (X) > R'Last then
59 X := R' (4.0) * (R'(Left.Re / 2.0) * R'(Right.Re / 2.0)
60 - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
63 if abs (Y) > R'Last then
64 Y := R' (4.0) * (R'(Left.Re / 2.0) * R'(Right.Im / 2.0)
65 - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
71 function "*" (Left, Right : Imaginary) return Real'Base is
73 return -R (Left) * R (Right);
76 function "*" (Left : Complex; Right : Real'Base) return Complex is
78 return Complex'(Left.Re * Right, Left.Im * Right);
81 function "*" (Left : Real'Base; Right : Complex) return Complex is
83 return (Left * Right.Re, Left * Right.Im);
86 function "*" (Left : Complex; Right : Imaginary) return Complex is
88 return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
91 function "*" (Left : Imaginary; Right : Complex) return Complex is
93 return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
96 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
98 return Left * Imaginary (Right);
101 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
103 return Imaginary (Left * R (Right));
110 function "**" (Left : Complex; Right : Integer) return Complex is
111 Result : Complex := (1.0, 0.0);
112 Factor : Complex := Left;
113 Exp : Integer := Right;
116 -- We use the standard logarithmic approach, Exp gets shifted right
117 -- testing successive low order bits and Factor is the value of the
118 -- base raised to the next power of 2. For positive exponents we
119 -- multiply the result by this factor, for negative exponents, we
120 -- divide by this factor.
124 -- For a positive exponent, if we get a constraint error during
125 -- this loop, it is an overflow, and the constraint error will
126 -- simply be passed on to the caller.
129 if Exp rem 2 /= 0 then
130 Result := Result * Factor;
133 Factor := Factor * Factor;
141 -- For the negative exponent case, a constraint error during this
142 -- calculation happens if Factor gets too large, and the proper
143 -- response is to return 0.0, since what we essentially have is
144 -- 1.0 / infinity, and the closest model number will be zero.
149 if Exp rem 2 /= 0 then
150 Result := Result * Factor;
153 Factor := Factor * Factor;
157 return R ' (1.0) / Result;
161 when Constraint_Error =>
167 function "**" (Left : Imaginary; Right : Integer) return Complex is
168 M : R := R (Left) ** Right;
171 when 0 => return (M, 0.0);
172 when 1 => return (0.0, M);
173 when 2 => return (-M, 0.0);
174 when 3 => return (0.0, -M);
175 when others => raise Program_Error;
183 function "+" (Right : Complex) return Complex is
188 function "+" (Left, Right : Complex) return Complex is
190 return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
193 function "+" (Right : Imaginary) return Imaginary is
198 function "+" (Left, Right : Imaginary) return Imaginary is
200 return Imaginary (R (Left) + R (Right));
203 function "+" (Left : Complex; Right : Real'Base) return Complex is
205 return Complex'(Left.Re + Right, Left.Im);
208 function "+" (Left : Real'Base; Right : Complex) return Complex is
210 return Complex'(Left + Right.Re, Right.Im);
213 function "+" (Left : Complex; Right : Imaginary) return Complex is
215 return Complex'(Left.Re, Left.Im + R (Right));
218 function "+" (Left : Imaginary; Right : Complex) return Complex is
220 return Complex'(Right.Re, R (Left) + Right.Im);
223 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
225 return Complex'(Right, R (Left));
228 function "+" (Left : Real'Base; Right : Imaginary) return Complex is
230 return Complex'(Left, R (Right));
237 function "-" (Right : Complex) return Complex is
239 return (-Right.Re, -Right.Im);
242 function "-" (Left, Right : Complex) return Complex is
244 return (Left.Re - Right.Re, Left.Im - Right.Im);
247 function "-" (Right : Imaginary) return Imaginary is
249 return Imaginary (-R (Right));
252 function "-" (Left, Right : Imaginary) return Imaginary is
254 return Imaginary (R (Left) - R (Right));
257 function "-" (Left : Complex; Right : Real'Base) return Complex is
259 return Complex'(Left.Re - Right, Left.Im);
262 function "-" (Left : Real'Base; Right : Complex) return Complex is
264 return Complex'(Left - Right.Re, -Right.Im);
267 function "-" (Left : Complex; Right : Imaginary) return Complex is
269 return Complex'(Left.Re, Left.Im - R (Right));
272 function "-" (Left : Imaginary; Right : Complex) return Complex is
274 return Complex'(-Right.Re, R (Left) - Right.Im);
277 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
279 return Complex'(-Right, R (Left));
282 function "-" (Left : Real'Base; Right : Imaginary) return Complex is
284 return Complex'(Left, -R (Right));
291 function "/" (Left, Right : Complex) return Complex is
292 a : constant R := Left.Re;
293 b : constant R := Left.Im;
294 c : constant R := Right.Re;
295 d : constant R := Right.Im;
298 if c = 0.0 and then d = 0.0 then
299 raise Constraint_Error;
301 return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
302 Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
306 function "/" (Left, Right : Imaginary) return Real'Base is
308 return R (Left) / R (Right);
311 function "/" (Left : Complex; Right : Real'Base) return Complex is
313 return Complex'(Left.Re / Right, Left.Im / Right);
316 function "/" (Left : Real'Base; Right : Complex) return Complex is
317 a : constant R := Left;
318 c : constant R := Right.Re;
319 d : constant R := Right.Im;
321 return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
322 Im => -(a * d) / (c ** 2 + d ** 2));
325 function "/" (Left : Complex; Right : Imaginary) return Complex is
326 a : constant R := Left.Re;
327 b : constant R := Left.Im;
328 d : constant R := R (Right);
331 return (b / d, -a / d);
334 function "/" (Left : Imaginary; Right : Complex) return Complex is
335 b : constant R := R (Left);
336 c : constant R := Right.Re;
337 d : constant R := Right.Im;
340 return (Re => b * d / (c ** 2 + d ** 2),
341 Im => b * c / (c ** 2 + d ** 2));
344 function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
346 return Imaginary (R (Left) / Right);
349 function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
351 return Imaginary (-Left / R (Right));
358 function "<" (Left, Right : Imaginary) return Boolean is
360 return R (Left) < R (Right);
367 function "<=" (Left, Right : Imaginary) return Boolean is
369 return R (Left) <= R (Right);
376 function ">" (Left, Right : Imaginary) return Boolean is
378 return R (Left) > R (Right);
385 function ">=" (Left, Right : Imaginary) return Boolean is
387 return R (Left) >= R (Right);
394 function "abs" (Right : Imaginary) return Real'Base is
396 return abs R (Right);
403 function Argument (X : Complex) return Real'Base is
404 a : constant R := X.Re;
405 b : constant R := X.Im;
414 return R'Copy_Sign (Pi, b);
426 arg := R (Atan (Double (abs (b / a))));
445 when Constraint_Error =>
453 function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
456 return Argument (X) * Cycle / Two_Pi;
458 raise Argument_Error;
462 ----------------------------
463 -- Compose_From_Cartesian --
464 ----------------------------
466 function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
469 end Compose_From_Cartesian;
471 function Compose_From_Cartesian (Re : Real'Base) return Complex is
474 end Compose_From_Cartesian;
476 function Compose_From_Cartesian (Im : Imaginary) return Complex is
478 return (0.0, R (Im));
479 end Compose_From_Cartesian;
481 ------------------------
482 -- Compose_From_Polar --
483 ------------------------
485 function Compose_From_Polar (
486 Modulus, Argument : Real'Base)
490 if Modulus = 0.0 then
493 return (Modulus * R (Cos (Double (Argument))),
494 Modulus * R (Sin (Double (Argument))));
496 end Compose_From_Polar;
498 function Compose_From_Polar (
499 Modulus, Argument, Cycle : Real'Base)
505 if Modulus = 0.0 then
508 elsif Cycle > 0.0 then
509 if Argument = 0.0 then
510 return (Modulus, 0.0);
512 elsif Argument = Cycle / 4.0 then
513 return (0.0, Modulus);
515 elsif Argument = Cycle / 2.0 then
516 return (-Modulus, 0.0);
518 elsif Argument = 3.0 * Cycle / R (4.0) then
519 return (0.0, -Modulus);
521 Arg := Two_Pi * Argument / Cycle;
522 return (Modulus * R (Cos (Double (Arg))),
523 Modulus * R (Sin (Double (Arg))));
526 raise Argument_Error;
528 end Compose_From_Polar;
534 function Conjugate (X : Complex) return Complex is
536 return Complex'(X.Re, -X.Im);
543 function Im (X : Complex) return Real'Base is
548 function Im (X : Imaginary) return Real'Base is
557 function Modulus (X : Complex) return Real'Base is
565 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
566 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
567 -- squaring does not raise constraint_error but generates infinity,
568 -- we can use an explicit comparison to determine whether to use
569 -- the scaling expression.
572 raise Constraint_Error;
576 when Constraint_Error =>
578 * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
585 raise Constraint_Error;
589 when Constraint_Error =>
591 * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
594 -- Now deal with cases of underflow. If only one of the squares
595 -- underflows, return the modulus of the other component. If both
596 -- squares underflow, use scaling as above.
609 if abs (X.Re) > abs (X.Im) then
612 * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
616 * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
628 -- in all other cases, the naive computation will do.
631 return R (Sqrt (Double (Re2 + Im2)));
639 function Re (X : Complex) return Real'Base is
648 procedure Set_Im (X : in out Complex; Im : in Real'Base) is
653 procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
662 procedure Set_Re (X : in out Complex; Re : in Real'Base) is
667 end Ada.Numerics.Generic_Complex_Types;