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 --
10 -- Copyright (C) 1992-2001 Free Software Foundation, Inc. --
12 -- GNAT is free software; you can redistribute it and/or modify it under --
13 -- terms of the GNU General Public License as published by the Free Soft- --
14 -- ware Foundation; either version 2, or (at your option) any later ver- --
15 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
16 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
17 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
18 -- for more details. You should have received a copy of the GNU General --
19 -- Public License distributed with GNAT; see file COPYING. If not, write --
20 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
21 -- MA 02111-1307, USA. --
23 -- As a special exception, if other files instantiate generics from this --
24 -- unit, or you link this unit with other files to produce an executable, --
25 -- this unit does not by itself cause the resulting executable to be --
26 -- covered by the GNU General Public License. This exception does not --
27 -- however invalidate any other reasons why the executable file might be --
28 -- covered by the GNU Public License. --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- Extensive contributions were provided by Ada Core Technologies Inc. --
33 ------------------------------------------------------------------------------
35 with Ada
.Numerics
.Aux
; use Ada
.Numerics
.Aux
;
36 package body Ada
.Numerics
.Generic_Complex_Types
is
38 subtype R
is Real
'Base;
40 Two_Pi
: constant R
:= R
(2.0) * Pi
;
41 Half_Pi
: constant R
:= Pi
/ R
(2.0);
47 function "*" (Left
, Right
: Complex
) return Complex
is
52 X
:= Left
.Re
* Right
.Re
- Left
.Im
* Right
.Im
;
53 Y
:= Left
.Re
* Right
.Im
+ Left
.Im
* Right
.Re
;
55 -- If either component overflows, try to scale.
57 if abs (X
) > R
'Last then
58 X
:= R
' (4.0) * (R'(Left
.Re
/ 2.0) * R
'(Right.Re / 2.0)
59 - R'(Left
.Im
/ 2.0) * R
'(Right.Im / 2.0));
62 if abs (Y) > R'Last then
63 Y := R' (4.0) * (R
'(Left.Re / 2.0) * R'(Right
.Im
/ 2.0)
64 - 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.
148 if Exp rem 2 /= 0 then
149 Result := Result * Factor;
152 Factor := Factor * Factor;
156 return R ' (1.0) / Result
;
160 when Constraint_Error
=>
166 function "**" (Left
: Imaginary
; Right
: Integer) return Complex
is
167 M
: R
:= R
(Left
) ** Right
;
170 when 0 => return (M
, 0.0);
171 when 1 => return (0.0, M
);
172 when 2 => return (-M
, 0.0);
173 when 3 => return (0.0, -M
);
174 when others => raise Program_Error
;
182 function "+" (Right
: Complex
) return Complex
is
187 function "+" (Left
, Right
: Complex
) return Complex
is
189 return Complex
'(Left.Re + Right.Re, Left.Im + Right.Im);
192 function "+" (Right : Imaginary) return Imaginary is
197 function "+" (Left, Right : Imaginary) return Imaginary is
199 return Imaginary (R (Left) + R (Right));
202 function "+" (Left : Complex; Right : Real'Base) return Complex is
204 return Complex'(Left
.Re
+ Right
, Left
.Im
);
207 function "+" (Left
: Real
'Base; Right
: Complex
) return Complex
is
209 return Complex
'(Left + Right.Re, Right.Im);
212 function "+" (Left : Complex; Right : Imaginary) return Complex is
214 return Complex'(Left
.Re
, Left
.Im
+ R
(Right
));
217 function "+" (Left
: Imaginary
; Right
: Complex
) return Complex
is
219 return Complex
'(Right.Re, R (Left) + Right.Im);
222 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
224 return Complex'(Right
, R
(Left
));
227 function "+" (Left
: Real
'Base; Right
: Imaginary
) return Complex
is
229 return Complex
'(Left, R (Right));
236 function "-" (Right : Complex) return Complex is
238 return (-Right.Re, -Right.Im);
241 function "-" (Left, Right : Complex) return Complex is
243 return (Left.Re - Right.Re, Left.Im - Right.Im);
246 function "-" (Right : Imaginary) return Imaginary is
248 return Imaginary (-R (Right));
251 function "-" (Left, Right : Imaginary) return Imaginary is
253 return Imaginary (R (Left) - R (Right));
256 function "-" (Left : Complex; Right : Real'Base) return Complex is
258 return Complex'(Left
.Re
- Right
, Left
.Im
);
261 function "-" (Left
: Real
'Base; Right
: Complex
) return Complex
is
263 return Complex
'(Left - Right.Re, -Right.Im);
266 function "-" (Left : Complex; Right : Imaginary) return Complex is
268 return Complex'(Left
.Re
, Left
.Im
- R
(Right
));
271 function "-" (Left
: Imaginary
; Right
: Complex
) return Complex
is
273 return Complex
'(-Right.Re, R (Left) - Right.Im);
276 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
278 return Complex'(-Right
, R
(Left
));
281 function "-" (Left
: Real
'Base; Right
: Imaginary
) return Complex
is
283 return Complex
'(Left, -R (Right));
290 function "/" (Left, Right : Complex) return Complex is
291 a : constant R := Left.Re;
292 b : constant R := Left.Im;
293 c : constant R := Right.Re;
294 d : constant R := Right.Im;
297 if c = 0.0 and then d = 0.0 then
298 raise Constraint_Error;
300 return Complex'(Re
=> ((a
* c
) + (b
* d
)) / (c
** 2 + d
** 2),
301 Im
=> ((b
* c
) - (a
* d
)) / (c
** 2 + d
** 2));
305 function "/" (Left
, Right
: Imaginary
) return Real
'Base is
307 return R
(Left
) / R
(Right
);
310 function "/" (Left
: Complex
; Right
: Real
'Base) return Complex
is
312 return Complex
'(Left.Re / Right, Left.Im / Right);
315 function "/" (Left : Real'Base; Right : Complex) return Complex is
316 a : constant R := Left;
317 c : constant R := Right.Re;
318 d : constant R := Right.Im;
320 return Complex'(Re
=> (a
* c
) / (c
** 2 + d
** 2),
321 Im
=> -(a
* d
) / (c
** 2 + d
** 2));
324 function "/" (Left
: Complex
; Right
: Imaginary
) return Complex
is
325 a
: constant R
:= Left
.Re
;
326 b
: constant R
:= Left
.Im
;
327 d
: constant R
:= R
(Right
);
330 return (b
/ d
, -a
/ d
);
333 function "/" (Left
: Imaginary
; Right
: Complex
) return Complex
is
334 b
: constant R
:= R
(Left
);
335 c
: constant R
:= Right
.Re
;
336 d
: constant R
:= Right
.Im
;
339 return (Re
=> b
* d
/ (c
** 2 + d
** 2),
340 Im
=> b
* c
/ (c
** 2 + d
** 2));
343 function "/" (Left
: Imaginary
; Right
: Real
'Base) return Imaginary
is
345 return Imaginary
(R
(Left
) / Right
);
348 function "/" (Left
: Real
'Base; Right
: Imaginary
) return Imaginary
is
350 return Imaginary
(-Left
/ R
(Right
));
357 function "<" (Left
, Right
: Imaginary
) return Boolean is
359 return R
(Left
) < R
(Right
);
366 function "<=" (Left
, Right
: Imaginary
) return Boolean is
368 return R
(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 "abs" (Right
: Imaginary
) return Real
'Base is
395 return abs R
(Right
);
402 function Argument
(X
: Complex
) return Real
'Base is
403 a
: constant R
:= X
.Re
;
404 b
: constant R
:= X
.Im
;
413 return R
'Copy_Sign (Pi
, b
);
425 arg
:= R
(Atan
(Double
(abs (b
/ a
))));
444 when Constraint_Error
=>
452 function Argument
(X
: Complex
; Cycle
: Real
'Base) return Real
'Base is
455 return Argument
(X
) * Cycle
/ Two_Pi
;
457 raise Argument_Error
;
461 ----------------------------
462 -- Compose_From_Cartesian --
463 ----------------------------
465 function Compose_From_Cartesian
(Re
, Im
: Real
'Base) return Complex
is
468 end Compose_From_Cartesian
;
470 function Compose_From_Cartesian
(Re
: Real
'Base) return Complex
is
473 end Compose_From_Cartesian
;
475 function Compose_From_Cartesian
(Im
: Imaginary
) return Complex
is
477 return (0.0, R
(Im
));
478 end Compose_From_Cartesian
;
480 ------------------------
481 -- Compose_From_Polar --
482 ------------------------
484 function Compose_From_Polar
(
485 Modulus
, Argument
: Real
'Base)
489 if Modulus
= 0.0 then
492 return (Modulus
* R
(Cos
(Double
(Argument
))),
493 Modulus
* R
(Sin
(Double
(Argument
))));
495 end Compose_From_Polar
;
497 function Compose_From_Polar
(
498 Modulus
, Argument
, Cycle
: Real
'Base)
504 if Modulus
= 0.0 then
507 elsif Cycle
> 0.0 then
508 if Argument
= 0.0 then
509 return (Modulus
, 0.0);
511 elsif Argument
= Cycle
/ 4.0 then
512 return (0.0, Modulus
);
514 elsif Argument
= Cycle
/ 2.0 then
515 return (-Modulus
, 0.0);
517 elsif Argument
= 3.0 * Cycle
/ R
(4.0) then
518 return (0.0, -Modulus
);
520 Arg
:= Two_Pi
* Argument
/ Cycle
;
521 return (Modulus
* R
(Cos
(Double
(Arg
))),
522 Modulus
* R
(Sin
(Double
(Arg
))));
525 raise Argument_Error
;
527 end Compose_From_Polar
;
533 function Conjugate
(X
: Complex
) return Complex
is
535 return Complex
'(X.Re, -X.Im);
542 function Im (X : Complex) return Real'Base is
547 function Im (X : Imaginary) return Real'Base is
556 function Modulus (X : Complex) return Real'Base is
564 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
565 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
566 -- squaring does not raise constraint_error but generates infinity,
567 -- we can use an explicit comparison to determine whether to use
568 -- the scaling expression.
571 raise Constraint_Error;
575 when Constraint_Error =>
577 * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
584 raise Constraint_Error;
588 when Constraint_Error =>
590 * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
593 -- Now deal with cases of underflow. If only one of the squares
594 -- underflows, return the modulus of the other component. If both
595 -- squares underflow, use scaling as above.
608 if abs (X.Re) > abs (X.Im) then
611 * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
615 * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
627 -- in all other cases, the naive computation will do.
630 return R (Sqrt (Double (Re2 + Im2)));
638 function Re (X : Complex) return Real'Base is
647 procedure Set_Im (X : in out Complex; Im : in Real'Base) is
652 procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
661 procedure Set_Re (X : in out Complex; Re : in Real'Base) is
666 end Ada.Numerics.Generic_Complex_Types;