gcc/ada/a-ngcoty.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUNTIME 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 --                            $Revision: 1.16 $
  10 --                                                                          --
  11 --          Copyright (C) 1992-2001 Free Software Foundation, Inc.          --
  12 --                                                                          --
  13 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14 -- terms of the  GNU General Public License as published  by the Free Soft- --
  15 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19 -- for  more details.  You should have  received  a copy of the GNU General --
  20 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  22 -- MA 02111-1307, USA.                                                      --
  23 --                                                                          --
  24 -- As a special exception,  if other files  instantiate  generics from this --
  25 -- unit, or you link  this unit with other files  to produce an executable, --
  26 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  27 -- covered  by the  GNU  General  Public  License.  This exception does not --
  28 -- however invalidate  any other reasons why  the executable file  might be --
  29 -- covered by the  GNU Public License.                                      --
  30 --                                                                          --
  31 -- GNAT was originally developed  by the GNAT team at  New York University. --
  32 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
  33 --                                                                          --
  34 ------------------------------------------------------------------------------
  35
  36 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
  37 package body Ada.Numerics.Generic_Complex_Types is
  38
  39    subtype R is Real'Base;
  40
  41    Two_Pi  : constant R := R (2.0) * Pi;
  42    Half_Pi : constant R := Pi / R (2.0);
  43
  44    ---------
  45    -- "*" --
  46    ---------
  47
  48    function "*" (Left, Right : Complex) return Complex is
  49       X : R;
  50       Y : R;
  51
  52    begin
  53       X := Left.Re * Right.Re - Left.Im * Right.Im;
  54       Y := Left.Re * Right.Im + Left.Im * Right.Re;
  55
  56       --  If either component overflows, try to scale.
  57
  58       if abs (X) > R'Last then
  59          X := R' (4.0) * (R'(Left.Re / 2.0)  * R'(Right.Re / 2.0)
  60                 - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
  61       end if;
  62
  63       if abs (Y) > R'Last then
  64          Y := R' (4.0) * (R'(Left.Re / 2.0)  * R'(Right.Im / 2.0)
  65                 - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
  66       end if;
  67
  68       return (X, Y);
  69    end "*";
  70
  71    function "*" (Left, Right : Imaginary) return Real'Base is
  72    begin
  73       return -R (Left) * R (Right);
  74    end "*";
  75
  76    function "*" (Left : Complex; Right : Real'Base) return Complex is
  77    begin
  78       return Complex'(Left.Re * Right, Left.Im * Right);
  79    end "*";
  80
  81    function "*" (Left : Real'Base; Right : Complex) return Complex is
  82    begin
  83       return (Left * Right.Re, Left * Right.Im);
  84    end "*";
  85
  86    function "*" (Left : Complex; Right : Imaginary) return Complex is
  87    begin
  88       return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
  89    end "*";
  90
  91    function "*" (Left : Imaginary; Right : Complex) return Complex is
  92    begin
  93       return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
  94    end "*";
  95
  96    function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
  97    begin
  98       return Left * Imaginary (Right);
  99    end "*";
 100
 101    function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
 102    begin
 103       return Imaginary (Left * R (Right));
 104    end "*";
 105
 106    ----------
 107    -- "**" --
 108    ----------
 109
 110    function "**" (Left : Complex; Right : Integer) return Complex is
 111       Result : Complex := (1.0, 0.0);
 112       Factor : Complex := Left;
 113       Exp    : Integer := Right;
 114
 115    begin
 116       --  We use the standard logarithmic approach, Exp gets shifted right
 117       --  testing successive low order bits and Factor is the value of the
 118       --  base raised to the next power of 2. For positive exponents we
 119       --  multiply the result by this factor, for negative exponents, we
 120       --  divide by this factor.
 121
 122       if Exp >= 0 then
 123
 124          --  For a positive exponent, if we get a constraint error during
 125          --  this loop, it is an overflow, and the constraint error will
 126          --  simply be passed on to the caller.
 127
 128          while Exp /= 0 loop
 129             if Exp rem 2 /= 0 then
 130                Result := Result * Factor;
 131             end if;
 132
 133             Factor := Factor * Factor;
 134             Exp := Exp / 2;
 135          end loop;
 136
 137          return Result;
 138
 139       else -- Exp < 0 then
 140
 141          --  For the negative exponent case, a constraint error during this
 142          --  calculation happens if Factor gets too large, and the proper
 143          --  response is to return 0.0, since what we essentially have is
 144          --  1.0 / infinity, and the closest model number will be zero.
 145
 146          begin
 147
 148             while Exp /= 0 loop
 149                if Exp rem 2 /= 0 then
 150                   Result := Result * Factor;
 151                end if;
 152
 153                Factor := Factor * Factor;
 154                Exp := Exp / 2;
 155             end loop;
 156
 157             return R ' (1.0) / Result;
 158
 159          exception
 160
 161             when Constraint_Error =>
 162                return (0.0, 0.0);
 163          end;
 164       end if;
 165    end "**";
 166
 167    function "**" (Left : Imaginary; Right : Integer) return Complex is
 168       M : R := R (Left) ** Right;
 169    begin
 170       case Right mod 4 is
 171          when 0 => return (M,   0.0);
 172          when 1 => return (0.0, M);
 173          when 2 => return (-M,  0.0);
 174          when 3 => return (0.0, -M);
 175          when others => raise Program_Error;
 176       end case;
 177    end "**";
 178
 179    ---------
 180    -- "+" --
 181    ---------
 182
 183    function "+" (Right : Complex) return Complex is
 184    begin
 185       return Right;
 186    end "+";
 187
 188    function "+" (Left, Right : Complex) return Complex is
 189    begin
 190       return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
 191    end "+";
 192
 193    function "+" (Right : Imaginary) return Imaginary is
 194    begin
 195       return Right;
 196    end "+";
 197
 198    function "+" (Left, Right : Imaginary) return Imaginary is
 199    begin
 200       return Imaginary (R (Left) + R (Right));
 201    end "+";
 202
 203    function "+" (Left : Complex; Right : Real'Base) return Complex is
 204    begin
 205       return Complex'(Left.Re + Right, Left.Im);
 206    end "+";
 207
 208    function "+" (Left : Real'Base; Right : Complex) return Complex is
 209    begin
 210       return Complex'(Left + Right.Re, Right.Im);
 211    end "+";
 212
 213    function "+" (Left : Complex; Right : Imaginary) return Complex is
 214    begin
 215       return Complex'(Left.Re, Left.Im + R (Right));
 216    end "+";
 217
 218    function "+" (Left : Imaginary; Right : Complex) return Complex is
 219    begin
 220       return Complex'(Right.Re, R (Left) + Right.Im);
 221    end "+";
 222
 223    function "+" (Left : Imaginary; Right : Real'Base) return Complex is
 224    begin
 225       return Complex'(Right, R (Left));
 226    end "+";
 227
 228    function "+" (Left : Real'Base; Right : Imaginary) return Complex is
 229    begin
 230       return Complex'(Left, R (Right));
 231    end "+";
 232
 233    ---------
 234    -- "-" --
 235    ---------
 236
 237    function "-" (Right : Complex) return Complex is
 238    begin
 239       return (-Right.Re, -Right.Im);
 240    end "-";
 241
 242    function "-" (Left, Right : Complex) return Complex is
 243    begin
 244       return (Left.Re - Right.Re, Left.Im - Right.Im);
 245    end "-";
 246
 247    function "-" (Right : Imaginary) return Imaginary is
 248    begin
 249       return Imaginary (-R (Right));
 250    end "-";
 251
 252    function "-" (Left, Right : Imaginary) return Imaginary is
 253    begin
 254       return Imaginary (R (Left) - R (Right));
 255    end "-";
 256
 257    function "-" (Left : Complex; Right : Real'Base) return Complex is
 258    begin
 259       return Complex'(Left.Re - Right, Left.Im);
 260    end "-";
 261
 262    function "-" (Left : Real'Base; Right : Complex) return Complex is
 263    begin
 264       return Complex'(Left - Right.Re, -Right.Im);
 265    end "-";
 266
 267    function "-" (Left : Complex; Right : Imaginary) return Complex is
 268    begin
 269       return Complex'(Left.Re, Left.Im - R (Right));
 270    end "-";
 271
 272    function "-" (Left : Imaginary; Right : Complex) return Complex is
 273    begin
 274       return Complex'(-Right.Re, R (Left) - Right.Im);
 275    end "-";
 276
 277    function "-" (Left : Imaginary; Right : Real'Base) return Complex is
 278    begin
 279       return Complex'(-Right, R (Left));
 280    end "-";
 281
 282    function "-" (Left : Real'Base; Right : Imaginary) return Complex is
 283    begin
 284       return Complex'(Left, -R (Right));
 285    end "-";
 286
 287    ---------
 288    -- "/" --
 289    ---------
 290
 291    function "/" (Left, Right : Complex) return Complex is
 292       a : constant R := Left.Re;
 293       b : constant R := Left.Im;
 294       c : constant R := Right.Re;
 295       d : constant R := Right.Im;
 296
 297    begin
 298       if c = 0.0 and then d = 0.0 then
 299          raise Constraint_Error;
 300       else
 301          return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
 302                          Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
 303       end if;
 304    end "/";
 305
 306    function "/" (Left, Right : Imaginary) return Real'Base is
 307    begin
 308       return R (Left) / R (Right);
 309    end "/";
 310
 311    function "/" (Left : Complex; Right : Real'Base) return Complex is
 312    begin
 313       return Complex'(Left.Re / Right, Left.Im / Right);
 314    end "/";
 315
 316    function "/" (Left : Real'Base; Right : Complex) return Complex is
 317       a : constant R := Left;
 318       c : constant R := Right.Re;
 319       d : constant R := Right.Im;
 320    begin
 321       return Complex'(Re =>  (a * c) / (c ** 2 + d ** 2),
 322                       Im => -(a * d) / (c ** 2 + d ** 2));
 323    end "/";
 324
 325    function "/" (Left : Complex; Right : Imaginary) return Complex is
 326       a : constant R := Left.Re;
 327       b : constant R := Left.Im;
 328       d : constant R := R (Right);
 329
 330    begin
 331       return (b / d,  -a / d);
 332    end "/";
 333
 334    function "/" (Left : Imaginary; Right : Complex) return Complex is
 335       b : constant R := R (Left);
 336       c : constant R := Right.Re;
 337       d : constant R := Right.Im;
 338
 339    begin
 340       return (Re => b * d / (c ** 2 + d ** 2),
 341               Im => b * c / (c ** 2 + d ** 2));
 342    end "/";
 343
 344    function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
 345    begin
 346       return Imaginary (R (Left) / Right);
 347    end "/";
 348
 349    function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
 350    begin
 351       return Imaginary (-Left / R (Right));
 352    end "/";
 353
 354    ---------
 355    -- "<" --
 356    ---------
 357
 358    function "<" (Left, Right : Imaginary) return Boolean is
 359    begin
 360       return R (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    -- "abs" --
 392    -----------
 393
 394    function "abs" (Right : Imaginary) return Real'Base is
 395    begin
 396       return abs R (Right);
 397    end "abs";
 398
 399    --------------
 400    -- Argument --
 401    --------------
 402
 403    function Argument (X : Complex) return Real'Base is
 404       a   : constant R := X.Re;
 405       b   : constant R := X.Im;
 406       arg : R;
 407
 408    begin
 409       if b = 0.0 then
 410
 411          if a >= 0.0 then
 412             return 0.0;
 413          else
 414             return R'Copy_Sign (Pi, b);
 415          end if;
 416
 417       elsif a = 0.0 then
 418
 419          if b >= 0.0 then
 420             return Half_Pi;
 421          else
 422             return -Half_Pi;
 423          end if;
 424
 425       else
 426          arg := R (Atan (Double (abs (b / a))));
 427
 428          if a > 0.0 then
 429             if b > 0.0 then
 430                return arg;
 431             else                  --  b < 0.0
 432                return -arg;
 433             end if;
 434
 435          else                     --  a < 0.0
 436             if b >= 0.0 then
 437                return Pi - arg;
 438             else                  --  b < 0.0
 439                return -(Pi - arg);
 440             end if;
 441          end if;
 442       end if;
 443
 444    exception
 445       when Constraint_Error =>
 446          if b > 0.0 then
 447             return Half_Pi;
 448          else
 449             return -Half_Pi;
 450          end if;
 451    end Argument;
 452
 453    function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
 454    begin
 455       if Cycle > 0.0 then
 456          return Argument (X) * Cycle / Two_Pi;
 457       else
 458          raise Argument_Error;
 459       end if;
 460    end Argument;
 461
 462    ----------------------------
 463    -- Compose_From_Cartesian --
 464    ----------------------------
 465
 466    function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
 467    begin
 468       return (Re, Im);
 469    end Compose_From_Cartesian;
 470
 471    function Compose_From_Cartesian (Re : Real'Base) return Complex is
 472    begin
 473       return (Re, 0.0);
 474    end Compose_From_Cartesian;
 475
 476    function Compose_From_Cartesian (Im : Imaginary) return Complex is
 477    begin
 478       return (0.0, R (Im));
 479    end Compose_From_Cartesian;
 480
 481    ------------------------
 482    -- Compose_From_Polar --
 483    ------------------------
 484
 485    function Compose_From_Polar (
 486      Modulus, Argument : Real'Base)
 487      return Complex
 488    is
 489    begin
 490       if Modulus = 0.0 then
 491          return (0.0, 0.0);
 492       else
 493          return (Modulus * R (Cos (Double (Argument))),
 494                  Modulus * R (Sin (Double (Argument))));
 495       end if;
 496    end Compose_From_Polar;
 497
 498    function Compose_From_Polar (
 499      Modulus, Argument, Cycle : Real'Base)
 500      return Complex
 501    is
 502       Arg : Real'Base;
 503
 504    begin
 505       if Modulus = 0.0 then
 506          return (0.0, 0.0);
 507
 508       elsif Cycle > 0.0 then
 509          if Argument = 0.0 then
 510             return (Modulus, 0.0);
 511
 512          elsif Argument = Cycle / 4.0 then
 513             return (0.0, Modulus);
 514
 515          elsif Argument = Cycle / 2.0 then
 516             return (-Modulus, 0.0);
 517
 518          elsif Argument = 3.0 * Cycle / R (4.0) then
 519             return (0.0, -Modulus);
 520          else
 521             Arg := Two_Pi * Argument / Cycle;
 522             return (Modulus * R (Cos (Double (Arg))),
 523                     Modulus * R (Sin (Double (Arg))));
 524          end if;
 525       else
 526          raise Argument_Error;
 527       end if;
 528    end Compose_From_Polar;
 529
 530    ---------------
 531    -- Conjugate --
 532    ---------------
 533
 534    function Conjugate (X : Complex) return Complex is
 535    begin
 536       return Complex'(X.Re, -X.Im);
 537    end Conjugate;
 538
 539    --------
 540    -- Im --
 541    --------
 542
 543    function Im (X : Complex) return Real'Base is
 544    begin
 545       return X.Im;
 546    end Im;
 547
 548    function Im (X : Imaginary) return Real'Base is
 549    begin
 550       return R (X);
 551    end Im;
 552
 553    -------------
 554    -- Modulus --
 555    -------------
 556
 557    function Modulus (X : Complex) return Real'Base is
 558       Re2, Im2 : R;
 559
 560    begin
 561
 562       begin
 563          Re2 := X.Re ** 2;
 564
 565          --  To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
 566          --  compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
 567          --  squaring does not raise constraint_error but generates infinity,
 568          --  we can use an explicit comparison to determine whether to use
 569          --  the scaling expression.
 570
 571          if Re2 > R'Last then
 572             raise Constraint_Error;
 573          end if;
 574
 575       exception
 576          when Constraint_Error =>
 577             return abs (X.Re)
 578               * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
 579       end;
 580
 581       begin
 582          Im2 := X.Im ** 2;
 583
 584          if Im2 > R'Last then
 585             raise Constraint_Error;
 586          end if;
 587
 588       exception
 589          when Constraint_Error =>
 590             return abs (X.Im)
 591               * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
 592       end;
 593
 594       --  Now deal with cases of underflow. If only one of the squares
 595       --  underflows, return the modulus of the other component. If both
 596       --  squares underflow, use scaling as above.
 597
 598       if Re2 = 0.0 then
 599
 600          if X.Re = 0.0 then
 601             return abs (X.Im);
 602
 603          elsif Im2 = 0.0 then
 604
 605             if X.Im = 0.0 then
 606                return abs (X.Re);
 607
 608             else
 609                if abs (X.Re) > abs (X.Im) then
 610                   return
 611                     abs (X.Re)
 612                       * R (Sqrt (Double (R (1.0) + (X.Im / X.Re) ** 2)));
 613                else
 614                   return
 615                     abs (X.Im)
 616                       * R (Sqrt (Double (R (1.0) + (X.Re / X.Im) ** 2)));
 617                end if;
 618             end if;
 619
 620          else
 621             return abs (X.Im);
 622          end if;
 623
 624
 625       elsif Im2 = 0.0 then
 626          return abs (X.Re);
 627
 628          --  in all other cases, the naive computation will do.
 629
 630       else
 631          return R (Sqrt (Double (Re2 + Im2)));
 632       end if;
 633    end Modulus;
 634
 635    --------
 636    -- Re --
 637    --------
 638
 639    function Re (X : Complex) return Real'Base is
 640    begin
 641       return X.Re;
 642    end Re;
 643
 644    ------------
 645    -- Set_Im --
 646    ------------
 647
 648    procedure Set_Im (X : in out Complex; Im : in Real'Base) is
 649    begin
 650       X.Im := Im;
 651    end Set_Im;
 652
 653    procedure Set_Im (X : out Imaginary; Im : in Real'Base) is
 654    begin
 655       X := Imaginary (Im);
 656    end Set_Im;
 657
 658    ------------
 659    -- Set_Re --
 660    ------------
 661
 662    procedure Set_Re (X : in out Complex; Re : in Real'Base) is
 663    begin
 664       X.Re := Re;
 665    end Set_Re;
 666
 667 end Ada.Numerics.Generic_Complex_Types;