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-2024, 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_Generic_Float
;
34 package body Ada
.Numerics
.Generic_Complex_Types
is
36 package Aux
is new Ada
.Numerics
.Aux_Generic_Float
(Real
);
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
49 Scale
: constant R
:= R
(R
'Machine_Radix) ** ((R
'Machine_Emax - 1) / 2);
50 -- In case of overflow, scale the operands by the largest power of the
51 -- radix (to avoid rounding error), so that the square of the scale does
52 -- not overflow itself.
58 X
:= Left
.Re
* Right
.Re
- Left
.Im
* Right
.Im
;
59 Y
:= Left
.Re
* Right
.Im
+ Left
.Im
* Right
.Re
;
61 -- If either component overflows, try to scale (skip in fast math mode)
63 if not Standard
'Fast_Math then
65 -- Note that the test below is written as a negation. This is to
66 -- account for the fact that X and Y may be NaNs, because both of
67 -- their operands could overflow. Given that all operations on NaNs
68 -- return false, the test can only be written thus.
70 if not (abs (X
) <= R
'Last) then
72 (CodePeer
, Intentional
,
73 "test always false", "test for infinity");
75 X
:= Scale
**2 * ((Left
.Re
/ Scale
) * (Right
.Re
/ Scale
) -
76 (Left
.Im
/ Scale
) * (Right
.Im
/ Scale
));
79 if not (abs (Y
) <= R
'Last) then
81 (CodePeer
, Intentional
,
82 "test always false", "test for infinity");
84 Y
:= Scale
**2 * ((Left
.Re
/ Scale
) * (Right
.Im
/ Scale
)
85 + (Left
.Im
/ Scale
) * (Right
.Re
/ Scale
));
92 function "*" (Left
, Right
: Imaginary
) return Real
'Base is
94 return -(R
(Left
) * R
(Right
));
97 function "*" (Left
: Complex
; Right
: Real
'Base) return Complex
is
99 return Complex
'(Left.Re * Right, Left.Im * Right);
102 function "*" (Left : Real'Base; Right : Complex) return Complex is
104 return (Left * Right.Re, Left * Right.Im);
107 function "*" (Left : Complex; Right : Imaginary) return Complex is
109 return Complex'(-(Left
.Im
* R
(Right
)), Left
.Re
* R
(Right
));
112 function "*" (Left
: Imaginary
; Right
: Complex
) return Complex
is
114 return Complex
'(-(R (Left) * Right.Im), R (Left) * Right.Re);
117 function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
119 return Left * Imaginary (Right);
122 function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
124 return Imaginary (Left * R (Right));
131 function "**" (Left : Complex; Right : Integer) return Complex is
132 Result : Complex := (1.0, 0.0);
133 Factor : Complex := Left;
134 Exp : Integer := Right;
137 -- We use the standard logarithmic approach, Exp gets shifted right
138 -- testing successive low order bits and Factor is the value of the
139 -- base raised to the next power of 2. For positive exponents we
140 -- multiply the result by this factor, for negative exponents, we
141 -- divide by this factor.
145 -- For a positive exponent, if we get a constraint error during
146 -- this loop, it is an overflow, and the constraint error will
147 -- simply be passed on to the caller.
150 if Exp rem 2 /= 0 then
151 Result := Result * Factor;
154 Factor := Factor * Factor;
162 -- For the negative exponent case, a constraint error during this
163 -- calculation happens if Factor gets too large, and the proper
164 -- response is to return 0.0, since what we essentially have is
165 -- 1.0 / infinity, and the closest model number will be zero.
169 if Exp rem 2 /= 0 then
170 Result := Result * Factor;
173 Factor := Factor * Factor;
177 return R'(1.0) / Result
;
180 when Constraint_Error
=>
186 function "**" (Left
: Imaginary
; Right
: Integer) return Complex
is
187 M
: constant R
:= R
(Left
) ** Right
;
190 when 0 => return (M
, 0.0);
191 when 1 => return (0.0, M
);
192 when 2 => return (-M
, 0.0);
193 when 3 => return (0.0, -M
);
194 when others => raise Program_Error
;
202 function "+" (Right
: Complex
) return Complex
is
207 function "+" (Left
, Right
: Complex
) return Complex
is
209 return Complex
'(Left.Re + Right.Re, Left.Im + Right.Im);
212 function "+" (Right : Imaginary) return Imaginary is
217 function "+" (Left, Right : Imaginary) return Imaginary is
219 return Imaginary (R (Left) + R (Right));
222 function "+" (Left : Complex; Right : Real'Base) return Complex is
224 return Complex'(Left
.Re
+ Right
, Left
.Im
);
227 function "+" (Left
: Real
'Base; Right
: Complex
) return Complex
is
229 return Complex
'(Left + Right.Re, Right.Im);
232 function "+" (Left : Complex; Right : Imaginary) return Complex is
234 return Complex'(Left
.Re
, Left
.Im
+ R
(Right
));
237 function "+" (Left
: Imaginary
; Right
: Complex
) return Complex
is
239 return Complex
'(Right.Re, R (Left) + Right.Im);
242 function "+" (Left : Imaginary; Right : Real'Base) return Complex is
244 return Complex'(Right
, R
(Left
));
247 function "+" (Left
: Real
'Base; Right
: Imaginary
) return Complex
is
249 return Complex
'(Left, R (Right));
256 function "-" (Right : Complex) return Complex is
258 return (-Right.Re, -Right.Im);
261 function "-" (Left, Right : Complex) return Complex is
263 return (Left.Re - Right.Re, Left.Im - Right.Im);
266 function "-" (Right : Imaginary) return Imaginary is
268 return Imaginary (-R (Right));
271 function "-" (Left, Right : Imaginary) return Imaginary is
273 return Imaginary (R (Left) - R (Right));
276 function "-" (Left : Complex; Right : Real'Base) return Complex is
278 return Complex'(Left
.Re
- Right
, Left
.Im
);
281 function "-" (Left
: Real
'Base; Right
: Complex
) return Complex
is
283 return Complex
'(Left - Right.Re, -Right.Im);
286 function "-" (Left : Complex; Right : Imaginary) return Complex is
288 return Complex'(Left
.Re
, Left
.Im
- R
(Right
));
291 function "-" (Left
: Imaginary
; Right
: Complex
) return Complex
is
293 return Complex
'(-Right.Re, R (Left) - Right.Im);
296 function "-" (Left : Imaginary; Right : Real'Base) return Complex is
298 return Complex'(-Right
, R
(Left
));
301 function "-" (Left
: Real
'Base; Right
: Imaginary
) return Complex
is
303 return Complex
'(Left, -R (Right));
310 function "/" (Left, Right : Complex) return Complex is
311 a : constant R := Left.Re;
312 b : constant R := Left.Im;
313 c : constant R := Right.Re;
314 d : constant R := Right.Im;
317 if c = 0.0 and then d = 0.0 then
318 raise Constraint_Error;
320 return Complex'(Re
=> ((a
* c
) + (b
* d
)) / (c
** 2 + d
** 2),
321 Im
=> ((b
* c
) - (a
* d
)) / (c
** 2 + d
** 2));
325 function "/" (Left
, Right
: Imaginary
) return Real
'Base is
327 return R
(Left
) / R
(Right
);
330 function "/" (Left
: Complex
; Right
: Real
'Base) return Complex
is
332 return Complex
'(Left.Re / Right, Left.Im / Right);
335 function "/" (Left : Real'Base; Right : Complex) return Complex is
336 a : constant R := Left;
337 c : constant R := Right.Re;
338 d : constant R := Right.Im;
340 return Complex'(Re
=> (a
* c
) / (c
** 2 + d
** 2),
341 Im
=> -((a
* d
) / (c
** 2 + d
** 2)));
344 function "/" (Left
: Complex
; Right
: Imaginary
) return Complex
is
345 a
: constant R
:= Left
.Re
;
346 b
: constant R
:= Left
.Im
;
347 d
: constant R
:= R
(Right
);
350 return (b
/ d
, -(a
/ d
));
353 function "/" (Left
: Imaginary
; Right
: Complex
) return Complex
is
354 b
: constant R
:= R
(Left
);
355 c
: constant R
:= Right
.Re
;
356 d
: constant R
:= Right
.Im
;
359 return (Re
=> b
* d
/ (c
** 2 + d
** 2),
360 Im
=> b
* c
/ (c
** 2 + d
** 2));
363 function "/" (Left
: Imaginary
; Right
: Real
'Base) return Imaginary
is
365 return Imaginary
(R
(Left
) / Right
);
368 function "/" (Left
: Real
'Base; Right
: Imaginary
) return Imaginary
is
370 return Imaginary
(-(Left
/ R
(Right
)));
377 function "<" (Left
, Right
: Imaginary
) return Boolean is
379 return R
(Left
) < R
(Right
);
386 function "<=" (Left
, Right
: Imaginary
) return Boolean is
388 return R
(Left
) <= R
(Right
);
395 function ">" (Left
, Right
: Imaginary
) return Boolean is
397 return R
(Left
) > R
(Right
);
404 function ">=" (Left
, Right
: Imaginary
) return Boolean is
406 return R
(Left
) >= R
(Right
);
413 function "abs" (Right
: Imaginary
) return Real
'Base is
415 return abs R
(Right
);
422 function Argument
(X
: Complex
) return Real
'Base is
423 a
: constant R
:= X
.Re
;
424 b
: constant R
:= X
.Im
;
433 return R
'Copy_Sign (Pi
, b
);
445 arg
:= Aux
.Atan
(abs (b
/ a
));
464 when Constraint_Error
=>
472 function Argument
(X
: Complex
; Cycle
: Real
'Base) return Real
'Base is
475 return Argument
(X
) * Cycle
/ Two_Pi
;
477 raise Argument_Error
;
481 ----------------------------
482 -- Compose_From_Cartesian --
483 ----------------------------
485 function Compose_From_Cartesian
(Re
, Im
: Real
'Base) return Complex
is
488 end Compose_From_Cartesian
;
490 function Compose_From_Cartesian
(Re
: Real
'Base) return Complex
is
493 end Compose_From_Cartesian
;
495 function Compose_From_Cartesian
(Im
: Imaginary
) return Complex
is
497 return (0.0, R
(Im
));
498 end Compose_From_Cartesian
;
500 ------------------------
501 -- Compose_From_Polar --
502 ------------------------
504 function Compose_From_Polar
(
505 Modulus
, Argument
: Real
'Base)
509 if Modulus
= 0.0 then
512 return (Modulus
* Aux
.Cos
(Argument
),
513 Modulus
* Aux
.Sin
(Argument
));
515 end Compose_From_Polar
;
517 function Compose_From_Polar
(
518 Modulus
, Argument
, Cycle
: Real
'Base)
524 if Modulus
= 0.0 then
527 elsif Cycle
> 0.0 then
528 if Argument
= 0.0 then
529 return (Modulus
, 0.0);
531 elsif Argument
= Cycle
/ 4.0 then
532 return (0.0, Modulus
);
534 elsif Argument
= Cycle
/ 2.0 then
535 return (-Modulus
, 0.0);
537 elsif Argument
= 3.0 * Cycle
/ R
(4.0) then
538 return (0.0, -Modulus
);
540 Arg
:= Two_Pi
* Argument
/ Cycle
;
541 return (Modulus
* Aux
.Cos
(Arg
),
542 Modulus
* Aux
.Sin
(Arg
));
545 raise Argument_Error
;
547 end Compose_From_Polar
;
553 function Conjugate
(X
: Complex
) return Complex
is
555 return Complex
'(X.Re, -X.Im);
562 function Im (X : Complex) return Real'Base is
567 function Im (X : Imaginary) return Real'Base is
576 function Modulus (X : Complex) return Real'Base is
584 -- To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
585 -- compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
586 -- squaring does not raise constraint_error but generates infinity,
587 -- we can use an explicit comparison to determine whether to use
588 -- the scaling expression.
590 -- The scaling expression is computed in double format throughout
591 -- in order to prevent inaccuracies on machines where not all
592 -- immediate expressions are rounded, such as PowerPC.
594 -- ??? same weird test, why not Re2 > R'Last ???
595 if not (Re2 <= R'Last) then
596 raise Constraint_Error;
600 when Constraint_Error =>
601 pragma Assert (X.Re /= 0.0);
602 return R (abs (X.Re))
603 * Aux.Sqrt (1.0 + (R (X.Im) / R (X.Re)) ** 2);
609 -- ??? same weird test
610 if not (Im2 <= R'Last) then
611 raise Constraint_Error;
615 when Constraint_Error =>
616 pragma Assert (X.Im /= 0.0);
617 return R (abs (X.Im))
618 * Aux.Sqrt (1.0 + (R (X.Re) / R (X.Im)) ** 2);
621 -- Now deal with cases of underflow. If only one of the squares
622 -- underflows, return the modulus of the other component. If both
623 -- squares underflow, use scaling as above.
636 if abs (X.Re) > abs (X.Im) then
637 return R (abs (X.Re))
638 * Aux.Sqrt (1.0 + (R (X.Im) / R (X.Re)) ** 2);
640 return R (abs (X.Im))
641 * Aux.Sqrt (1.0 + (R (X.Re) / R (X.Im)) ** 2);
652 -- In all other cases, the naive computation will do
655 return Aux.Sqrt (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;