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

Subversion Repositories openrisc

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

Go to most recent revision | Details | Compare with Previous | View Log


------------------------------------------------------------------------------
--                                                                          --
--                         GNAT RUN-TIME COMPONENTS                         --
--                                                                          --
--            ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS             --
--                                                                          --
--                                 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.Generic_Elementary_Functions;
 
package body Ada.Numerics.Generic_Complex_Elementary_Functions is
 
   package Elementary_Functions is new
      Ada.Numerics.Generic_Elementary_Functions (Real'Base);
   use Elementary_Functions;
 
   PI      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
   PI_2    : constant := PI / 2.0;
   Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
   Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
 
   subtype T is Real'Base;
 
   Epsilon                 : constant T := 2.0      ** (1 - T'Model_Mantissa);
   Square_Root_Epsilon     : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
   Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
   Root_Root_Epsilon       : constant T := Sqrt_Two **
                                                 ((1 - T'Model_Mantissa) / 2);
   Log_Inverse_Epsilon_2   : constant T := T (T'Model_Mantissa - 1) / 2.0;
 
   Complex_Zero : constant Complex := (0.0,  0.0);
   Complex_One  : constant Complex := (1.0,  0.0);
   Complex_I    : constant Complex := (0.0,  1.0);
   Half_Pi      : constant Complex := (PI_2, 0.0);
 
   --------
   -- ** --
   --------
 
   function "**" (Left : Complex; Right : Complex) return Complex is
   begin
      if Re (Right) = 0.0
        and then Im (Right) = 0.0
        and then Re (Left)  = 0.0
        and then Im (Left)  = 0.0
      then
         raise Argument_Error;
 
      elsif Re (Left) = 0.0
        and then Im (Left) = 0.0
        and then Re (Right) < 0.0
      then
         raise Constraint_Error;
 
      elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
         return Left;
 
      elsif Right = (0.0, 0.0)  then
         return Complex_One;
 
      elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
         return 1.0 + Right;
 
      elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
         return Left;
 
      else
         return Exp (Right * Log (Left));
      end if;
   end "**";
 
   function "**" (Left : Real'Base; Right : Complex) return Complex is
   begin
      if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
         raise Argument_Error;
 
      elsif Left = 0.0 and then Re (Right) < 0.0 then
         raise Constraint_Error;
 
      elsif Left = 0.0 then
         return Compose_From_Cartesian (Left, 0.0);
 
      elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
         return Complex_One;
 
      elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
         return Compose_From_Cartesian (Left, 0.0);
 
      else
         return Exp (Log (Left) * Right);
      end if;
   end "**";
 
   function "**" (Left : Complex; Right : Real'Base) return Complex is
   begin
      if Right = 0.0
        and then Re (Left) = 0.0
        and then Im (Left) = 0.0
      then
         raise Argument_Error;
 
      elsif Re (Left) = 0.0
        and then Im (Left) = 0.0
        and then Right < 0.0
      then
         raise Constraint_Error;
 
      elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
         return Left;
 
      elsif Right = 0.0 then
         return Complex_One;
 
      elsif Right = 1.0 then
         return Left;
 
      else
         return Exp (Right * Log (Left));
      end if;
   end "**";
 
   ------------
   -- Arccos --
   ------------
 
   function Arccos (X : Complex) return Complex is
      Result : Complex;
 
   begin
      if X = Complex_One then
         return Complex_Zero;
 
      elsif abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Half_Pi - X;
 
      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
                            Complex_I * Sqrt ((1.0 - X) / 2.0));
      end if;
 
      Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
 
      if Im (X) = 0.0
        and then abs Re (X) <= 1.00
      then
         Set_Im (Result, Im (X));
      end if;
 
      return Result;
   end Arccos;
 
   -------------
   -- Arccosh --
   -------------
 
   function Arccosh (X : Complex) return Complex is
      Result : Complex;
 
   begin
      if X = Complex_One then
         return Complex_Zero;
 
      elsif abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
 
      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := Log_Two + Log (X);
 
      else
         Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
                              Sqrt ((X - 1.0) / 2.0));
      end if;
 
      if Re (Result) <= 0.0 then
         Result := -Result;
      end if;
 
      return Result;
   end Arccosh;
 
   ------------
   -- Arccot --
   ------------
 
   function Arccot (X : Complex) return Complex is
      Xt : Complex;
 
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Half_Pi - X;
 
      elsif abs Re (X) > 1.0 / Epsilon or else
            abs Im (X) > 1.0 / Epsilon
      then
         Xt := Complex_One  /  X;
 
         if Re (X) < 0.0 then
            Set_Re (Xt, PI - Re (Xt));
            return Xt;
         else
            return Xt;
         end if;
      end if;
 
      Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
 
      if Re (Xt) < 0.0 then
         Xt := PI + Xt;
      end if;
 
      return Xt;
   end Arccot;
 
   --------------
   -- Arccoth --
   --------------
 
   function Arccoth (X : Complex) return Complex is
      R : Complex;
 
   begin
      if X = (0.0, 0.0) then
         return Compose_From_Cartesian (0.0, PI_2);
 
      elsif abs Re (X) < Square_Root_Epsilon
         and then abs Im (X) < Square_Root_Epsilon
      then
         return PI_2 * Complex_I + X;
 
      elsif abs Re (X) > 1.0 / Epsilon or else
            abs Im (X) > 1.0 / Epsilon
      then
         if Im (X) > 0.0 then
            return (0.0, 0.0);
         else
            return PI * Complex_I;
         end if;
 
      elsif Im (X) = 0.0 and then Re (X) = 1.0 then
         raise Constraint_Error;
 
      elsif Im (X) = 0.0 and then Re (X) = -1.0 then
         raise Constraint_Error;
      end if;
 
      begin
         R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
 
      exception
         when Constraint_Error =>
            R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
      end;
 
      if Im (R) < 0.0 then
         Set_Im (R, PI + Im (R));
      end if;
 
      if Re (X) = 0.0 then
         Set_Re (R, Re (X));
      end if;
 
      return R;
   end Arccoth;
 
   ------------
   -- Arcsin --
   ------------
 
   function Arcsin (X : Complex) return Complex is
      Result : Complex;
 
   begin
      --  For very small argument, sin (x) = x
 
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
 
         if Im (Result) > PI_2 then
            Set_Im (Result, PI - Im (X));
 
         elsif Im (Result) < -PI_2 then
            Set_Im (Result, -(PI + Im (X)));
         end if;
 
         return Result;
      end if;
 
      Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
 
      if Re (X) = 0.0 then
         Set_Re (Result, Re (X));
 
      elsif Im (X) = 0.0
        and then abs Re (X) <= 1.00
      then
         Set_Im (Result, Im (X));
      end if;
 
      return Result;
   end Arcsin;
 
   -------------
   -- Arcsinh --
   -------------
 
   function Arcsinh (X : Complex) return Complex is
      Result : Complex;
 
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      elsif abs Re (X) > Inv_Square_Root_Epsilon or else
            abs Im (X) > Inv_Square_Root_Epsilon
      then
         Result := Log_Two + Log (X); -- may have wrong sign
 
         if (Re (X) < 0.0 and then Re (Result) > 0.0)
           or else (Re (X) > 0.0 and then Re (Result) < 0.0)
         then
            Set_Re (Result, -Re (Result));
         end if;
 
         return Result;
      end if;
 
      Result := Log (X + Sqrt (1.0 + X * X));
 
      if Re (X) = 0.0 then
         Set_Re (Result, Re (X));
      elsif Im  (X) = 0.0 then
         Set_Im (Result, Im  (X));
      end if;
 
      return Result;
   end Arcsinh;
 
   ------------
   -- Arctan --
   ------------
 
   function Arctan (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      else
         return -Complex_I * (Log (1.0 + Complex_I * X)
                            - Log (1.0 - Complex_I * X)) / 2.0;
      end if;
   end Arctan;
 
   -------------
   -- Arctanh --
   -------------
 
   function Arctanh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
      else
         return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
      end if;
   end Arctanh;
 
   ---------
   -- Cos --
   ---------
 
   function Cos (X : Complex) return Complex is
   begin
      return
        Compose_From_Cartesian
          (Cos (Re (X))  * Cosh (Im (X)),
           -(Sin (Re (X)) * Sinh (Im (X))));
   end Cos;
 
   ----------
   -- Cosh --
   ----------
 
   function Cosh (X : Complex) return Complex is
   begin
      return
        Compose_From_Cartesian
          (Cosh (Re (X)) * Cos (Im (X)),
           Sinh (Re (X)) * Sin (Im (X)));
   end Cosh;
 
   ---------
   -- Cot --
   ---------
 
   function Cot (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Complex_One  /  X;
 
      elsif Im (X) > Log_Inverse_Epsilon_2 then
         return -Complex_I;
 
      elsif Im (X) < -Log_Inverse_Epsilon_2 then
         return Complex_I;
      end if;
 
      return Cos (X) / Sin (X);
   end Cot;
 
   ----------
   -- Coth --
   ----------
 
   function Coth (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return Complex_One  /  X;
 
      elsif Re (X) > Log_Inverse_Epsilon_2 then
         return Complex_One;
 
      elsif Re (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_One;
 
      else
         return Cosh (X) / Sinh (X);
      end if;
   end Coth;
 
   ---------
   -- Exp --
   ---------
 
   function Exp (X : Complex) return Complex is
      EXP_RE_X : constant Real'Base := Exp (Re (X));
 
   begin
      return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
                                     EXP_RE_X * Sin (Im (X)));
   end Exp;
 
   function Exp (X : Imaginary) return Complex is
      ImX : constant Real'Base := Im (X);
 
   begin
      return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
   end Exp;
 
   ---------
   -- Log --
   ---------
 
   function Log (X : Complex) return Complex is
      ReX : Real'Base;
      ImX : Real'Base;
      Z   : Complex;
 
   begin
      if Re (X) = 0.0 and then Im (X) = 0.0 then
         raise Constraint_Error;
 
      elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
        and then abs Im (X) < Root_Root_Epsilon
      then
         Z := X;
         Set_Re (Z, Re (Z) - 1.0);
 
         return (1.0 - (1.0 / 2.0 -
                       (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
      end if;
 
      begin
         ReX := Log (Modulus (X));
 
      exception
         when Constraint_Error =>
            ReX := Log (Modulus (X / 2.0)) - Log_Two;
      end;
 
      ImX := Arctan (Im (X), Re (X));
 
      if ImX > PI then
         ImX := ImX - 2.0 * PI;
      end if;
 
      return Compose_From_Cartesian (ReX, ImX);
   end Log;
 
   ---------
   -- Sin --
   ---------
 
   function Sin (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon then
         return X;
      end if;
 
      return
        Compose_From_Cartesian
          (Sin (Re (X)) * Cosh (Im (X)),
           Cos (Re (X)) * Sinh (Im (X)));
   end Sin;
 
   ----------
   -- Sinh --
   ----------
 
   function Sinh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      else
         return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
                                        Cosh (Re (X)) * Sin (Im (X)));
      end if;
   end Sinh;
 
   ----------
   -- Sqrt --
   ----------
 
   function Sqrt (X : Complex) return Complex is
      ReX : constant Real'Base := Re (X);
      ImX : constant Real'Base := Im (X);
      XR  : constant Real'Base := abs Re (X);
      YR  : constant Real'Base := abs Im (X);
      R   : Real'Base;
      R_X : Real'Base;
      R_Y : Real'Base;
 
   begin
      --  Deal with pure real case, see (RM G.1.2(39))
 
      if ImX = 0.0 then
         if ReX > 0.0 then
            return
              Compose_From_Cartesian
                (Sqrt (ReX), 0.0);
 
         elsif ReX = 0.0 then
            return X;
 
         else
            return
              Compose_From_Cartesian
                (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
         end if;
 
      elsif ReX = 0.0 then
         R_X := Sqrt (YR / 2.0);
 
         if ImX > 0.0 then
            return Compose_From_Cartesian (R_X, R_X);
         else
            return Compose_From_Cartesian (R_X, -R_X);
         end if;
 
      else
         R  := Sqrt (XR ** 2 + YR ** 2);
 
         --  If the square of the modulus overflows, try rescaling the
         --  real and imaginary parts. We cannot depend on an exception
         --  being raised on all targets.
 
         if R > Real'Base'Last then
            raise Constraint_Error;
         end if;
 
         --  We are solving the system
 
         --  XR = R_X ** 2 - Y_R ** 2      (1)
         --  YR = 2.0 * R_X * R_Y          (2)
         --
         --  The symmetric solution involves square roots for both R_X and
         --  R_Y, but it is more accurate to use the square root with the
         --  larger argument for either R_X or R_Y, and equation (2) for the
         --  other.
 
         if ReX < 0.0 then
            R_Y := Sqrt (0.5 * (R - ReX));
            R_X := YR / (2.0 * R_Y);
 
         else
            R_X := Sqrt (0.5 * (R + ReX));
            R_Y := YR / (2.0 * R_X);
         end if;
      end if;
 
      if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
         R_Y := -R_Y;
      end if;
      return Compose_From_Cartesian (R_X, R_Y);
 
   exception
      when Constraint_Error =>
 
         --  Rescale and try again
 
         R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
         R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
         R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
 
         if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
            R_Y := -R_Y;
         end if;
 
         return Compose_From_Cartesian (R_X, R_Y);
   end Sqrt;
 
   ---------
   -- Tan --
   ---------
 
   function Tan (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      elsif Im (X) > Log_Inverse_Epsilon_2 then
         return Complex_I;
 
      elsif Im (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_I;
 
      else
         return Sin (X) / Cos (X);
      end if;
   end Tan;
 
   ----------
   -- Tanh --
   ----------
 
   function Tanh (X : Complex) return Complex is
   begin
      if abs Re (X) < Square_Root_Epsilon and then
         abs Im (X) < Square_Root_Epsilon
      then
         return X;
 
      elsif Re (X) > Log_Inverse_Epsilon_2 then
         return Complex_One;
 
      elsif Re (X) < -Log_Inverse_Epsilon_2 then
         return -Complex_One;
 
      else
         return Sinh (X) / Cosh (X);
      end if;
   end Tanh;
 
end Ada.Numerics.Generic_Complex_Elementary_Functions;

Browse

Tools

Subversion Repositories openrisc

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

Line No.	Rev	Author	Line
1	706	jeremybenn	`------------------------------------------------------------------------------`
2			`-- --`
3			`-- GNAT RUN-TIME COMPONENTS --`
4			`-- --`
5			`-- ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS --`
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.Generic_Elementary_Functions;`
33
34			`package body Ada.Numerics.Generic_Complex_Elementary_Functions is`
35
36			`package Elementary_Functions is new`
37			`Ada.Numerics.Generic_Elementary_Functions (Real'Base);`
38			`use Elementary_Functions;`
39
40			`PI : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;`
41			`PI_2 : constant := PI / 2.0;`
42			`Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;`
43			`Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;`
44
45			`subtype T is Real'Base;`
46
47			`Epsilon : constant T := 2.0 ** (1 - T'Model_Mantissa);`
48			`Square_Root_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);`
49			`Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);`
50			`Root_Root_Epsilon : constant T := Sqrt_Two **`
51			`((1 - T'Model_Mantissa) / 2);`
52			`Log_Inverse_Epsilon_2 : constant T := T (T'Model_Mantissa - 1) / 2.0;`
53
54			`Complex_Zero : constant Complex := (0.0, 0.0);`
55			`Complex_One : constant Complex := (1.0, 0.0);`
56			`Complex_I : constant Complex := (0.0, 1.0);`
57			`Half_Pi : constant Complex := (PI_2, 0.0);`
58
59			`--------`
60			`-- ** --`
61			`--------`
62
63			`function "**" (Left : Complex; Right : Complex) return Complex is`
64			`begin`
65			`if Re (Right) = 0.0`
66			`and then Im (Right) = 0.0`
67			`and then Re (Left) = 0.0`
68			`and then Im (Left) = 0.0`
69			`then`
70			`raise Argument_Error;`
71
72			`elsif Re (Left) = 0.0`
73			`and then Im (Left) = 0.0`
74			`and then Re (Right) < 0.0`
75			`then`
76			`raise Constraint_Error;`
77
78			`elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then`
79			`return Left;`
80
81			`elsif Right = (0.0, 0.0) then`
82			`return Complex_One;`
83
84			`elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then`
85			`return 1.0 + Right;`
86
87			`elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then`
88			`return Left;`
89
90			`else`
91			`return Exp (Right * Log (Left));`
92			`end if;`
93			`end "**";`
94
95			`function "**" (Left : Real'Base; Right : Complex) return Complex is`
96			`begin`
97			`if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then`
98			`raise Argument_Error;`
99
100			`elsif Left = 0.0 and then Re (Right) < 0.0 then`
101			`raise Constraint_Error;`
102
103			`elsif Left = 0.0 then`
104			`return Compose_From_Cartesian (Left, 0.0);`
105
106			`elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then`
107			`return Complex_One;`
108
109			`elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then`
110			`return Compose_From_Cartesian (Left, 0.0);`
111
112			`else`
113			`return Exp (Log (Left) * Right);`
114			`end if;`
115			`end "**";`
116
117			`function "**" (Left : Complex; Right : Real'Base) return Complex is`
118			`begin`
119			`if Right = 0.0`
120			`and then Re (Left) = 0.0`
121			`and then Im (Left) = 0.0`
122			`then`
123			`raise Argument_Error;`
124
125			`elsif Re (Left) = 0.0`
126			`and then Im (Left) = 0.0`
127			`and then Right < 0.0`
128			`then`
129			`raise Constraint_Error;`
130
131			`elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then`
132			`return Left;`
133
134			`elsif Right = 0.0 then`
135			`return Complex_One;`
136
137			`elsif Right = 1.0 then`
138			`return Left;`
139
140			`else`
141			`return Exp (Right * Log (Left));`
142			`end if;`
143			`end "**";`
144
145			`------------`
146			`-- Arccos --`
147			`------------`
148
149			`function Arccos (X : Complex) return Complex is`
150			`Result : Complex;`
151
152			`begin`
153			`if X = Complex_One then`
154			`return Complex_Zero;`
155
156			`elsif abs Re (X) < Square_Root_Epsilon and then`
157			`abs Im (X) < Square_Root_Epsilon`
158			`then`
159			`return Half_Pi - X;`
160
161			`elsif abs Re (X) > Inv_Square_Root_Epsilon or else`
162			`abs Im (X) > Inv_Square_Root_Epsilon`
163			`then`
164			`return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +`
165			`Complex_I * Sqrt ((1.0 - X) / 2.0));`
166			`end if;`
167
168			`Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));`
169
170			`if Im (X) = 0.0`
171			`and then abs Re (X) <= 1.00`
172			`then`
173			`Set_Im (Result, Im (X));`
174			`end if;`
175
176			`return Result;`
177			`end Arccos;`
178
179			`-------------`
180			`-- Arccosh --`
181			`-------------`
182
183			`function Arccosh (X : Complex) return Complex is`
184			`Result : Complex;`
185
186			`begin`
187			`if X = Complex_One then`
188			`return Complex_Zero;`
189
190			`elsif abs Re (X) < Square_Root_Epsilon and then`
191			`abs Im (X) < Square_Root_Epsilon`
192			`then`
193			`Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));`
194
195			`elsif abs Re (X) > Inv_Square_Root_Epsilon or else`
196			`abs Im (X) > Inv_Square_Root_Epsilon`
197			`then`
198			`Result := Log_Two + Log (X);`
199
200			`else`
201			`Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +`
202			`Sqrt ((X - 1.0) / 2.0));`
203			`end if;`
204
205			`if Re (Result) <= 0.0 then`
206			`Result := -Result;`
207			`end if;`
208
209			`return Result;`
210			`end Arccosh;`
211
212			`------------`
213			`-- Arccot --`
214			`------------`
215
216			`function Arccot (X : Complex) return Complex is`
217			`Xt : Complex;`
218
219			`begin`
220			`if abs Re (X) < Square_Root_Epsilon and then`
221			`abs Im (X) < Square_Root_Epsilon`
222			`then`
223			`return Half_Pi - X;`
224
225			`elsif abs Re (X) > 1.0 / Epsilon or else`
226			`abs Im (X) > 1.0 / Epsilon`
227			`then`
228			`Xt := Complex_One / X;`
229
230			`if Re (X) < 0.0 then`
231			`Set_Re (Xt, PI - Re (Xt));`
232			`return Xt;`
233			`else`
234			`return Xt;`
235			`end if;`
236			`end if;`
237
238			`Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;`
239
240			`if Re (Xt) < 0.0 then`
241			`Xt := PI + Xt;`
242			`end if;`
243
244			`return Xt;`
245			`end Arccot;`
246
247			`--------------`
248			`-- Arccoth --`
249			`--------------`
250
251			`function Arccoth (X : Complex) return Complex is`
252			`R : Complex;`
253
254			`begin`
255			`if X = (0.0, 0.0) then`
256			`return Compose_From_Cartesian (0.0, PI_2);`
257
258			`elsif abs Re (X) < Square_Root_Epsilon`
259			`and then abs Im (X) < Square_Root_Epsilon`
260			`then`
261			`return PI_2 * Complex_I + X;`
262
263			`elsif abs Re (X) > 1.0 / Epsilon or else`
264			`abs Im (X) > 1.0 / Epsilon`
265			`then`
266			`if Im (X) > 0.0 then`
267			`return (0.0, 0.0);`
268			`else`
269			`return PI * Complex_I;`
270			`end if;`
271
272			`elsif Im (X) = 0.0 and then Re (X) = 1.0 then`
273			`raise Constraint_Error;`
274
275			`elsif Im (X) = 0.0 and then Re (X) = -1.0 then`
276			`raise Constraint_Error;`
277			`end if;`
278
279			`begin`
280			`R := Log ((1.0 + X) / (X - 1.0)) / 2.0;`
281
282			`exception`
283			`when Constraint_Error =>`
284			`R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;`
285			`end;`
286
287			`if Im (R) < 0.0 then`
288			`Set_Im (R, PI + Im (R));`
289			`end if;`
290
291			`if Re (X) = 0.0 then`
292			`Set_Re (R, Re (X));`
293			`end if;`
294
295			`return R;`
296			`end Arccoth;`
297
298			`------------`
299			`-- Arcsin --`
300			`------------`
301
302			`function Arcsin (X : Complex) return Complex is`
303			`Result : Complex;`
304
305			`begin`
306			`-- For very small argument, sin (x) = x`
307
308			`if abs Re (X) < Square_Root_Epsilon and then`
309			`abs Im (X) < Square_Root_Epsilon`
310			`then`
311			`return X;`
312
313			`elsif abs Re (X) > Inv_Square_Root_Epsilon or else`
314			`abs Im (X) > Inv_Square_Root_Epsilon`
315			`then`
316			`Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));`
317
318			`if Im (Result) > PI_2 then`
319			`Set_Im (Result, PI - Im (X));`
320
321			`elsif Im (Result) < -PI_2 then`
322			`Set_Im (Result, -(PI + Im (X)));`
323			`end if;`
324
325			`return Result;`
326			`end if;`
327
328			`Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));`
329
330			`if Re (X) = 0.0 then`
331			`Set_Re (Result, Re (X));`
332
333			`elsif Im (X) = 0.0`
334			`and then abs Re (X) <= 1.00`
335			`then`
336			`Set_Im (Result, Im (X));`
337			`end if;`
338
339			`return Result;`
340			`end Arcsin;`
341
342			`-------------`
343			`-- Arcsinh --`
344			`-------------`
345
346			`function Arcsinh (X : Complex) return Complex is`
347			`Result : Complex;`
348
349			`begin`
350			`if abs Re (X) < Square_Root_Epsilon and then`
351			`abs Im (X) < Square_Root_Epsilon`
352			`then`
353			`return X;`
354
355			`elsif abs Re (X) > Inv_Square_Root_Epsilon or else`
356			`abs Im (X) > Inv_Square_Root_Epsilon`
357			`then`
358			`Result := Log_Two + Log (X); -- may have wrong sign`
359
360			`if (Re (X) < 0.0 and then Re (Result) > 0.0)`
361			`or else (Re (X) > 0.0 and then Re (Result) < 0.0)`
362			`then`
363			`Set_Re (Result, -Re (Result));`
364			`end if;`
365
366			`return Result;`
367			`end if;`
368
369			`Result := Log (X + Sqrt (1.0 + X * X));`
370
371			`if Re (X) = 0.0 then`
372			`Set_Re (Result, Re (X));`
373			`elsif Im (X) = 0.0 then`
374			`Set_Im (Result, Im (X));`
375			`end if;`
376
377			`return Result;`
378			`end Arcsinh;`
379
380			`------------`
381			`-- Arctan --`
382			`------------`
383
384			`function Arctan (X : Complex) return Complex is`
385			`begin`
386			`if abs Re (X) < Square_Root_Epsilon and then`
387			`abs Im (X) < Square_Root_Epsilon`
388			`then`
389			`return X;`
390
391			`else`
392			`return -Complex_I * (Log (1.0 + Complex_I * X)`
393			`- Log (1.0 - Complex_I * X)) / 2.0;`
394			`end if;`
395			`end Arctan;`
396
397			`-------------`
398			`-- Arctanh --`
399			`-------------`
400
401			`function Arctanh (X : Complex) return Complex is`
402			`begin`
403			`if abs Re (X) < Square_Root_Epsilon and then`
404			`abs Im (X) < Square_Root_Epsilon`
405			`then`
406			`return X;`
407			`else`
408			`return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;`
409			`end if;`
410			`end Arctanh;`
411
412			`---------`
413			`-- Cos --`
414			`---------`
415
416			`function Cos (X : Complex) return Complex is`
417			`begin`
418			`return`
419			`Compose_From_Cartesian`
420			`(Cos (Re (X)) * Cosh (Im (X)),`
421			`-(Sin (Re (X)) * Sinh (Im (X))));`
422			`end Cos;`
423
424			`----------`
425			`-- Cosh --`
426			`----------`
427
428			`function Cosh (X : Complex) return Complex is`
429			`begin`
430			`return`
431			`Compose_From_Cartesian`
432			`(Cosh (Re (X)) * Cos (Im (X)),`
433			`Sinh (Re (X)) * Sin (Im (X)));`
434			`end Cosh;`
435
436			`---------`
437			`-- Cot --`
438			`---------`
439
440			`function Cot (X : Complex) return Complex is`
441			`begin`
442			`if abs Re (X) < Square_Root_Epsilon and then`
443			`abs Im (X) < Square_Root_Epsilon`
444			`then`
445			`return Complex_One / X;`
446
447			`elsif Im (X) > Log_Inverse_Epsilon_2 then`
448			`return -Complex_I;`
449
450			`elsif Im (X) < -Log_Inverse_Epsilon_2 then`
451			`return Complex_I;`
452			`end if;`
453
454			`return Cos (X) / Sin (X);`
455			`end Cot;`
456
457			`----------`
458			`-- Coth --`
459			`----------`
460
461			`function Coth (X : Complex) return Complex is`
462			`begin`
463			`if abs Re (X) < Square_Root_Epsilon and then`
464			`abs Im (X) < Square_Root_Epsilon`
465			`then`
466			`return Complex_One / X;`
467
468			`elsif Re (X) > Log_Inverse_Epsilon_2 then`
469			`return Complex_One;`
470
471			`elsif Re (X) < -Log_Inverse_Epsilon_2 then`
472			`return -Complex_One;`
473
474			`else`
475			`return Cosh (X) / Sinh (X);`
476			`end if;`
477			`end Coth;`
478
479			`---------`
480			`-- Exp --`
481			`---------`
482
483			`function Exp (X : Complex) return Complex is`
484			`EXP_RE_X : constant Real'Base := Exp (Re (X));`
485
486			`begin`
487			`return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),`
488			`EXP_RE_X * Sin (Im (X)));`
489			`end Exp;`
490
491			`function Exp (X : Imaginary) return Complex is`
492			`ImX : constant Real'Base := Im (X);`
493
494			`begin`
495			`return Compose_From_Cartesian (Cos (ImX), Sin (ImX));`
496			`end Exp;`
497
498			`---------`
499			`-- Log --`
500			`---------`
501
502			`function Log (X : Complex) return Complex is`
503			`ReX : Real'Base;`
504			`ImX : Real'Base;`
505			`Z : Complex;`
506
507			`begin`
508			`if Re (X) = 0.0 and then Im (X) = 0.0 then`
509			`raise Constraint_Error;`
510
511			`elsif abs (1.0 - Re (X)) < Root_Root_Epsilon`
512			`and then abs Im (X) < Root_Root_Epsilon`
513			`then`
514			`Z := X;`
515			`Set_Re (Z, Re (Z) - 1.0);`
516
517			`return (1.0 - (1.0 / 2.0 -`
518			`(1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;`
519			`end if;`
520
521			`begin`
522			`ReX := Log (Modulus (X));`
523
524			`exception`
525			`when Constraint_Error =>`
526			`ReX := Log (Modulus (X / 2.0)) - Log_Two;`
527			`end;`
528
529			`ImX := Arctan (Im (X), Re (X));`
530
531			`if ImX > PI then`
532			`ImX := ImX - 2.0 * PI;`
533			`end if;`
534
535			`return Compose_From_Cartesian (ReX, ImX);`
536			`end Log;`
537
538			`---------`
539			`-- Sin --`
540			`---------`
541
542			`function Sin (X : Complex) return Complex is`
543			`begin`
544			`if abs Re (X) < Square_Root_Epsilon and then`
545			`abs Im (X) < Square_Root_Epsilon then`
546			`return X;`
547			`end if;`
548
549			`return`
550			`Compose_From_Cartesian`
551			`(Sin (Re (X)) * Cosh (Im (X)),`
552			`Cos (Re (X)) * Sinh (Im (X)));`
553			`end Sin;`
554
555			`----------`
556			`-- Sinh --`
557			`----------`
558
559			`function Sinh (X : Complex) return Complex is`
560			`begin`
561			`if abs Re (X) < Square_Root_Epsilon and then`
562			`abs Im (X) < Square_Root_Epsilon`
563			`then`
564			`return X;`
565
566			`else`
567			`return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),`
568			`Cosh (Re (X)) * Sin (Im (X)));`
569			`end if;`
570			`end Sinh;`
571
572			`----------`
573			`-- Sqrt --`
574			`----------`
575
576			`function Sqrt (X : Complex) return Complex is`
577			`ReX : constant Real'Base := Re (X);`
578			`ImX : constant Real'Base := Im (X);`
579			`XR : constant Real'Base := abs Re (X);`
580			`YR : constant Real'Base := abs Im (X);`
581			`R : Real'Base;`
582			`R_X : Real'Base;`
583			`R_Y : Real'Base;`
584
585			`begin`
586			`-- Deal with pure real case, see (RM G.1.2(39))`
587
588			`if ImX = 0.0 then`
589			`if ReX > 0.0 then`
590			`return`
591			`Compose_From_Cartesian`
592			`(Sqrt (ReX), 0.0);`
593
594			`elsif ReX = 0.0 then`
595			`return X;`
596
597			`else`
598			`return`
599			`Compose_From_Cartesian`
600			`(0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));`
601			`end if;`
602
603			`elsif ReX = 0.0 then`
604			`R_X := Sqrt (YR / 2.0);`
605
606			`if ImX > 0.0 then`
607			`return Compose_From_Cartesian (R_X, R_X);`
608			`else`
609			`return Compose_From_Cartesian (R_X, -R_X);`
610			`end if;`
611
612			`else`
613			`R := Sqrt (XR 2 + YR 2);`
614
615			`-- If the square of the modulus overflows, try rescaling the`
616			`-- real and imaginary parts. We cannot depend on an exception`
617			`-- being raised on all targets.`
618
619			`if R > Real'Base'Last then`
620			`raise Constraint_Error;`
621			`end if;`
622
623			`-- We are solving the system`
624
625			`-- XR = R_X 2 - Y_R 2 (1)`
626			`-- YR = 2.0 * R_X * R_Y (2)`
627			`--`
628			`-- The symmetric solution involves square roots for both R_X and`
629			`-- R_Y, but it is more accurate to use the square root with the`
630			`-- larger argument for either R_X or R_Y, and equation (2) for the`
631			`-- other.`
632
633			`if ReX < 0.0 then`
634			`R_Y := Sqrt (0.5 * (R - ReX));`
635			`R_X := YR / (2.0 * R_Y);`
636
637			`else`
638			`R_X := Sqrt (0.5 * (R + ReX));`
639			`R_Y := YR / (2.0 * R_X);`
640			`end if;`
641			`end if;`
642
643			`if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude`
644			`R_Y := -R_Y;`
645			`end if;`
646			`return Compose_From_Cartesian (R_X, R_Y);`
647
648			`exception`
649			`when Constraint_Error =>`
650
651			`-- Rescale and try again`
652
653			`R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));`
654			`R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));`
655			`R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));`
656
657			`if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude`
658			`R_Y := -R_Y;`
659			`end if;`
660
661			`return Compose_From_Cartesian (R_X, R_Y);`
662			`end Sqrt;`
663
664			`---------`
665			`-- Tan --`
666			`---------`
667
668			`function Tan (X : Complex) return Complex is`
669			`begin`
670			`if abs Re (X) < Square_Root_Epsilon and then`
671			`abs Im (X) < Square_Root_Epsilon`
672			`then`
673			`return X;`
674
675			`elsif Im (X) > Log_Inverse_Epsilon_2 then`
676			`return Complex_I;`
677
678			`elsif Im (X) < -Log_Inverse_Epsilon_2 then`
679			`return -Complex_I;`
680
681			`else`
682			`return Sin (X) / Cos (X);`
683			`end if;`
684			`end Tan;`
685
686			`----------`
687			`-- Tanh --`
688			`----------`
689
690			`function Tanh (X : Complex) return Complex is`
691			`begin`
692			`if abs Re (X) < Square_Root_Epsilon and then`
693			`abs Im (X) < Square_Root_Epsilon`
694			`then`
695			`return X;`
696
697			`elsif Re (X) > Log_Inverse_Epsilon_2 then`
698			`return Complex_One;`
699
700			`elsif Re (X) < -Log_Inverse_Epsilon_2 then`
701			`return -Complex_One;`
702
703			`else`
704			`return Sinh (X) / Cosh (X);`
705			`end if;`
706			`end Tanh;`
707
708			`end Ada.Numerics.Generic_Complex_Elementary_Functions;`