From 554fd8c5195424bdbcabf5de30fdc183aba391bd Mon Sep 17 00:00:00 2001 From: upstream source tree Date: Sun, 15 Mar 2015 20:14:05 -0400 Subject: obtained gcc-4.6.4.tar.bz2 from upstream website; verified gcc-4.6.4.tar.bz2.sig; imported gcc-4.6.4 source tree from verified upstream tarball. downloading a git-generated archive based on the 'upstream' tag should provide you with a source tree that is binary identical to the one extracted from the above tarball. if you have obtained the source via the command 'git clone', however, do note that line-endings of files in your working directory might differ from line-endings of the respective files in the upstream repository. --- gcc/ada/a-ngcoty.adb | 681 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 681 insertions(+) create mode 100644 gcc/ada/a-ngcoty.adb (limited to 'gcc/ada/a-ngcoty.adb') diff --git a/gcc/ada/a-ngcoty.adb b/gcc/ada/a-ngcoty.adb new file mode 100644 index 000000000..7cf48713a --- /dev/null +++ b/gcc/ada/a-ngcoty.adb @@ -0,0 +1,681 @@ +------------------------------------------------------------------------------ +-- -- +-- 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 -- +-- . -- +-- -- +-- 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; -- cgit v1.2.3