URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [gcc/] [ada/] [a-ngcoty.adb] - Blame information for rev 826

Details | Compare with Previous | View Log


------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--   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    --
--                                                                          --
--                                 B o d y                                  --
--                                                                          --
--          Copyright (C) 1992-2009, Free Software Foundation, Inc.         --
--                                                                          --
-- GNAT is free software;  you can  redistribute it  and/or modify it under --
-- terms of the  GNU General Public License as published  by the Free Soft- --
-- ware  Foundation;  either version 3,  or (at your option) any later ver- --
-- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
-- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE.                                     --
--                                                                          --
-- As a special exception under Section 7 of GPL version 3, you are granted --
-- additional permissions described in the GCC Runtime Library Exception,   --
-- version 3.1, as published by the Free Software Foundation.               --
--                                                                          --
-- You should have received a copy of the GNU General Public License and    --
-- a copy of the GCC Runtime Library Exception along with this program;     --
-- see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see    --
-- <http://www.gnu.org/licenses/>.                                          --
--                                                                          --
-- GNAT was originally developed  by the GNAT team at  New York University. --
-- Extensive contributions were provided by Ada Core Technologies Inc.      --
--                                                                          --
------------------------------------------------------------------------------
 
with Ada.Numerics.Aux; use Ada.Numerics.Aux;
 
package body Ada.Numerics.Generic_Complex_Types is
 
   subtype R is Real'Base;
 
   Two_Pi  : constant R := R (2.0) * Pi;
   Half_Pi : constant R := Pi / R (2.0);
 
   ---------
   -- "*" --
   ---------
 
   function "*" (Left, Right : Complex) return Complex is
      X : R;
      Y : R;
 
   begin
      X := Left.Re * Right.Re - Left.Im * Right.Im;
      Y := Left.Re * Right.Im + Left.Im * Right.Re;
 
      --  If either component overflows, try to scale (skip in fast math mode)
 
      if not Standard'Fast_Math then
         if abs (X) > R'Last then
            X := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Re / 2.0)
                            - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
         end if;
 
         if abs (Y) > R'Last then
            Y := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Im / 2.0)
                            - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
         end if;
      end if;
 
      return (X, Y);
   end "*";
 
   function "*" (Left, Right : Imaginary) return Real'Base is
   begin
      return -(R (Left) * R (Right));
   end "*";
 
   function "*" (Left : Complex; Right : Real'Base) return Complex is
   begin
      return Complex'(Left.Re * Right, Left.Im * Right);
   end "*";
 
   function "*" (Left : Real'Base; Right : Complex) return Complex is
   begin
      return (Left * Right.Re, Left * Right.Im);
   end "*";
 
   function "*" (Left : Complex; Right : Imaginary) return Complex is
   begin
      return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
   end "*";
 
   function "*" (Left : Imaginary; Right : Complex) return Complex is
   begin
      return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
   end "*";
 
   function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
   begin
      return Left * Imaginary (Right);
   end "*";
 
   function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
   begin
      return Imaginary (Left * R (Right));
   end "*";
 
   ----------
   -- "**" --
   ----------
 
   function "**" (Left : Complex; Right : Integer) return Complex is
      Result : Complex := (1.0, 0.0);
      Factor : Complex := Left;
      Exp    : Integer := Right;
 
   begin
      --  We use the standard logarithmic approach, Exp gets shifted right
      --  testing successive low order bits and Factor is the value of the
      --  base raised to the next power of 2. For positive exponents we
      --  multiply the result by this factor, for negative exponents, we
      --  divide by this factor.
 
      if Exp >= 0 then
 
         --  For a positive exponent, if we get a constraint error during
         --  this loop, it is an overflow, and the constraint error will
         --  simply be passed on to the caller.
 
         while Exp /= 0 loop
            if Exp rem 2 /= 0 then
               Result := Result * Factor;
            end if;
 
            Factor := Factor * Factor;
            Exp := Exp / 2;
         end loop;
 
         return Result;
 
      else -- Exp < 0 then
 
         --  For the negative exponent case, a constraint error during this
         --  calculation happens if Factor gets too large, and the proper
         --  response is to return 0.0, since what we essentially have is
         --  1.0 / infinity, and the closest model number will be zero.
 
         begin
            while Exp /= 0 loop
               if Exp rem 2 /= 0 then
                  Result := Result * Factor;
               end if;
 
               Factor := Factor * Factor;
               Exp := Exp / 2;
            end loop;
 
            return R'(1.0) / Result;
 
         exception
            when Constraint_Error =>
               return (0.0, 0.0);
         end;
      end if;
   end "**";
 
   function "**" (Left : Imaginary; Right : Integer) return Complex is
      M : constant R := R (Left) ** Right;
   begin
      case Right mod 4 is
         when 0 => return (M,   0.0);
         when 1 => return (0.0, M);
         when 2 => return (-M,  0.0);
         when 3 => return (0.0, -M);
         when others => raise Program_Error;
      end case;
   end "**";
 
   ---------
   -- "+" --
   ---------
 
   function "+" (Right : Complex) return Complex is
   begin
      return Right;
   end "+";
 
   function "+" (Left, Right : Complex) return Complex is
   begin
      return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
   end "+";
 
   function "+" (Right : Imaginary) return Imaginary is
   begin
      return Right;
   end "+";
 
   function "+" (Left, Right : Imaginary) return Imaginary is
   begin
      return Imaginary (R (Left) + R (Right));
   end "+";
 
   function "+" (Left : Complex; Right : Real'Base) return Complex is
   begin
      return Complex'(Left.Re + Right, Left.Im);
   end "+";
 
   function "+" (Left : Real'Base; Right : Complex) return Complex is
   begin
      return Complex'(Left + Right.Re, Right.Im);
   end "+";
 
   function "+" (Left : Complex; Right : Imaginary) return Complex is
   begin
      return Complex'(Left.Re, Left.Im + R (Right));
   end "+";
 
   function "+" (Left : Imaginary; Right : Complex) return Complex is
   begin
      return Complex'(Right.Re, R (Left) + Right.Im);
   end "+";
 
   function "+" (Left : Imaginary; Right : Real'Base) return Complex is
   begin
      return Complex'(Right, R (Left));
   end "+";
 
   function "+" (Left : Real'Base; Right : Imaginary) return Complex is
   begin
      return Complex'(Left, R (Right));
   end "+";
 
   ---------
   -- "-" --
   ---------
 
   function "-" (Right : Complex) return Complex is
   begin
      return (-Right.Re, -Right.Im);
   end "-";
 
   function "-" (Left, Right : Complex) return Complex is
   begin
      return (Left.Re - Right.Re, Left.Im - Right.Im);
   end "-";
 
   function "-" (Right : Imaginary) return Imaginary is
   begin
      return Imaginary (-R (Right));
   end "-";
 
   function "-" (Left, Right : Imaginary) return Imaginary is
   begin
      return Imaginary (R (Left) - R (Right));
   end "-";
 
   function "-" (Left : Complex; Right : Real'Base) return Complex is
   begin
      return Complex'(Left.Re - Right, Left.Im);
   end "-";
 
   function "-" (Left : Real'Base; Right : Complex) return Complex is
   begin
      return Complex'(Left - Right.Re, -Right.Im);
   end "-";
 
   function "-" (Left : Complex; Right : Imaginary) return Complex is
   begin
      return Complex'(Left.Re, Left.Im - R (Right));
   end "-";
 
   function "-" (Left : Imaginary; Right : Complex) return Complex is
   begin
      return Complex'(-Right.Re, R (Left) - Right.Im);
   end "-";
 
   function "-" (Left : Imaginary; Right : Real'Base) return Complex is
   begin
      return Complex'(-Right, R (Left));
   end "-";
 
   function "-" (Left : Real'Base; Right : Imaginary) return Complex is
   begin
      return Complex'(Left, -R (Right));
   end "-";
 
   ---------
   -- "/" --
   ---------
 
   function "/" (Left, Right : Complex) return Complex is
      a : constant R := Left.Re;
      b : constant R := Left.Im;
      c : constant R := Right.Re;
      d : constant R := Right.Im;
 
   begin
      if c = 0.0 and then d = 0.0 then
         raise Constraint_Error;
      else
         return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
                         Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
      end if;
   end "/";
 
   function "/" (Left, Right : Imaginary) return Real'Base is
   begin
      return R (Left) / R (Right);
   end "/";
 
   function "/" (Left : Complex; Right : Real'Base) return Complex is
   begin
      return Complex'(Left.Re / Right, Left.Im / Right);
   end "/";
 
   function "/" (Left : Real'Base; Right : Complex) return Complex is
      a : constant R := Left;
      c : constant R := Right.Re;
      d : constant R := Right.Im;
   begin
      return Complex'(Re =>   (a * c) / (c ** 2 + d ** 2),
                      Im => -((a * d) / (c ** 2 + d ** 2)));
   end "/";
 
   function "/" (Left : Complex; Right : Imaginary) return Complex is
      a : constant R := Left.Re;
      b : constant R := Left.Im;
      d : constant R := R (Right);
 
   begin
      return (b / d,  -(a / d));
   end "/";
 
   function "/" (Left : Imaginary; Right : Complex) return Complex is
      b : constant R := R (Left);
      c : constant R := Right.Re;
      d : constant R := Right.Im;
 
   begin
      return (Re => b * d / (c ** 2 + d ** 2),
              Im => b * c / (c ** 2 + d ** 2));
   end "/";
 
   function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
   begin
      return Imaginary (R (Left) / Right);
   end "/";
 
   function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
   begin
      return Imaginary (-(Left / R (Right)));
   end "/";
 
   ---------
   -- "<" --
   ---------
 
   function "<" (Left, Right : Imaginary) return Boolean is
   begin
      return R (Left) < R (Right);
   end "<";
 
   ----------
   -- "<=" --
   ----------
 
   function "<=" (Left, Right : Imaginary) return Boolean is
   begin
      return R (Left) <= R (Right);
   end "<=";
 
   ---------
   -- ">" --
   ---------
 
   function ">" (Left, Right : Imaginary) return Boolean is
   begin
      return R (Left) > R (Right);
   end ">";
 
   ----------
   -- ">=" --
   ----------
 
   function ">=" (Left, Right : Imaginary) return Boolean is
   begin
      return R (Left) >= R (Right);
   end ">=";
 
   -----------
   -- "abs" --
   -----------
 
   function "abs" (Right : Imaginary) return Real'Base is
   begin
      return abs R (Right);
   end "abs";
 
   --------------
   -- Argument --
   --------------
 
   function Argument (X : Complex) return Real'Base is
      a   : constant R := X.Re;
      b   : constant R := X.Im;
      arg : R;
 
   begin
      if b = 0.0 then
 
         if a >= 0.0 then
            return 0.0;
         else
            return R'Copy_Sign (Pi, b);
         end if;
 
      elsif a = 0.0 then
 
         if b >= 0.0 then
            return Half_Pi;
         else
            return -Half_Pi;
         end if;
 
      else
         arg := R (Atan (Double (abs (b / a))));
 
         if a > 0.0 then
            if b > 0.0 then
               return arg;
            else                  --  b < 0.0
               return -arg;
            end if;
 
         else                     --  a < 0.0
            if b >= 0.0 then
               return Pi - arg;
            else                  --  b < 0.0
               return -(Pi - arg);
            end if;
         end if;
      end if;
 
   exception
      when Constraint_Error =>
         if b > 0.0 then
            return Half_Pi;
         else
            return -Half_Pi;
         end if;
   end Argument;
 
   function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
   begin
      if Cycle > 0.0 then
         return Argument (X) * Cycle / Two_Pi;
      else
         raise Argument_Error;
      end if;
   end Argument;
 
   ----------------------------
   -- Compose_From_Cartesian --
   ----------------------------
 
   function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
   begin
      return (Re, Im);
   end Compose_From_Cartesian;
 
   function Compose_From_Cartesian (Re : Real'Base) return Complex is
   begin
      return (Re, 0.0);
   end Compose_From_Cartesian;
 
   function Compose_From_Cartesian (Im : Imaginary) return Complex is
   begin
      return (0.0, R (Im));
   end Compose_From_Cartesian;
 
   ------------------------
   -- Compose_From_Polar --
   ------------------------
 
   function Compose_From_Polar (
     Modulus, Argument : Real'Base)
     return Complex
   is
   begin
      if Modulus = 0.0 then
         return (0.0, 0.0);
      else
         return (Modulus * R (Cos (Double (Argument))),
                 Modulus * R (Sin (Double (Argument))));
      end if;
   end Compose_From_Polar;
 
   function Compose_From_Polar (
     Modulus, Argument, Cycle : Real'Base)
     return Complex
   is
      Arg : Real'Base;
 
   begin
      if Modulus = 0.0 then
         return (0.0, 0.0);
 
      elsif Cycle > 0.0 then
         if Argument = 0.0 then
            return (Modulus, 0.0);
 
         elsif Argument = Cycle / 4.0 then
            return (0.0, Modulus);
 
         elsif Argument = Cycle / 2.0 then
            return (-Modulus, 0.0);
 
         elsif Argument = 3.0 * Cycle / R (4.0) then
            return (0.0, -Modulus);
         else
            Arg := Two_Pi * Argument / Cycle;
            return (Modulus * R (Cos (Double (Arg))),
                    Modulus * R (Sin (Double (Arg))));
         end if;
      else
         raise Argument_Error;
      end if;
   end Compose_From_Polar;
 
   ---------------
   -- Conjugate --
   ---------------
 
   function Conjugate (X : Complex) return Complex is
   begin
      return Complex'(X.Re, -X.Im);
   end Conjugate;
 
   --------
   -- Im --
   --------
 
   function Im (X : Complex) return Real'Base is
   begin
      return X.Im;
   end Im;
 
   function Im (X : Imaginary) return Real'Base is
   begin
      return R (X);
   end Im;
 
   -------------
   -- Modulus --
   -------------
 
   function Modulus (X : Complex) return Real'Base is
      Re2, Im2 : R;
 
   begin
 
      begin
         Re2 := X.Re ** 2;
 
         --  To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
         --  compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
         --  squaring does not raise constraint_error but generates infinity,
         --  we can use an explicit comparison to determine whether to use
         --  the scaling expression.
 
         --  The scaling expression is computed in double format throughout
         --  in order to prevent inaccuracies on machines where not all
         --  immediate expressions are rounded, such as PowerPC.
 
         if Re2 > R'Last then
            raise Constraint_Error;
         end if;
 
      exception
         when Constraint_Error =>
            return R (Double (abs (X.Re))
              * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
      end;
 
      begin
         Im2 := X.Im ** 2;
 
         if Im2 > R'Last then
            raise Constraint_Error;
         end if;
 
      exception
         when Constraint_Error =>
            return R (Double (abs (X.Im))
              * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
      end;
 
      --  Now deal with cases of underflow. If only one of the squares
      --  underflows, return the modulus of the other component. If both
      --  squares underflow, use scaling as above.
 
      if Re2 = 0.0 then
 
         if X.Re = 0.0 then
            return abs (X.Im);
 
         elsif Im2 = 0.0 then
 
            if X.Im = 0.0 then
               return abs (X.Re);
 
            else
               if abs (X.Re) > abs (X.Im) then
                  return
                    R (Double (abs (X.Re))
                      * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
               else
                  return
                    R (Double (abs (X.Im))
                      * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
               end if;
            end if;
 
         else
            return abs (X.Im);
         end if;
 
      elsif Im2 = 0.0 then
         return abs (X.Re);
 
      --  In all other cases, the naive computation will do
 
      else
         return R (Sqrt (Double (Re2 + Im2)));
      end if;
   end Modulus;
 
   --------
   -- Re --
   --------
 
   function Re (X : Complex) return Real'Base is
   begin
      return X.Re;
   end Re;
 
   ------------
   -- Set_Im --
   ------------
 
   procedure Set_Im (X : in out Complex; Im : Real'Base) is
   begin
      X.Im := Im;
   end Set_Im;
 
   procedure Set_Im (X : out Imaginary; Im : Real'Base) is
   begin
      X := Imaginary (Im);
   end Set_Im;
 
   ------------
   -- Set_Re --
   ------------
 
   procedure Set_Re (X : in out Complex; Re : Real'Base) is
   begin
      X.Re := Re;
   end Set_Re;
 
end Ada.Numerics.Generic_Complex_Types;

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [gcc/] [ada/] [a-ngcoty.adb] - Blame information for rev 826

Line No.	Rev	Author	Line
1	281	jeremybenn	`------------------------------------------------------------------------------`
2			`-- --`
3			`-- GNAT RUN-TIME COMPONENTS --`
4			`-- --`
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 --`
6			`-- --`
7			`-- B o d y --`
8			`-- --`
9			`-- Copyright (C) 1992-2009, Free Software Foundation, Inc. --`
10			`-- --`
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. --`
17			`-- --`
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. --`
21			`-- --`
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/>. --`
26			`-- --`
27			`-- GNAT was originally developed by the GNAT team at New York University. --`
28			`-- Extensive contributions were provided by Ada Core Technologies Inc. --`
29			`-- --`
30			`------------------------------------------------------------------------------`
31
32			`with Ada.Numerics.Aux; use Ada.Numerics.Aux;`
33
34			`package body Ada.Numerics.Generic_Complex_Types is`
35
36			`subtype R is Real'Base;`
37
38			`Two_Pi : constant R := R (2.0) * Pi;`
39			`Half_Pi : constant R := Pi / R (2.0);`
40
41			`---------`
42			`-- "*" --`
43			`---------`
44
45			`function "*" (Left, Right : Complex) return Complex is`
46			`X : R;`
47			`Y : R;`
48
49			`begin`
50			`X := Left.Re * Right.Re - Left.Im * Right.Im;`
51			`Y := Left.Re * Right.Im + Left.Im * Right.Re;`
52
53			`-- If either component overflows, try to scale (skip in fast math mode)`
54
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));`
59			`end if;`
60
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));`
64			`end if;`
65			`end if;`
66
67			`return (X, Y);`
68			`end "*";`
69
70			`function "*" (Left, Right : Imaginary) return Real'Base is`
71			`begin`
72			`return -(R (Left) * R (Right));`
73			`end "*";`
74
75			`function "*" (Left : Complex; Right : Real'Base) return Complex is`
76			`begin`
77			`return Complex'(Left.Re * Right, Left.Im * Right);`
78			`end "*";`
79
80			`function "*" (Left : Real'Base; Right : Complex) return Complex is`
81			`begin`
82			`return (Left * Right.Re, Left * Right.Im);`
83			`end "*";`
84
85			`function "*" (Left : Complex; Right : Imaginary) return Complex is`
86			`begin`
87			`return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));`
88			`end "*";`
89
90			`function "*" (Left : Imaginary; Right : Complex) return Complex is`
91			`begin`
92			`return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);`
93			`end "*";`
94
95			`function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is`
96			`begin`
97			`return Left * Imaginary (Right);`
98			`end "*";`
99
100			`function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is`
101			`begin`
102			`return Imaginary (Left * R (Right));`
103			`end "*";`
104
105			`----------`
106			`-- "**" --`
107			`----------`
108
109			`function "**" (Left : Complex; Right : Integer) return Complex is`
110			`Result : Complex := (1.0, 0.0);`
111			`Factor : Complex := Left;`
112			`Exp : Integer := Right;`
113
114			`begin`
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.`
120
121			`if Exp >= 0 then`
122
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.`
126
127			`while Exp /= 0 loop`
128			`if Exp rem 2 /= 0 then`
129			`Result := Result * Factor;`
130			`end if;`
131
132			`Factor := Factor * Factor;`
133			`Exp := Exp / 2;`
134			`end loop;`
135
136			`return Result;`
137
138			`else -- Exp < 0 then`
139
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.`
144
145			`begin`
146			`while Exp /= 0 loop`
147			`if Exp rem 2 /= 0 then`
148			`Result := Result * Factor;`
149			`end if;`
150
151			`Factor := Factor * Factor;`
152			`Exp := Exp / 2;`
153			`end loop;`
154
155			`return R'(1.0) / Result;`
156
157			`exception`
158			`when Constraint_Error =>`
159			`return (0.0, 0.0);`
160			`end;`
161			`end if;`
162			`end "**";`
163
164			`function "**" (Left : Imaginary; Right : Integer) return Complex is`
165			`M : constant R := R (Left) ** Right;`
166			`begin`
167			`case Right mod 4 is`
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;`
173			`end case;`
174			`end "**";`
175
176			`---------`
177			`-- "+" --`
178			`---------`
179
180			`function "+" (Right : Complex) return Complex is`
181			`begin`
182			`return Right;`
183			`end "+";`
184
185			`function "+" (Left, Right : Complex) return Complex is`
186			`begin`
187			`return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);`
188			`end "+";`
189
190			`function "+" (Right : Imaginary) return Imaginary is`
191			`begin`
192			`return Right;`
193			`end "+";`
194
195			`function "+" (Left, Right : Imaginary) return Imaginary is`
196			`begin`
197			`return Imaginary (R (Left) + R (Right));`
198			`end "+";`
199
200			`function "+" (Left : Complex; Right : Real'Base) return Complex is`
201			`begin`
202			`return Complex'(Left.Re + Right, Left.Im);`
203			`end "+";`
204
205			`function "+" (Left : Real'Base; Right : Complex) return Complex is`
206			`begin`
207			`return Complex'(Left + Right.Re, Right.Im);`
208			`end "+";`
209
210			`function "+" (Left : Complex; Right : Imaginary) return Complex is`
211			`begin`
212			`return Complex'(Left.Re, Left.Im + R (Right));`
213			`end "+";`
214
215			`function "+" (Left : Imaginary; Right : Complex) return Complex is`
216			`begin`
217			`return Complex'(Right.Re, R (Left) + Right.Im);`
218			`end "+";`
219
220			`function "+" (Left : Imaginary; Right : Real'Base) return Complex is`
221			`begin`
222			`return Complex'(Right, R (Left));`
223			`end "+";`
224
225			`function "+" (Left : Real'Base; Right : Imaginary) return Complex is`
226			`begin`
227			`return Complex'(Left, R (Right));`
228			`end "+";`
229
230			`---------`
231			`-- "-" --`
232			`---------`
233
234			`function "-" (Right : Complex) return Complex is`
235			`begin`
236			`return (-Right.Re, -Right.Im);`
237			`end "-";`
238
239			`function "-" (Left, Right : Complex) return Complex is`
240			`begin`
241			`return (Left.Re - Right.Re, Left.Im - Right.Im);`
242			`end "-";`
243
244			`function "-" (Right : Imaginary) return Imaginary is`
245			`begin`
246			`return Imaginary (-R (Right));`
247			`end "-";`
248
249			`function "-" (Left, Right : Imaginary) return Imaginary is`
250			`begin`
251			`return Imaginary (R (Left) - R (Right));`
252			`end "-";`
253
254			`function "-" (Left : Complex; Right : Real'Base) return Complex is`
255			`begin`
256			`return Complex'(Left.Re - Right, Left.Im);`
257			`end "-";`
258
259			`function "-" (Left : Real'Base; Right : Complex) return Complex is`
260			`begin`
261			`return Complex'(Left - Right.Re, -Right.Im);`
262			`end "-";`
263
264			`function "-" (Left : Complex; Right : Imaginary) return Complex is`
265			`begin`
266			`return Complex'(Left.Re, Left.Im - R (Right));`
267			`end "-";`
268
269			`function "-" (Left : Imaginary; Right : Complex) return Complex is`
270			`begin`
271			`return Complex'(-Right.Re, R (Left) - Right.Im);`
272			`end "-";`
273
274			`function "-" (Left : Imaginary; Right : Real'Base) return Complex is`
275			`begin`
276			`return Complex'(-Right, R (Left));`
277			`end "-";`
278
279			`function "-" (Left : Real'Base; Right : Imaginary) return Complex is`
280			`begin`
281			`return Complex'(Left, -R (Right));`
282			`end "-";`
283
284			`---------`
285			`-- "/" --`
286			`---------`
287
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;`
293
294			`begin`
295			`if c = 0.0 and then d = 0.0 then`
296			`raise Constraint_Error;`
297			`else`
298			`return Complex'(Re => ((a * c) + (b * d)) / (c 2 + d 2),`
299			`Im => ((b * c) - (a * d)) / (c 2 + d 2));`
300			`end if;`
301			`end "/";`
302
303			`function "/" (Left, Right : Imaginary) return Real'Base is`
304			`begin`
305			`return R (Left) / R (Right);`
306			`end "/";`
307
308			`function "/" (Left : Complex; Right : Real'Base) return Complex is`
309			`begin`
310			`return Complex'(Left.Re / Right, Left.Im / Right);`
311			`end "/";`
312
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;`
317			`begin`
318			`return Complex'(Re => (a * c) / (c 2 + d 2),`
319			`Im => -((a * d) / (c 2 + d 2)));`
320			`end "/";`
321
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);`
326
327			`begin`
328			`return (b / d, -(a / d));`
329			`end "/";`
330
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;`
335
336			`begin`
337			`return (Re => b * d / (c 2 + d 2),`
338			`Im => b * c / (c 2 + d 2));`
339			`end "/";`
340
341			`function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is`
342			`begin`
343			`return Imaginary (R (Left) / Right);`
344			`end "/";`
345
346			`function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is`
347			`begin`
348			`return Imaginary (-(Left / R (Right)));`
349			`end "/";`
350
351			`---------`
352			`-- "<" --`
353			`---------`
354
355			`function "<" (Left, Right : Imaginary) return Boolean is`
356			`begin`
357			`return R (Left) < R (Right);`
358			`end "<";`
359
360			`----------`
361			`-- "<=" --`
362			`----------`
363
364			`function "<=" (Left, Right : Imaginary) return Boolean is`
365			`begin`
366			`return R (Left) <= R (Right);`
367			`end "<=";`
368
369			`---------`
370			`-- ">" --`
371			`---------`
372
373			`function ">" (Left, Right : Imaginary) return Boolean is`
374			`begin`
375			`return R (Left) > R (Right);`
376			`end ">";`
377
378			`----------`
379			`-- ">=" --`
380			`----------`
381
382			`function ">=" (Left, Right : Imaginary) return Boolean is`
383			`begin`
384			`return R (Left) >= R (Right);`
385			`end ">=";`
386
387			`-----------`
388			`-- "abs" --`
389			`-----------`
390
391			`function "abs" (Right : Imaginary) return Real'Base is`
392			`begin`
393			`return abs R (Right);`
394			`end "abs";`
395
396			`--------------`
397			`-- Argument --`
398			`--------------`
399
400			`function Argument (X : Complex) return Real'Base is`
401			`a : constant R := X.Re;`
402			`b : constant R := X.Im;`
403			`arg : R;`
404
405			`begin`
406			`if b = 0.0 then`
407
408			`if a >= 0.0 then`
409			`return 0.0;`
410			`else`
411			`return R'Copy_Sign (Pi, b);`
412			`end if;`
413
414			`elsif a = 0.0 then`
415
416			`if b >= 0.0 then`
417			`return Half_Pi;`
418			`else`
419			`return -Half_Pi;`
420			`end if;`
421
422			`else`
423			`arg := R (Atan (Double (abs (b / a))));`
424
425			`if a > 0.0 then`
426			`if b > 0.0 then`
427			`return arg;`
428			`else -- b < 0.0`
429			`return -arg;`
430			`end if;`
431
432			`else -- a < 0.0`
433			`if b >= 0.0 then`
434			`return Pi - arg;`
435			`else -- b < 0.0`
436			`return -(Pi - arg);`
437			`end if;`
438			`end if;`
439			`end if;`
440
441			`exception`
442			`when Constraint_Error =>`
443			`if b > 0.0 then`
444			`return Half_Pi;`
445			`else`
446			`return -Half_Pi;`
447			`end if;`
448			`end Argument;`
449
450			`function Argument (X : Complex; Cycle : Real'Base) return Real'Base is`
451			`begin`
452			`if Cycle > 0.0 then`
453			`return Argument (X) * Cycle / Two_Pi;`
454			`else`
455			`raise Argument_Error;`
456			`end if;`
457			`end Argument;`
458
459			`----------------------------`
460			`-- Compose_From_Cartesian --`
461			`----------------------------`
462
463			`function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is`
464			`begin`
465			`return (Re, Im);`
466			`end Compose_From_Cartesian;`
467
468			`function Compose_From_Cartesian (Re : Real'Base) return Complex is`
469			`begin`
470			`return (Re, 0.0);`
471			`end Compose_From_Cartesian;`
472
473			`function Compose_From_Cartesian (Im : Imaginary) return Complex is`
474			`begin`
475			`return (0.0, R (Im));`
476			`end Compose_From_Cartesian;`
477
478			`------------------------`
479			`-- Compose_From_Polar --`
480			`------------------------`
481
482			`function Compose_From_Polar (`
483			`Modulus, Argument : Real'Base)`
484			`return Complex`
485			`is`
486			`begin`
487			`if Modulus = 0.0 then`
488			`return (0.0, 0.0);`
489			`else`
490			`return (Modulus * R (Cos (Double (Argument))),`
491			`Modulus * R (Sin (Double (Argument))));`
492			`end if;`
493			`end Compose_From_Polar;`
494
495			`function Compose_From_Polar (`
496			`Modulus, Argument, Cycle : Real'Base)`
497			`return Complex`
498			`is`
499			`Arg : Real'Base;`
500
501			`begin`
502			`if Modulus = 0.0 then`
503			`return (0.0, 0.0);`
504
505			`elsif Cycle > 0.0 then`
506			`if Argument = 0.0 then`
507			`return (Modulus, 0.0);`
508
509			`elsif Argument = Cycle / 4.0 then`
510			`return (0.0, Modulus);`
511
512			`elsif Argument = Cycle / 2.0 then`
513			`return (-Modulus, 0.0);`
514
515			`elsif Argument = 3.0 * Cycle / R (4.0) then`
516			`return (0.0, -Modulus);`
517			`else`
518			`Arg := Two_Pi * Argument / Cycle;`
519			`return (Modulus * R (Cos (Double (Arg))),`
520			`Modulus * R (Sin (Double (Arg))));`
521			`end if;`
522			`else`
523			`raise Argument_Error;`
524			`end if;`
525			`end Compose_From_Polar;`
526
527			`---------------`
528			`-- Conjugate --`
529			`---------------`
530
531			`function Conjugate (X : Complex) return Complex is`
532			`begin`
533			`return Complex'(X.Re, -X.Im);`
534			`end Conjugate;`
535
536			`--------`
537			`-- Im --`
538			`--------`
539
540			`function Im (X : Complex) return Real'Base is`
541			`begin`
542			`return X.Im;`
543			`end Im;`
544
545			`function Im (X : Imaginary) return Real'Base is`
546			`begin`
547			`return R (X);`
548			`end Im;`
549
550			`-------------`
551			`-- Modulus --`
552			`-------------`
553
554			`function Modulus (X : Complex) return Real'Base is`
555			`Re2, Im2 : R;`
556
557			`begin`
558
559			`begin`
560			`Re2 := X.Re ** 2;`
561
562			`-- To compute (a2 + b2) (0.5) when a2 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.`
567
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.`
571
572			`if Re2 > R'Last then`
573			`raise Constraint_Error;`
574			`end if;`
575
576			`exception`
577			`when Constraint_Error =>`
578			`return R (Double (abs (X.Re))`
579			`* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));`
580			`end;`
581
582			`begin`
583			`Im2 := X.Im ** 2;`
584
585			`if Im2 > R'Last then`
586			`raise Constraint_Error;`
587			`end if;`
588
589			`exception`
590			`when Constraint_Error =>`
591			`return R (Double (abs (X.Im))`
592			`* Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));`
593			`end;`
594
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.`
598
599			`if Re2 = 0.0 then`
600
601			`if X.Re = 0.0 then`
602			`return abs (X.Im);`
603
604			`elsif Im2 = 0.0 then`
605
606			`if X.Im = 0.0 then`
607			`return abs (X.Re);`
608
609			`else`
610			`if abs (X.Re) > abs (X.Im) then`
611			`return`
612			`R (Double (abs (X.Re))`
613			`* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));`
614			`else`
615			`return`
616			`R (Double (abs (X.Im))`
617			`* Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));`
618			`end if;`
619			`end if;`
620
621			`else`
622			`return abs (X.Im);`
623			`end if;`
624
625			`elsif Im2 = 0.0 then`
626			`return abs (X.Re);`
627
628			`-- In all other cases, the naive computation will do`
629
630			`else`
631			`return R (Sqrt (Double (Re2 + Im2)));`
632			`end if;`
633			`end Modulus;`
634
635			`--------`
636			`-- Re --`
637			`--------`
638
639			`function Re (X : Complex) return Real'Base is`
640			`begin`
641			`return X.Re;`
642			`end Re;`
643
644			`------------`
645			`-- Set_Im --`
646			`------------`
647
648			`procedure Set_Im (X : in out Complex; Im : Real'Base) is`
649			`begin`
650			`X.Im := Im;`
651			`end Set_Im;`
652
653			`procedure Set_Im (X : out Imaginary; Im : Real'Base) is`
654			`begin`
655			`X := Imaginary (Im);`
656			`end Set_Im;`
657
658			`------------`
659			`-- Set_Re --`
660			`------------`
661
662			`procedure Set_Re (X : in out Complex; Re : Real'Base) is`
663			`begin`
664			`X.Re := Re;`
665			`end Set_Re;`
666
667			`end Ada.Numerics.Generic_Complex_Types;`