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-2009, 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
50 X := Left.Re * Right.Re - Left.Im * Right.Im;
51 Y := Left.Re * Right.Im + Left.Im * Right.Re;
53 -- If either component overflows, try to scale (skip in fast math mode)
55 if not Standard'Fast_Math then
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));
70 function "*" (Left, Right : Imaginary) return Real'Base is
72 return -(R (Left) * R (Right));
75 function "*" (Left : Complex; Right : Real'Base) return Complex is
77 return Complex'(Left.Re * Right, Left.Im * Right);
80 function "*" (Left : Real'Base; Right : Complex) return Complex is
82 return (Left * Right.Re, Left * Right.Im);
85 function "*" (Left : Complex; Right : Imaginary) return Complex is
87 return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
90 function "*" (Left : Imaginary; Right : Complex) return Complex is
92 return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
95 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
97 return Left * Imaginary (Right);
100 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
102 return Imaginary (Left * R (Right));
109 function "**" (Left : Complex; Right : Integer) return Complex is
110 Result : Complex := (1.0, 0.0);
111 Factor : Complex := Left;
112 Exp : Integer := Right;
115 -- We use the standard logarithmic approach, Exp gets shifted right
116 -- testing successive low order bits and Factor is the value of the
117 -- base raised to the next power of 2. For positive exponents we
118 -- multiply the result by this factor, for negative exponents, we
119 -- divide by this factor.
123 -- For a positive exponent, if we get a constraint error during
124 -- this loop, it is an overflow, and the constraint error will
125 -- simply be passed on to the caller.
128 if Exp rem 2 /= 0 then
129 Result := Result * Factor;
132 Factor := Factor * Factor;
140 -- For the negative exponent case, a constraint error during this
141 -- calculation happens if Factor gets too large, and the proper
142 -- response is to return 0.0, since what we essentially have is
143 -- 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;
158 when Constraint_Error =>
164 function "**" (Left : Imaginary; Right : Integer) return Complex is
165 M : constant R := R (Left) ** Right;
168 when 0 => return (M, 0.0);
169 when 1 => return (0.0, M);
170 when 2 => return (-M, 0.0);
171 when 3 => return (0.0, -M);
172 when others => raise Program_Error;
180 function "+" (Right : Complex) return Complex is
185 function "+" (Left, Right : Complex) return Complex is
187 return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
190 function "+" (Right : Imaginary) return Imaginary is
195 function "+" (Left, Right : Imaginary) return Imaginary is
197 return Imaginary (R (Left) + R (Right));
200 function "+" (Left : Complex; Right : Real'Base) return Complex is
202 return Complex'(Left.Re + Right, Left.Im);
205 function "+" (Left : Real'Base; Right : Complex) return Complex is
207 return Complex'(Left + Right.Re, Right.Im);
210 function "+" (Left : Complex; Right : Imaginary) return Complex is
212 return Complex'(Left.Re, Left.Im + R (Right));
215 function "+" (Left : Imaginary; Right : Complex) return Complex is
217 return Complex'(Right.Re, R (Left) + Right.Im);
220 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
222 return Complex'(Right, R (Left));
225 function "+" (Left : Real'Base; Right : Imaginary) return Complex is
227 return Complex'(Left, R (Right));
234 function "-" (Right : Complex) return Complex is
236 return (-Right.Re, -Right.Im);
239 function "-" (Left, Right : Complex) return Complex is
241 return (Left.Re - Right.Re, Left.Im - Right.Im);
244 function "-" (Right : Imaginary) return Imaginary is
246 return Imaginary (-R (Right));
249 function "-" (Left, Right : Imaginary) return Imaginary is
251 return Imaginary (R (Left) - R (Right));
254 function "-" (Left : Complex; Right : Real'Base) return Complex is
256 return Complex'(Left.Re - Right, Left.Im);
259 function "-" (Left : Real'Base; Right : Complex) return Complex is
261 return Complex'(Left - Right.Re, -Right.Im);
264 function "-" (Left : Complex; Right : Imaginary) return Complex is
266 return Complex'(Left.Re, Left.Im - R (Right));
269 function "-" (Left : Imaginary; Right : Complex) return Complex is
271 return Complex'(-Right.Re, R (Left) - Right.Im);
274 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
276 return Complex'(-Right, R (Left));
279 function "-" (Left : Real'Base; Right : Imaginary) return Complex is
281 return Complex'(Left, -R (Right));
288 function "/" (Left, Right : Complex) return Complex is
289 a : constant R := Left.Re;
290 b : constant R := Left.Im;
291 c : constant R := Right.Re;
292 d : constant R := Right.Im;
295 if c = 0.0 and then d = 0.0 then
296 raise Constraint_Error;
298 return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
299 Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
303 function "/" (Left, Right : Imaginary) return Real'Base is
305 return R (Left) / R (Right);
308 function "/" (Left : Complex; Right : Real'Base) return Complex is
310 return Complex'(Left.Re / Right, Left.Im / Right);
313 function "/" (Left : Real'Base; Right : Complex) return Complex is
314 a : constant R := Left;
315 c : constant R := Right.Re;
316 d : constant R := Right.Im;
318 return Complex'(Re => (a * c) / (c ** 2 + d ** 2),
319 Im => -((a * d) / (c ** 2 + d ** 2)));
322 function "/" (Left : Complex; Right : Imaginary) return Complex is
323 a : constant R := Left.Re;
324 b : constant R := Left.Im;
325 d : constant R := R (Right);
328 return (b / d, -(a / d));
331 function "/" (Left : Imaginary; Right : Complex) return Complex is
332 b : constant R := R (Left);
333 c : constant R := Right.Re;
334 d : constant R := Right.Im;
337 return (Re => b * d / (c ** 2 + d ** 2),
338 Im => b * c / (c ** 2 + d ** 2));
341 function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
343 return Imaginary (R (Left) / Right);
346 function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
348 return Imaginary (-(Left / R (Right)));
355 function "<" (Left, Right : Imaginary) return Boolean is
357 return R (Left) < R (Right);
364 function "<=" (Left, Right : Imaginary) return Boolean is
366 return R (Left) <= R (Right);
373 function ">" (Left, Right : Imaginary) return Boolean is
375 return R (Left) > R (Right);
382 function ">=" (Left, Right : Imaginary) return Boolean is
384 return R (Left) >= R (Right);
391 function "abs" (Right : Imaginary) return Real'Base is
393 return abs R (Right);
400 function Argument (X : Complex) return Real'Base is
401 a : constant R := X.Re;
402 b : constant R := X.Im;
411 return R'Copy_Sign (Pi, b);
423 arg := R (Atan (Double (abs (b / a))));
442 when Constraint_Error =>
450 function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
453 return Argument (X) * Cycle / Two_Pi;
455 raise Argument_Error;
459 ----------------------------
460 -- Compose_From_Cartesian --
461 ----------------------------
463 function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
466 end Compose_From_Cartesian;
468 function Compose_From_Cartesian (Re : Real'Base) return Complex is
471 end Compose_From_Cartesian;
473 function Compose_From_Cartesian (Im : Imaginary) return Complex is
475 return (0.0, R (Im));
476 end Compose_From_Cartesian;
478 ------------------------
479 -- Compose_From_Polar --
480 ------------------------
482 function Compose_From_Polar (
483 Modulus, Argument : Real'Base)
487 if Modulus = 0.0 then
490 return (Modulus * R (Cos (Double (Argument))),
491 Modulus * R (Sin (Double (Argument))));
493 end Compose_From_Polar;
495 function Compose_From_Polar (
496 Modulus, Argument, Cycle : Real'Base)
502 if Modulus = 0.0 then
505 elsif Cycle > 0.0 then
506 if Argument = 0.0 then
507 return (Modulus, 0.0);
509 elsif Argument = Cycle / 4.0 then
510 return (0.0, Modulus);
512 elsif Argument = Cycle / 2.0 then
513 return (-Modulus, 0.0);
515 elsif Argument = 3.0 * Cycle / R (4.0) then
516 return (0.0, -Modulus);
518 Arg := Two_Pi * Argument / Cycle;
519 return (Modulus * R (Cos (Double (Arg))),
520 Modulus * R (Sin (Double (Arg))));
523 raise Argument_Error;
525 end Compose_From_Polar;
531 function Conjugate (X : Complex) return Complex is
533 return Complex'(X.Re, -X.Im);
540 function Im (X : Complex) return Real'Base is
545 function Im (X : Imaginary) return Real'Base is
554 function Modulus (X : Complex) return Real'Base is
562 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
563 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
564 -- squaring does not raise constraint_error but generates infinity,
565 -- we can use an explicit comparison to determine whether to use
566 -- the scaling expression.
568 -- The scaling expression is computed in double format throughout
569 -- in order to prevent inaccuracies on machines where not all
570 -- immediate expressions are rounded, such as PowerPC.
573 raise Constraint_Error;
577 when Constraint_Error =>
578 return R (Double (abs (X.Re))
579 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
586 raise Constraint_Error;
590 when Constraint_Error =>
591 return R (Double (abs (X.Im))
592 * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
595 -- Now deal with cases of underflow. If only one of the squares
596 -- underflows, return the modulus of the other component. If both
597 -- squares underflow, use scaling as above.
610 if abs (X.Re) > abs (X.Im) then
612 R (Double (abs (X.Re))
613 * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
616 R (Double (abs (X.Im))
617 * Sqrt (1.0 + (Double (X.Re) / Double (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 : Real'Base) is
653 procedure Set_Im (X : out Imaginary; Im : Real'Base) is
662 procedure Set_Re (X : in out Complex; Re : Real'Base) is
667 end Ada.Numerics.Generic_Complex_Types;