gcc/ada/a-ngcoty.adb

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