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-2005 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 2, 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. See the GNU General Public License --
17 -- for more details. You should have received a copy of the GNU General --
18 -- Public License distributed with GNAT; see file COPYING. If not, write --
19 -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, USA. --
22 -- As a special exception, if other files instantiate generics from this --
23 -- unit, or you link this unit with other files to produce an executable, --
24 -- this unit does not by itself cause the resulting executable to be --
25 -- covered by the GNU General Public License. This exception does not --
26 -- however invalidate any other reasons why the executable file might be --
27 -- covered by the GNU Public License. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
34 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
35 package body Ada.Numerics.Generic_Complex_Types is
37 subtype R is Real'Base;
39 Two_Pi : constant R := R (2.0) * Pi;
40 Half_Pi : constant R := Pi / R (2.0);
46 function "*" (Left, Right : Complex) return Complex is
51 X := Left.Re * Right.Re - Left.Im * Right.Im;
52 Y := Left.Re * Right.Im + Left.Im * Right.Re;
54 -- If either component overflows, try to scale.
56 if abs (X) > R'Last then
57 X := R'(4.0) * (R'(Left.Re / 2.0) * R'(Right.Re / 2.0)
58 - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
61 if abs (Y) > R'Last then
62 Y := R'(4.0) * (R'(Left.Re / 2.0) * R'(Right.Im / 2.0)
63 - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
69 function "*" (Left, Right : Imaginary) return Real'Base is
71 return -R (Left) * R (Right);
74 function "*" (Left : Complex; Right : Real'Base) return Complex is
76 return Complex'(Left.Re * Right, Left.Im * Right);
79 function "*" (Left : Real'Base; Right : Complex) return Complex is
81 return (Left * Right.Re, Left * Right.Im);
84 function "*" (Left : Complex; Right : Imaginary) return Complex is
86 return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
89 function "*" (Left : Imaginary; Right : Complex) return Complex is
91 return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
94 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
96 return Left * Imaginary (Right);
99 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
101 return Imaginary (Left * R (Right));
108 function "**" (Left : Complex; Right : Integer) return Complex is
109 Result : Complex := (1.0, 0.0);
110 Factor : Complex := Left;
111 Exp : Integer := Right;
114 -- We use the standard logarithmic approach, Exp gets shifted right
115 -- testing successive low order bits and Factor is the value of the
116 -- base raised to the next power of 2. For positive exponents we
117 -- multiply the result by this factor, for negative exponents, we
118 -- divide by this factor.
122 -- For a positive exponent, if we get a constraint error during
123 -- this loop, it is an overflow, and the constraint error will
124 -- simply be passed on to the caller.
127 if Exp rem 2 /= 0 then
128 Result := Result * Factor;
131 Factor := Factor * Factor;
139 -- For the negative exponent case, a constraint error during this
140 -- calculation happens if Factor gets too large, and the proper
141 -- response is to return 0.0, since what we essentially have is
142 -- 1.0 / infinity, and the closest model number will be zero.
147 if Exp rem 2 /= 0 then
148 Result := Result * Factor;
151 Factor := Factor * Factor;
155 return R'(1.0) / Result;
159 when Constraint_Error =>
165 function "**" (Left : Imaginary; Right : Integer) return Complex is
166 M : constant R := R (Left) ** Right;
169 when 0 => return (M, 0.0);
170 when 1 => return (0.0, M);
171 when 2 => return (-M, 0.0);
172 when 3 => return (0.0, -M);
173 when others => raise Program_Error;
181 function "+" (Right : Complex) return Complex is
186 function "+" (Left, Right : Complex) return Complex is
188 return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
191 function "+" (Right : Imaginary) return Imaginary is
196 function "+" (Left, Right : Imaginary) return Imaginary is
198 return Imaginary (R (Left) + R (Right));
201 function "+" (Left : Complex; Right : Real'Base) return Complex is
203 return Complex'(Left.Re + Right, Left.Im);
206 function "+" (Left : Real'Base; Right : Complex) return Complex is
208 return Complex'(Left + Right.Re, Right.Im);
211 function "+" (Left : Complex; Right : Imaginary) return Complex is
213 return Complex'(Left.Re, Left.Im + R (Right));
216 function "+" (Left : Imaginary; Right : Complex) return Complex is
218 return Complex'(Right.Re, R (Left) + Right.Im);
221 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
223 return Complex'(Right, R (Left));
226 function "+" (Left : Real'Base; Right : Imaginary) return Complex is
228 return Complex'(Left, R (Right));
235 function "-" (Right : Complex) return Complex is
237 return (-Right.Re, -Right.Im);
240 function "-" (Left, Right : Complex) return Complex is
242 return (Left.Re - Right.Re, Left.Im - Right.Im);
245 function "-" (Right : Imaginary) return Imaginary is
247 return Imaginary (-R (Right));
250 function "-" (Left, Right : Imaginary) return Imaginary is
252 return Imaginary (R (Left) - R (Right));
255 function "-" (Left : Complex; Right : Real'Base) return Complex is
257 return Complex'(Left.Re - Right, Left.Im);
260 function "-" (Left : Real'Base; Right : Complex) return Complex is
262 return Complex'(Left - Right.Re, -Right.Im);
265 function "-" (Left : Complex; Right : Imaginary) return Complex is
267 return Complex'(Left.Re, Left.Im - R (Right));
270 function "-" (Left : Imaginary; Right : Complex) return Complex is
272 return Complex'(-Right.Re, R (Left) - Right.Im);
275 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
277 return Complex'(-Right, R (Left));
280 function "-" (Left : Real'Base; Right : Imaginary) return Complex is
282 return Complex'(Left, -R (Right));
289 function "/" (Left, Right : Complex) return Complex is
290 a : constant R := Left.Re;
291 b : constant R := Left.Im;
292 c : constant R := Right.Re;
293 d : constant R := Right.Im;
296 if c = 0.0 and then d = 0.0 then
297 raise Constraint_Error;
299 return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
300 Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
304 function "/" (Left, Right : Imaginary) return Real'Base is
306 return R (Left) / R (Right);
309 function "/" (Left : Complex; Right : Real'Base) return Complex is
311 return Complex'(Left.Re / Right, Left.Im / Right);
314 function "/" (Left : Real'Base; Right : Complex) return Complex is
315 a : constant R := Left;
316 c : constant R := Right.Re;
317 d : constant R := Right.Im;
319 return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
320 Im => -(a * d) / (c ** 2 + d ** 2));
323 function "/" (Left : Complex; Right : Imaginary) return Complex is
324 a : constant R := Left.Re;
325 b : constant R := Left.Im;
326 d : constant R := R (Right);
329 return (b / d, -a / d);
332 function "/" (Left : Imaginary; Right : Complex) return Complex is
333 b : constant R := R (Left);
334 c : constant R := Right.Re;
335 d : constant R := Right.Im;
338 return (Re => b * d / (c ** 2 + d ** 2),
339 Im => b * c / (c ** 2 + d ** 2));
342 function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
344 return Imaginary (R (Left) / Right);
347 function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
349 return Imaginary (-Left / R (Right));
356 function "<" (Left, Right : Imaginary) return Boolean is
358 return R (Left) < R (Right);
365 function "<=" (Left, Right : Imaginary) return Boolean is
367 return R (Left) <= R (Right);
374 function ">" (Left, Right : Imaginary) return Boolean is
376 return R (Left) > R (Right);
383 function ">=" (Left, Right : Imaginary) return Boolean is
385 return R (Left) >= R (Right);
392 function "abs" (Right : Imaginary) return Real'Base is
394 return abs R (Right);
401 function Argument (X : Complex) return Real'Base is
402 a : constant R := X.Re;
403 b : constant R := X.Im;
412 return R'Copy_Sign (Pi, b);
424 arg := R (Atan (Double (abs (b / a))));
443 when Constraint_Error =>
451 function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
454 return Argument (X) * Cycle / Two_Pi;
456 raise Argument_Error;
460 ----------------------------
461 -- Compose_From_Cartesian --
462 ----------------------------
464 function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
467 end Compose_From_Cartesian;
469 function Compose_From_Cartesian (Re : Real'Base) return Complex is
472 end Compose_From_Cartesian;
474 function Compose_From_Cartesian (Im : Imaginary) return Complex is
476 return (0.0, R (Im));
477 end Compose_From_Cartesian;
479 ------------------------
480 -- Compose_From_Polar --
481 ------------------------
483 function Compose_From_Polar (
484 Modulus, Argument : Real'Base)
488 if Modulus = 0.0 then
491 return (Modulus * R (Cos (Double (Argument))),
492 Modulus * R (Sin (Double (Argument))));
494 end Compose_From_Polar;
496 function Compose_From_Polar (
497 Modulus, Argument, Cycle : Real'Base)
503 if Modulus = 0.0 then
506 elsif Cycle > 0.0 then
507 if Argument = 0.0 then
508 return (Modulus, 0.0);
510 elsif Argument = Cycle / 4.0 then
511 return (0.0, Modulus);
513 elsif Argument = Cycle / 2.0 then
514 return (-Modulus, 0.0);
516 elsif Argument = 3.0 * Cycle / R (4.0) then
517 return (0.0, -Modulus);
519 Arg := Two_Pi * Argument / Cycle;
520 return (Modulus * R (Cos (Double (Arg))),
521 Modulus * R (Sin (Double (Arg))));
524 raise Argument_Error;
526 end Compose_From_Polar;
532 function Conjugate (X : Complex) return Complex is
534 return Complex'(X.Re, -X.Im);
541 function Im (X : Complex) return Real'Base is
546 function Im (X : Imaginary) return Real'Base is
555 function Modulus (X : Complex) return Real'Base is
563 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
564 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
565 -- squaring does not raise constraint_error but generates infinity,
566 -- we can use an explicit comparison to determine whether to use
567 -- the scaling expression.
569 -- The scaling expression is computed in double format throughout
570 -- in order to prevent inaccuracies on machines where not all
571 -- immediate expressions are rounded, such as PowerPC.
574 raise Constraint_Error;
578 when Constraint_Error =>
579 return R (Double (abs (X.Re))
580 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
587 raise Constraint_Error;
591 when Constraint_Error =>
592 return R (Double (abs (X.Im))
593 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
596 -- Now deal with cases of underflow. If only one of the squares
597 -- underflows, return the modulus of the other component. If both
598 -- squares underflow, use scaling as above.
611 if abs (X.Re) > abs (X.Im) then
613 R (Double (abs (X.Re))
614 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
617 R (Double (abs (X.Im))
618 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
629 -- in all other cases, the naive computation will do.
632 return R (Sqrt (Double (Re2 + Im2)));
640 function Re (X : Complex) return Real'Base is
649 procedure Set_Im (X : in out Complex; Im : in Real'Base) is
654 procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
663 procedure Set_Re (X : in out Complex; Re : in Real'Base) is
668 end Ada.Numerics.Generic_Complex_Types;