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-2006, 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 2,  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.  See the GNU General Public License --
  17 -- for  more details.  You should have  received  a copy of the GNU General --
  18 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19 -- to  the  Free Software Foundation,  51  Franklin  Street,  Fifth  Floor, --
  20 -- Boston, MA 02110-1301, USA.                                              --
  21 --                                                                          --
  22 -- As a special exception,  if other files  instantiate  generics from this --
  23 -- unit, or you link  this unit with other files  to produce an executable, --
  24 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  25 -- covered  by the  GNU  General  Public  License.  This exception does not --
  26 -- however invalidate  any other reasons why  the executable file  might be --
  27 -- covered by the  GNU Public License.                                      --
  28 --                                                                          --
  29 -- GNAT was originally developed  by the GNAT team at  New York University. --
  30 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  31 --                                                                          --
  32 ------------------------------------------------------------------------------
  33
  34 with Ada.Numerics.Aux; use Ada.Numerics.Aux;
  35
  36 package body Ada.Numerics.Generic_Complex_Types is
  37
  38    subtype R is Real'Base;
  39
  40    Two_Pi  : constant R := R (2.0) * Pi;
  41    Half_Pi : constant R := Pi / R (2.0);
  42
  43    ---------
  44    -- "*" --
  45    ---------
  46
  47    function "*" (Left, Right : Complex) return Complex is
  48       X : R;
  49       Y : R;
  50
  51    begin
  52       X := Left.Re * Right.Re - Left.Im * Right.Im;
  53       Y := Left.Re * Right.Im + Left.Im * Right.Re;
  54
  55       --  If either component overflows, try to scale
  56
  57       if abs (X) > R'Last then
  58          X := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Re / 2.0)
  59                 - R'(Left.Im / 2.0) * R'(Right.Im / 2.0));
  60       end if;
  61
  62       if abs (Y) > R'Last then
  63          Y := R'(4.0) * (R'(Left.Re / 2.0)  * R'(Right.Im / 2.0)
  64                 - R'(Left.Im / 2.0) * R'(Right.Re / 2.0));
  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
 147             while Exp /= 0 loop
 148                if Exp rem 2 /= 0 then
 149                   Result := Result * Factor;
 150                end if;
 151
 152                Factor := Factor * Factor;
 153                Exp := Exp / 2;
 154             end loop;
 155
 156             return R'(1.0) / Result;
 157
 158          exception
 159
 160             when Constraint_Error =>
 161                return (0.0, 0.0);
 162          end;
 163       end if;
 164    end "**";
 165
 166    function "**" (Left : Imaginary; Right : Integer) return Complex is
 167       M : constant R := R (Left) ** Right;
 168    begin
 169       case Right mod 4 is
 170          when 0 => return (M,   0.0);
 171          when 1 => return (0.0, M);
 172          when 2 => return (-M,  0.0);
 173          when 3 => return (0.0, -M);
 174          when others => raise Program_Error;
 175       end case;
 176    end "**";
 177
 178    ---------
 179    -- "+" --
 180    ---------
 181
 182    function "+" (Right : Complex) return Complex is
 183    begin
 184       return Right;
 185    end "+";
 186
 187    function "+" (Left, Right : Complex) return Complex is
 188    begin
 189       return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
 190    end "+";
 191
 192    function "+" (Right : Imaginary) return Imaginary is
 193    begin
 194       return Right;
 195    end "+";
 196
 197    function "+" (Left, Right : Imaginary) return Imaginary is
 198    begin
 199       return Imaginary (R (Left) + R (Right));
 200    end "+";
 201
 202    function "+" (Left : Complex; Right : Real'Base) return Complex is
 203    begin
 204       return Complex'(Left.Re + Right, Left.Im);
 205    end "+";
 206
 207    function "+" (Left : Real'Base; Right : Complex) return Complex is
 208    begin
 209       return Complex'(Left + Right.Re, Right.Im);
 210    end "+";
 211
 212    function "+" (Left : Complex; Right : Imaginary) return Complex is
 213    begin
 214       return Complex'(Left.Re, Left.Im + R (Right));
 215    end "+";
 216
 217    function "+" (Left : Imaginary; Right : Complex) return Complex is
 218    begin
 219       return Complex'(Right.Re, R (Left) + Right.Im);
 220    end "+";
 221
 222    function "+" (Left : Imaginary; Right : Real'Base) return Complex is
 223    begin
 224       return Complex'(Right, R (Left));
 225    end "+";
 226
 227    function "+" (Left : Real'Base; Right : Imaginary) return Complex is
 228    begin
 229       return Complex'(Left, R (Right));
 230    end "+";
 231
 232    ---------
 233    -- "-" --
 234    ---------
 235
 236    function "-" (Right : Complex) return Complex is
 237    begin
 238       return (-Right.Re, -Right.Im);
 239    end "-";
 240
 241    function "-" (Left, Right : Complex) return Complex is
 242    begin
 243       return (Left.Re - Right.Re, Left.Im - Right.Im);
 244    end "-";
 245
 246    function "-" (Right : Imaginary) return Imaginary is
 247    begin
 248       return Imaginary (-R (Right));
 249    end "-";
 250
 251    function "-" (Left, Right : Imaginary) return Imaginary is
 252    begin
 253       return Imaginary (R (Left) - R (Right));
 254    end "-";
 255
 256    function "-" (Left : Complex; Right : Real'Base) return Complex is
 257    begin
 258       return Complex'(Left.Re - Right, Left.Im);
 259    end "-";
 260
 261    function "-" (Left : Real'Base; Right : Complex) return Complex is
 262    begin
 263       return Complex'(Left - Right.Re, -Right.Im);
 264    end "-";
 265
 266    function "-" (Left : Complex; Right : Imaginary) return Complex is
 267    begin
 268       return Complex'(Left.Re, Left.Im - R (Right));
 269    end "-";
 270
 271    function "-" (Left : Imaginary; Right : Complex) return Complex is
 272    begin
 273       return Complex'(-Right.Re, R (Left) - Right.Im);
 274    end "-";
 275
 276    function "-" (Left : Imaginary; Right : Real'Base) return Complex is
 277    begin
 278       return Complex'(-Right, R (Left));
 279    end "-";
 280
 281    function "-" (Left : Real'Base; Right : Imaginary) return Complex is
 282    begin
 283       return Complex'(Left, -R (Right));
 284    end "-";
 285
 286    ---------
 287    -- "/" --
 288    ---------
 289
 290    function "/" (Left, Right : Complex) return Complex is
 291       a : constant R := Left.Re;
 292       b : constant R := Left.Im;
 293       c : constant R := Right.Re;
 294       d : constant R := Right.Im;
 295
 296    begin
 297       if c = 0.0 and then d = 0.0 then
 298          raise Constraint_Error;
 299       else
 300          return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
 301                          Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
 302       end if;
 303    end "/";
 304
 305    function "/" (Left, Right : Imaginary) return Real'Base is
 306    begin
 307       return R (Left) / R (Right);
 308    end "/";
 309
 310    function "/" (Left : Complex; Right : Real'Base) return Complex is
 311    begin
 312       return Complex'(Left.Re / Right, Left.Im / Right);
 313    end "/";
 314
 315    function "/" (Left : Real'Base; Right : Complex) return Complex is
 316       a : constant R := Left;
 317       c : constant R := Right.Re;
 318       d : constant R := Right.Im;
 319    begin
 320       return Complex'(Re =>   (a * c) / (c ** 2 + d ** 2),
 321                       Im => -((a * d) / (c ** 2 + d ** 2)));
 322    end "/";
 323
 324    function "/" (Left : Complex; Right : Imaginary) return Complex is
 325       a : constant R := Left.Re;
 326       b : constant R := Left.Im;
 327       d : constant R := R (Right);
 328
 329    begin
 330       return (b / d,  -(a / d));
 331    end "/";
 332
 333    function "/" (Left : Imaginary; Right : Complex) return Complex is
 334       b : constant R := R (Left);
 335       c : constant R := Right.Re;
 336       d : constant R := Right.Im;
 337
 338    begin
 339       return (Re => b * d / (c ** 2 + d ** 2),
 340               Im => b * c / (c ** 2 + d ** 2));
 341    end "/";
 342
 343    function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
 344    begin
 345       return Imaginary (R (Left) / Right);
 346    end "/";
 347
 348    function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
 349    begin
 350       return Imaginary (-(Left / R (Right)));
 351    end "/";
 352
 353    ---------
 354    -- "<" --
 355    ---------
 356
 357    function "<" (Left, Right : Imaginary) return Boolean is
 358    begin
 359       return R (Left) < R (Right);
 360    end "<";
 361
 362    ----------
 363    -- "<=" --
 364    ----------
 365
 366    function "<=" (Left, Right : Imaginary) return Boolean is
 367    begin
 368       return R (Left) <= R (Right);
 369    end "<=";
 370
 371    ---------
 372    -- ">" --
 373    ---------
 374
 375    function ">" (Left, Right : Imaginary) return Boolean is
 376    begin
 377       return R (Left) > R (Right);
 378    end ">";
 379
 380    ----------
 381    -- ">=" --
 382    ----------
 383
 384    function ">=" (Left, Right : Imaginary) return Boolean is
 385    begin
 386       return R (Left) >= R (Right);
 387    end ">=";
 388
 389    -----------
 390    -- "abs" --
 391    -----------
 392
 393    function "abs" (Right : Imaginary) return Real'Base is
 394    begin
 395       return abs R (Right);
 396    end "abs";
 397
 398    --------------
 399    -- Argument --
 400    --------------
 401
 402    function Argument (X : Complex) return Real'Base is
 403       a   : constant R := X.Re;
 404       b   : constant R := X.Im;
 405       arg : R;
 406
 407    begin
 408       if b = 0.0 then
 409
 410          if a >= 0.0 then
 411             return 0.0;
 412          else
 413             return R'Copy_Sign (Pi, b);
 414          end if;
 415
 416       elsif a = 0.0 then
 417
 418          if b >= 0.0 then
 419             return Half_Pi;
 420          else
 421             return -Half_Pi;
 422          end if;
 423
 424       else
 425          arg := R (Atan (Double (abs (b / a))));
 426
 427          if a > 0.0 then
 428             if b > 0.0 then
 429                return arg;
 430             else                  --  b < 0.0
 431                return -arg;
 432             end if;
 433
 434          else                     --  a < 0.0
 435             if b >= 0.0 then
 436                return Pi - arg;
 437             else                  --  b < 0.0
 438                return -(Pi - arg);
 439             end if;
 440          end if;
 441       end if;
 442
 443    exception
 444       when Constraint_Error =>
 445          if b > 0.0 then
 446             return Half_Pi;
 447          else
 448             return -Half_Pi;
 449          end if;
 450    end Argument;
 451
 452    function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
 453    begin
 454       if Cycle > 0.0 then
 455          return Argument (X) * Cycle / Two_Pi;
 456       else
 457          raise Argument_Error;
 458       end if;
 459    end Argument;
 460
 461    ----------------------------
 462    -- Compose_From_Cartesian --
 463    ----------------------------
 464
 465    function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
 466    begin
 467       return (Re, Im);
 468    end Compose_From_Cartesian;
 469
 470    function Compose_From_Cartesian (Re : Real'Base) return Complex is
 471    begin
 472       return (Re, 0.0);
 473    end Compose_From_Cartesian;
 474
 475    function Compose_From_Cartesian (Im : Imaginary) return Complex is
 476    begin
 477       return (0.0, R (Im));
 478    end Compose_From_Cartesian;
 479
 480    ------------------------
 481    -- Compose_From_Polar --
 482    ------------------------
 483
 484    function Compose_From_Polar (
 485      Modulus, Argument : Real'Base)
 486      return Complex
 487    is
 488    begin
 489       if Modulus = 0.0 then
 490          return (0.0, 0.0);
 491       else
 492          return (Modulus * R (Cos (Double (Argument))),
 493                  Modulus * R (Sin (Double (Argument))));
 494       end if;
 495    end Compose_From_Polar;
 496
 497    function Compose_From_Polar (
 498      Modulus, Argument, Cycle : Real'Base)
 499      return Complex
 500    is
 501       Arg : Real'Base;
 502
 503    begin
 504       if Modulus = 0.0 then
 505          return (0.0, 0.0);
 506
 507       elsif Cycle > 0.0 then
 508          if Argument = 0.0 then
 509             return (Modulus, 0.0);
 510
 511          elsif Argument = Cycle / 4.0 then
 512             return (0.0, Modulus);
 513
 514          elsif Argument = Cycle / 2.0 then
 515             return (-Modulus, 0.0);
 516
 517          elsif Argument = 3.0 * Cycle / R (4.0) then
 518             return (0.0, -Modulus);
 519          else
 520             Arg := Two_Pi * Argument / Cycle;
 521             return (Modulus * R (Cos (Double (Arg))),
 522                     Modulus * R (Sin (Double (Arg))));
 523          end if;
 524       else
 525          raise Argument_Error;
 526       end if;
 527    end Compose_From_Polar;
 528
 529    ---------------
 530    -- Conjugate --
 531    ---------------
 532
 533    function Conjugate (X : Complex) return Complex is
 534    begin
 535       return Complex'(X.Re, -X.Im);
 536    end Conjugate;
 537
 538    --------
 539    -- Im --
 540    --------
 541
 542    function Im (X : Complex) return Real'Base is
 543    begin
 544       return X.Im;
 545    end Im;
 546
 547    function Im (X : Imaginary) return Real'Base is
 548    begin
 549       return R (X);
 550    end Im;
 551
 552    -------------
 553    -- Modulus --
 554    -------------
 555
 556    function Modulus (X : Complex) return Real'Base is
 557       Re2, Im2 : R;
 558
 559    begin
 560
 561       begin
 562          Re2 := X.Re ** 2;
 563
 564          --  To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
 565          --  compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
 566          --  squaring does not raise constraint_error but generates infinity,
 567          --  we can use an explicit comparison to determine whether to use
 568          --  the scaling expression.
 569
 570          --  The scaling expression is computed in double format throughout
 571          --  in order to prevent inaccuracies on machines where not all
 572          --  immediate expressions are rounded, such as PowerPC.
 573
 574          if Re2 > R'Last then
 575             raise Constraint_Error;
 576          end if;
 577
 578       exception
 579          when Constraint_Error =>
 580             return R (Double (abs (X.Re))
 581               * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
 582       end;
 583
 584       begin
 585          Im2 := X.Im ** 2;
 586
 587          if Im2 > R'Last then
 588             raise Constraint_Error;
 589          end if;
 590
 591       exception
 592          when Constraint_Error =>
 593             return R (Double (abs (X.Im))
 594               * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
 595       end;
 596
 597       --  Now deal with cases of underflow. If only one of the squares
 598       --  underflows, return the modulus of the other component. If both
 599       --  squares underflow, use scaling as above.
 600
 601       if Re2 = 0.0 then
 602
 603          if X.Re = 0.0 then
 604             return abs (X.Im);
 605
 606          elsif Im2 = 0.0 then
 607
 608             if X.Im = 0.0 then
 609                return abs (X.Re);
 610
 611             else
 612                if abs (X.Re) > abs (X.Im) then
 613                   return
 614                     R (Double (abs (X.Re))
 615                       * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
 616                else
 617                   return
 618                     R (Double (abs (X.Im))
 619                       * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
 620                end if;
 621             end if;
 622
 623          else
 624             return abs (X.Im);
 625          end if;
 626
 627       elsif Im2 = 0.0 then
 628          return abs (X.Re);
 629
 630       --  In all other cases, the naive computation will do
 631
 632       else
 633          return R (Sqrt (Double (Re2 + Im2)));
 634       end if;
 635    end Modulus;
 636
 637    --------
 638    -- Re --
 639    --------
 640
 641    function Re (X : Complex) return Real'Base is
 642    begin
 643       return X.Re;
 644    end Re;
 645
 646    ------------
 647    -- Set_Im --
 648    ------------
 649
 650    procedure Set_Im (X : in out Complex; Im : Real'Base) is
 651    begin
 652       X.Im := Im;
 653    end Set_Im;
 654
 655    procedure Set_Im (X : out Imaginary; Im : Real'Base) is
 656    begin
 657       X := Imaginary (Im);
 658    end Set_Im;
 659
 660    ------------
 661    -- Set_Re --
 662    ------------
 663
 664    procedure Set_Re (X : in out Complex; Re : Real'Base) is
 665    begin
 666       X.Re := Re;
 667    end Set_Re;
 668
 669 end Ada.Numerics.Generic_Complex_Types;