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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [gcc/] [ada/] [a-ngcoty.adb] - Blame information for rev 706

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-2010, 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
 
      Scale : constant R := R (R'Machine_Radix) ** ((R'Machine_Emax - 1) / 2);
      --  In case of overflow, scale the operands by the largest power of the
      --  radix (to avoid rounding error), so that the square of the scale does
      --  not overflow itself.
 
      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
 
         --  Note that the test below is written as a negation. This is to
         --  account for the fact that X and Y may be NaNs, because both of
         --  their operands could overflow. Given that all operations on NaNs
         --  return false, the test can only be written thus.
 
         if not (abs (X) <= R'Last) then
            X := Scale**2 * ((Left.Re / Scale) * (Right.Re / Scale) -
                             (Left.Im / Scale) * (Right.Im / Scale));
         end if;
 
         if not (abs (Y) <= R'Last) then
            Y := Scale**2 * ((Left.Re / Scale) * (Right.Im / Scale)
                           + (Left.Im / Scale) * (Right.Re / Scale));
         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.
 
         --  ??? same weird test, why not Re2 > R'Last ???
         if not (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;
 
         --  ??? same weird test
         if not (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-dev/] [or1k-gcc/] [gcc/] [ada/] [a-ngcoty.adb] - Blame information for rev 706

Line No.	Rev	Author	Line
1	706	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-2010, 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
47			`Scale : constant R := R (R'Machine_Radix) ** ((R'Machine_Emax - 1) / 2);`
48			`-- In case of overflow, scale the operands by the largest power of the`
49			`-- radix (to avoid rounding error), so that the square of the scale does`
50			`-- not overflow itself.`
51
52			`X : R;`
53			`Y : R;`
54
55			`begin`
56			`X := Left.Re * Right.Re - Left.Im * Right.Im;`
57			`Y := Left.Re * Right.Im + Left.Im * Right.Re;`
58
59			`-- If either component overflows, try to scale (skip in fast math mode)`
60
61			`if not Standard'Fast_Math then`
62
63			`-- Note that the test below is written as a negation. This is to`
64			`-- account for the fact that X and Y may be NaNs, because both of`
65			`-- their operands could overflow. Given that all operations on NaNs`
66			`-- return false, the test can only be written thus.`
67
68			`if not (abs (X) <= R'Last) then`
69			`X := Scale*2 ((Left.Re / Scale) * (Right.Re / Scale) -`
70			`(Left.Im / Scale) * (Right.Im / Scale));`
71			`end if;`
72
73			`if not (abs (Y) <= R'Last) then`
74			`Y := Scale*2 ((Left.Re / Scale) * (Right.Im / Scale)`
75			`+ (Left.Im / Scale) * (Right.Re / Scale));`
76			`end if;`
77			`end if;`
78
79			`return (X, Y);`
80			`end "*";`
81
82			`function "*" (Left, Right : Imaginary) return Real'Base is`
83			`begin`
84			`return -(R (Left) * R (Right));`
85			`end "*";`
86
87			`function "*" (Left : Complex; Right : Real'Base) return Complex is`
88			`begin`
89			`return Complex'(Left.Re * Right, Left.Im * Right);`
90			`end "*";`
91
92			`function "*" (Left : Real'Base; Right : Complex) return Complex is`
93			`begin`
94			`return (Left * Right.Re, Left * Right.Im);`
95			`end "*";`
96
97			`function "*" (Left : Complex; Right : Imaginary) return Complex is`
98			`begin`
99			`return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));`
100			`end "*";`
101
102			`function "*" (Left : Imaginary; Right : Complex) return Complex is`
103			`begin`
104			`return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);`
105			`end "*";`
106
107			`function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is`
108			`begin`
109			`return Left * Imaginary (Right);`
110			`end "*";`
111
112			`function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is`
113			`begin`
114			`return Imaginary (Left * R (Right));`
115			`end "*";`
116
117			`----------`
118			`-- "**" --`
119			`----------`
120
121			`function "**" (Left : Complex; Right : Integer) return Complex is`
122			`Result : Complex := (1.0, 0.0);`
123			`Factor : Complex := Left;`
124			`Exp : Integer := Right;`
125
126			`begin`
127			`-- We use the standard logarithmic approach, Exp gets shifted right`
128			`-- testing successive low order bits and Factor is the value of the`
129			`-- base raised to the next power of 2. For positive exponents we`
130			`-- multiply the result by this factor, for negative exponents, we`
131			`-- divide by this factor.`
132
133			`if Exp >= 0 then`
134
135			`-- For a positive exponent, if we get a constraint error during`
136			`-- this loop, it is an overflow, and the constraint error will`
137			`-- simply be passed on to the caller.`
138
139			`while Exp /= 0 loop`
140			`if Exp rem 2 /= 0 then`
141			`Result := Result * Factor;`
142			`end if;`
143
144			`Factor := Factor * Factor;`
145			`Exp := Exp / 2;`
146			`end loop;`
147
148			`return Result;`
149
150			`else -- Exp < 0 then`
151
152			`-- For the negative exponent case, a constraint error during this`
153			`-- calculation happens if Factor gets too large, and the proper`
154			`-- response is to return 0.0, since what we essentially have is`
155			`-- 1.0 / infinity, and the closest model number will be zero.`
156
157			`begin`
158			`while Exp /= 0 loop`
159			`if Exp rem 2 /= 0 then`
160			`Result := Result * Factor;`
161			`end if;`
162
163			`Factor := Factor * Factor;`
164			`Exp := Exp / 2;`
165			`end loop;`
166
167			`return R'(1.0) / Result;`
168
169			`exception`
170			`when Constraint_Error =>`
171			`return (0.0, 0.0);`
172			`end;`
173			`end if;`
174			`end "**";`
175
176			`function "**" (Left : Imaginary; Right : Integer) return Complex is`
177			`M : constant R := R (Left) ** Right;`
178			`begin`
179			`case Right mod 4 is`
180			`when 0 => return (M, 0.0);`
181			`when 1 => return (0.0, M);`
182			`when 2 => return (-M, 0.0);`
183			`when 3 => return (0.0, -M);`
184			`when others => raise Program_Error;`
185			`end case;`
186			`end "**";`
187
188			`---------`
189			`-- "+" --`
190			`---------`
191
192			`function "+" (Right : Complex) return Complex is`
193			`begin`
194			`return Right;`
195			`end "+";`
196
197			`function "+" (Left, Right : Complex) return Complex is`
198			`begin`
199			`return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);`
200			`end "+";`
201
202			`function "+" (Right : Imaginary) return Imaginary is`
203			`begin`
204			`return Right;`
205			`end "+";`
206
207			`function "+" (Left, Right : Imaginary) return Imaginary is`
208			`begin`
209			`return Imaginary (R (Left) + R (Right));`
210			`end "+";`
211
212			`function "+" (Left : Complex; Right : Real'Base) return Complex is`
213			`begin`
214			`return Complex'(Left.Re + Right, Left.Im);`
215			`end "+";`
216
217			`function "+" (Left : Real'Base; Right : Complex) return Complex is`
218			`begin`
219			`return Complex'(Left + Right.Re, Right.Im);`
220			`end "+";`
221
222			`function "+" (Left : Complex; Right : Imaginary) return Complex is`
223			`begin`
224			`return Complex'(Left.Re, Left.Im + R (Right));`
225			`end "+";`
226
227			`function "+" (Left : Imaginary; Right : Complex) return Complex is`
228			`begin`
229			`return Complex'(Right.Re, R (Left) + Right.Im);`
230			`end "+";`
231
232			`function "+" (Left : Imaginary; Right : Real'Base) return Complex is`
233			`begin`
234			`return Complex'(Right, R (Left));`
235			`end "+";`
236
237			`function "+" (Left : Real'Base; Right : Imaginary) return Complex is`
238			`begin`
239			`return Complex'(Left, R (Right));`
240			`end "+";`
241
242			`---------`
243			`-- "-" --`
244			`---------`
245
246			`function "-" (Right : Complex) return Complex is`
247			`begin`
248			`return (-Right.Re, -Right.Im);`
249			`end "-";`
250
251			`function "-" (Left, Right : Complex) return Complex is`
252			`begin`
253			`return (Left.Re - Right.Re, Left.Im - Right.Im);`
254			`end "-";`
255
256			`function "-" (Right : Imaginary) return Imaginary is`
257			`begin`
258			`return Imaginary (-R (Right));`
259			`end "-";`
260
261			`function "-" (Left, Right : Imaginary) return Imaginary is`
262			`begin`
263			`return Imaginary (R (Left) - R (Right));`
264			`end "-";`
265
266			`function "-" (Left : Complex; Right : Real'Base) return Complex is`
267			`begin`
268			`return Complex'(Left.Re - Right, Left.Im);`
269			`end "-";`
270
271			`function "-" (Left : Real'Base; Right : Complex) return Complex is`
272			`begin`
273			`return Complex'(Left - Right.Re, -Right.Im);`
274			`end "-";`
275
276			`function "-" (Left : Complex; Right : Imaginary) return Complex is`
277			`begin`
278			`return Complex'(Left.Re, Left.Im - R (Right));`
279			`end "-";`
280
281			`function "-" (Left : Imaginary; Right : Complex) return Complex is`
282			`begin`
283			`return Complex'(-Right.Re, R (Left) - Right.Im);`
284			`end "-";`
285
286			`function "-" (Left : Imaginary; Right : Real'Base) return Complex is`
287			`begin`
288			`return Complex'(-Right, R (Left));`
289			`end "-";`
290
291			`function "-" (Left : Real'Base; Right : Imaginary) return Complex is`
292			`begin`
293			`return Complex'(Left, -R (Right));`
294			`end "-";`
295
296			`---------`
297			`-- "/" --`
298			`---------`
299
300			`function "/" (Left, Right : Complex) return Complex is`
301			`a : constant R := Left.Re;`
302			`b : constant R := Left.Im;`
303			`c : constant R := Right.Re;`
304			`d : constant R := Right.Im;`
305
306			`begin`
307			`if c = 0.0 and then d = 0.0 then`
308			`raise Constraint_Error;`
309			`else`
310			`return Complex'(Re => ((a * c) + (b * d)) / (c 2 + d 2),`
311			`Im => ((b * c) - (a * d)) / (c 2 + d 2));`
312			`end if;`
313			`end "/";`
314
315			`function "/" (Left, Right : Imaginary) return Real'Base is`
316			`begin`
317			`return R (Left) / R (Right);`
318			`end "/";`
319
320			`function "/" (Left : Complex; Right : Real'Base) return Complex is`
321			`begin`
322			`return Complex'(Left.Re / Right, Left.Im / Right);`
323			`end "/";`
324
325			`function "/" (Left : Real'Base; Right : Complex) return Complex is`
326			`a : constant R := Left;`
327			`c : constant R := Right.Re;`
328			`d : constant R := Right.Im;`
329			`begin`
330			`return Complex'(Re => (a * c) / (c 2 + d 2),`
331			`Im => -((a * d) / (c 2 + d 2)));`
332			`end "/";`
333
334			`function "/" (Left : Complex; Right : Imaginary) return Complex is`
335			`a : constant R := Left.Re;`
336			`b : constant R := Left.Im;`
337			`d : constant R := R (Right);`
338
339			`begin`
340			`return (b / d, -(a / d));`
341			`end "/";`
342
343			`function "/" (Left : Imaginary; Right : Complex) return Complex is`
344			`b : constant R := R (Left);`
345			`c : constant R := Right.Re;`
346			`d : constant R := Right.Im;`
347
348			`begin`
349			`return (Re => b * d / (c 2 + d 2),`
350			`Im => b * c / (c 2 + d 2));`
351			`end "/";`
352
353			`function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is`
354			`begin`
355			`return Imaginary (R (Left) / Right);`
356			`end "/";`
357
358			`function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is`
359			`begin`
360			`return Imaginary (-(Left / R (Right)));`
361			`end "/";`
362
363			`---------`
364			`-- "<" --`
365			`---------`
366
367			`function "<" (Left, Right : Imaginary) return Boolean is`
368			`begin`
369			`return R (Left) < R (Right);`
370			`end "<";`
371
372			`----------`
373			`-- "<=" --`
374			`----------`
375
376			`function "<=" (Left, Right : Imaginary) return Boolean is`
377			`begin`
378			`return R (Left) <= R (Right);`
379			`end "<=";`
380
381			`---------`
382			`-- ">" --`
383			`---------`
384
385			`function ">" (Left, Right : Imaginary) return Boolean is`
386			`begin`
387			`return R (Left) > R (Right);`
388			`end ">";`
389
390			`----------`
391			`-- ">=" --`
392			`----------`
393
394			`function ">=" (Left, Right : Imaginary) return Boolean is`
395			`begin`
396			`return R (Left) >= R (Right);`
397			`end ">=";`
398
399			`-----------`
400			`-- "abs" --`
401			`-----------`
402
403			`function "abs" (Right : Imaginary) return Real'Base is`
404			`begin`
405			`return abs R (Right);`
406			`end "abs";`
407
408			`--------------`
409			`-- Argument --`
410			`--------------`
411
412			`function Argument (X : Complex) return Real'Base is`
413			`a : constant R := X.Re;`
414			`b : constant R := X.Im;`
415			`arg : R;`
416
417			`begin`
418			`if b = 0.0 then`
419
420			`if a >= 0.0 then`
421			`return 0.0;`
422			`else`
423			`return R'Copy_Sign (Pi, b);`
424			`end if;`
425
426			`elsif a = 0.0 then`
427
428			`if b >= 0.0 then`
429			`return Half_Pi;`
430			`else`
431			`return -Half_Pi;`
432			`end if;`
433
434			`else`
435			`arg := R (Atan (Double (abs (b / a))));`
436
437			`if a > 0.0 then`
438			`if b > 0.0 then`
439			`return arg;`
440			`else -- b < 0.0`
441			`return -arg;`
442			`end if;`
443
444			`else -- a < 0.0`
445			`if b >= 0.0 then`
446			`return Pi - arg;`
447			`else -- b < 0.0`
448			`return -(Pi - arg);`
449			`end if;`
450			`end if;`
451			`end if;`
452
453			`exception`
454			`when Constraint_Error =>`
455			`if b > 0.0 then`
456			`return Half_Pi;`
457			`else`
458			`return -Half_Pi;`
459			`end if;`
460			`end Argument;`
461
462			`function Argument (X : Complex; Cycle : Real'Base) return Real'Base is`
463			`begin`
464			`if Cycle > 0.0 then`
465			`return Argument (X) * Cycle / Two_Pi;`
466			`else`
467			`raise Argument_Error;`
468			`end if;`
469			`end Argument;`
470
471			`----------------------------`
472			`-- Compose_From_Cartesian --`
473			`----------------------------`
474
475			`function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is`
476			`begin`
477			`return (Re, Im);`
478			`end Compose_From_Cartesian;`
479
480			`function Compose_From_Cartesian (Re : Real'Base) return Complex is`
481			`begin`
482			`return (Re, 0.0);`
483			`end Compose_From_Cartesian;`
484
485			`function Compose_From_Cartesian (Im : Imaginary) return Complex is`
486			`begin`
487			`return (0.0, R (Im));`
488			`end Compose_From_Cartesian;`
489
490			`------------------------`
491			`-- Compose_From_Polar --`
492			`------------------------`
493
494			`function Compose_From_Polar (`
495			`Modulus, Argument : Real'Base)`
496			`return Complex`
497			`is`
498			`begin`
499			`if Modulus = 0.0 then`
500			`return (0.0, 0.0);`
501			`else`
502			`return (Modulus * R (Cos (Double (Argument))),`
503			`Modulus * R (Sin (Double (Argument))));`
504			`end if;`
505			`end Compose_From_Polar;`
506
507			`function Compose_From_Polar (`
508			`Modulus, Argument, Cycle : Real'Base)`
509			`return Complex`
510			`is`
511			`Arg : Real'Base;`
512
513			`begin`
514			`if Modulus = 0.0 then`
515			`return (0.0, 0.0);`
516
517			`elsif Cycle > 0.0 then`
518			`if Argument = 0.0 then`
519			`return (Modulus, 0.0);`
520
521			`elsif Argument = Cycle / 4.0 then`
522			`return (0.0, Modulus);`
523
524			`elsif Argument = Cycle / 2.0 then`
525			`return (-Modulus, 0.0);`
526
527			`elsif Argument = 3.0 * Cycle / R (4.0) then`
528			`return (0.0, -Modulus);`
529			`else`
530			`Arg := Two_Pi * Argument / Cycle;`
531			`return (Modulus * R (Cos (Double (Arg))),`
532			`Modulus * R (Sin (Double (Arg))));`
533			`end if;`
534			`else`
535			`raise Argument_Error;`
536			`end if;`
537			`end Compose_From_Polar;`
538
539			`---------------`
540			`-- Conjugate --`
541			`---------------`
542
543			`function Conjugate (X : Complex) return Complex is`
544			`begin`
545			`return Complex'(X.Re, -X.Im);`
546			`end Conjugate;`
547
548			`--------`
549			`-- Im --`
550			`--------`
551
552			`function Im (X : Complex) return Real'Base is`
553			`begin`
554			`return X.Im;`
555			`end Im;`
556
557			`function Im (X : Imaginary) return Real'Base is`
558			`begin`
559			`return R (X);`
560			`end Im;`
561
562			`-------------`
563			`-- Modulus --`
564			`-------------`
565
566			`function Modulus (X : Complex) return Real'Base is`
567			`Re2, Im2 : R;`
568
569			`begin`
570
571			`begin`
572			`Re2 := X.Re ** 2;`
573
574			`-- To compute (a2 + b2) (0.5) when a2 may be out of bounds,`
575			`-- compute a * (1 + (b/a) 2) (0.5). On a machine where the`
576			`-- squaring does not raise constraint_error but generates infinity,`
577			`-- we can use an explicit comparison to determine whether to use`
578			`-- the scaling expression.`
579
580			`-- The scaling expression is computed in double format throughout`
581			`-- in order to prevent inaccuracies on machines where not all`
582			`-- immediate expressions are rounded, such as PowerPC.`
583
584			`-- ??? same weird test, why not Re2 > R'Last ???`
585			`if not (Re2 <= R'Last) then`
586			`raise Constraint_Error;`
587			`end if;`
588
589			`exception`
590			`when Constraint_Error =>`
591			`return R (Double (abs (X.Re))`
592			`* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));`
593			`end;`
594
595			`begin`
596			`Im2 := X.Im ** 2;`
597
598			`-- ??? same weird test`
599			`if not (Im2 <= R'Last) then`
600			`raise Constraint_Error;`
601			`end if;`
602
603			`exception`
604			`when Constraint_Error =>`
605			`return R (Double (abs (X.Im))`
606			`* Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));`
607			`end;`
608
609			`-- Now deal with cases of underflow. If only one of the squares`
610			`-- underflows, return the modulus of the other component. If both`
611			`-- squares underflow, use scaling as above.`
612
613			`if Re2 = 0.0 then`
614
615			`if X.Re = 0.0 then`
616			`return abs (X.Im);`
617
618			`elsif Im2 = 0.0 then`
619
620			`if X.Im = 0.0 then`
621			`return abs (X.Re);`
622
623			`else`
624			`if abs (X.Re) > abs (X.Im) then`
625			`return`
626			`R (Double (abs (X.Re))`
627			`* Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));`
628			`else`
629			`return`
630			`R (Double (abs (X.Im))`
631			`* Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));`
632			`end if;`
633			`end if;`
634
635			`else`
636			`return abs (X.Im);`
637			`end if;`
638
639			`elsif Im2 = 0.0 then`
640			`return abs (X.Re);`
641
642			`-- In all other cases, the naive computation will do`
643
644			`else`
645			`return R (Sqrt (Double (Re2 + Im2)));`
646			`end if;`
647			`end Modulus;`
648
649			`--------`
650			`-- Re --`
651			`--------`
652
653			`function Re (X : Complex) return Real'Base is`
654			`begin`
655			`return X.Re;`
656			`end Re;`
657
658			`------------`
659			`-- Set_Im --`
660			`------------`
661
662			`procedure Set_Im (X : in out Complex; Im : Real'Base) is`
663			`begin`
664			`X.Im := Im;`
665			`end Set_Im;`
666
667			`procedure Set_Im (X : out Imaginary; Im : Real'Base) is`
668			`begin`
669			`X := Imaginary (Im);`
670			`end Set_Im;`
671
672			`------------`
673			`-- Set_Re --`
674			`------------`
675
676			`procedure Set_Re (X : in out Complex; Re : Real'Base) is`
677			`begin`
678			`X.Re := Re;`
679			`end Set_Re;`
680
681			`end Ada.Numerics.Generic_Complex_Types;`