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-2010, 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
  47       Scale : constant R := R (R'Machine_Radix) ** ((R'Machine_Emax - 1) / 2);
  48       --  In case of overflow, scale the operands by the largest power of the
  49       --  radix (to avoid rounding error), so that the square of the scale does
  50       --  not overflow itself.
  51
  52       X : R;
  53       Y : R;
  54
  55    begin
  56       X := Left.Re * Right.Re - Left.Im * Right.Im;
  57       Y := Left.Re * Right.Im + Left.Im * Right.Re;
  58
  59       --  If either component overflows, try to scale (skip in fast math mode)
  60
  61       if not Standard'Fast_Math then
  62
  63          --  Note that the test below is written as a negation. This is to
  64          --  account for the fact that X and Y may be NaNs, because both of
  65          --  their operands could overflow. Given that all operations on NaNs
  66          --  return false, the test can only be written thus.
  67
  68          if not (abs (X) <= R'Last) then
  69             X := Scale**2 * ((Left.Re / Scale) * (Right.Re / Scale) -
  70                              (Left.Im / Scale) * (Right.Im / Scale));
  71          end if;
  72
  73          if not (abs (Y) <= R'Last) then
  74             Y := Scale**2 * ((Left.Re / Scale) * (Right.Im / Scale)
  75                            + (Left.Im / Scale) * (Right.Re / Scale));
  76          end if;
  77       end if;
  78
  79       return (X, Y);
  80    end "*";
  81
  82    function "*" (Left, Right : Imaginary) return Real'Base is
  83    begin
  84       return -(R (Left) * R (Right));
  85    end "*";
  86
  87    function "*" (Left : Complex; Right : Real'Base) return Complex is
  88    begin
  89       return Complex'(Left.Re * Right, Left.Im * Right);
  90    end "*";
  91
  92    function "*" (Left : Real'Base; Right : Complex) return Complex is
  93    begin
  94       return (Left * Right.Re, Left * Right.Im);
  95    end "*";
  96
  97    function "*" (Left : Complex; Right : Imaginary) return Complex is
  98    begin
  99       return Complex'(-(Left.Im * R (Right)), Left.Re * R (Right));
 100    end "*";
 101
 102    function "*" (Left : Imaginary; Right : Complex) return Complex is
 103    begin
 104       return Complex'(-(R (Left) * Right.Im), R (Left) * Right.Re);
 105    end "*";
 106
 107    function "*" (Left : Imaginary; Right : Real'Base) return Imaginary is
 108    begin
 109       return Left * Imaginary (Right);
 110    end "*";
 111
 112    function "*" (Left : Real'Base; Right : Imaginary) return Imaginary is
 113    begin
 114       return Imaginary (Left * R (Right));
 115    end "*";
 116
 117    ----------
 118    -- "**" --
 119    ----------
 120
 121    function "**" (Left : Complex; Right : Integer) return Complex is
 122       Result : Complex := (1.0, 0.0);
 123       Factor : Complex := Left;
 124       Exp    : Integer := Right;
 125
 126    begin
 127       --  We use the standard logarithmic approach, Exp gets shifted right
 128       --  testing successive low order bits and Factor is the value of the
 129       --  base raised to the next power of 2. For positive exponents we
 130       --  multiply the result by this factor, for negative exponents, we
 131       --  divide by this factor.
 132
 133       if Exp >= 0 then
 134
 135          --  For a positive exponent, if we get a constraint error during
 136          --  this loop, it is an overflow, and the constraint error will
 137          --  simply be passed on to the caller.
 138
 139          while Exp /= 0 loop
 140             if Exp rem 2 /= 0 then
 141                Result := Result * Factor;
 142             end if;
 143
 144             Factor := Factor * Factor;
 145             Exp := Exp / 2;
 146          end loop;
 147
 148          return Result;
 149
 150       else -- Exp < 0 then
 151
 152          --  For the negative exponent case, a constraint error during this
 153          --  calculation happens if Factor gets too large, and the proper
 154          --  response is to return 0.0, since what we essentially have is
 155          --  1.0 / infinity, and the closest model number will be zero.
 156
 157          begin
 158             while Exp /= 0 loop
 159                if Exp rem 2 /= 0 then
 160                   Result := Result * Factor;
 161                end if;
 162
 163                Factor := Factor * Factor;
 164                Exp := Exp / 2;
 165             end loop;
 166
 167             return R'(1.0) / Result;
 168
 169          exception
 170             when Constraint_Error =>
 171                return (0.0, 0.0);
 172          end;
 173       end if;
 174    end "**";
 175
 176    function "**" (Left : Imaginary; Right : Integer) return Complex is
 177       M : constant R := R (Left) ** Right;
 178    begin
 179       case Right mod 4 is
 180          when 0 => return (M,   0.0);
 181          when 1 => return (0.0, M);
 182          when 2 => return (-M,  0.0);
 183          when 3 => return (0.0, -M);
 184          when others => raise Program_Error;
 185       end case;
 186    end "**";
 187
 188    ---------
 189    -- "+" --
 190    ---------
 191
 192    function "+" (Right : Complex) return Complex is
 193    begin
 194       return Right;
 195    end "+";
 196
 197    function "+" (Left, Right : Complex) return Complex is
 198    begin
 199       return Complex'(Left.Re + Right.Re, Left.Im + Right.Im);
 200    end "+";
 201
 202    function "+" (Right : Imaginary) return Imaginary is
 203    begin
 204       return Right;
 205    end "+";
 206
 207    function "+" (Left, Right : Imaginary) return Imaginary is
 208    begin
 209       return Imaginary (R (Left) + R (Right));
 210    end "+";
 211
 212    function "+" (Left : Complex; Right : Real'Base) return Complex is
 213    begin
 214       return Complex'(Left.Re + Right, Left.Im);
 215    end "+";
 216
 217    function "+" (Left : Real'Base; Right : Complex) return Complex is
 218    begin
 219       return Complex'(Left + Right.Re, Right.Im);
 220    end "+";
 221
 222    function "+" (Left : Complex; Right : Imaginary) return Complex is
 223    begin
 224       return Complex'(Left.Re, Left.Im + R (Right));
 225    end "+";
 226
 227    function "+" (Left : Imaginary; Right : Complex) return Complex is
 228    begin
 229       return Complex'(Right.Re, R (Left) + Right.Im);
 230    end "+";
 231
 232    function "+" (Left : Imaginary; Right : Real'Base) return Complex is
 233    begin
 234       return Complex'(Right, R (Left));
 235    end "+";
 236
 237    function "+" (Left : Real'Base; Right : Imaginary) return Complex is
 238    begin
 239       return Complex'(Left, R (Right));
 240    end "+";
 241
 242    ---------
 243    -- "-" --
 244    ---------
 245
 246    function "-" (Right : Complex) return Complex is
 247    begin
 248       return (-Right.Re, -Right.Im);
 249    end "-";
 250
 251    function "-" (Left, Right : Complex) return Complex is
 252    begin
 253       return (Left.Re - Right.Re, Left.Im - Right.Im);
 254    end "-";
 255
 256    function "-" (Right : Imaginary) return Imaginary is
 257    begin
 258       return Imaginary (-R (Right));
 259    end "-";
 260
 261    function "-" (Left, Right : Imaginary) return Imaginary is
 262    begin
 263       return Imaginary (R (Left) - R (Right));
 264    end "-";
 265
 266    function "-" (Left : Complex; Right : Real'Base) return Complex is
 267    begin
 268       return Complex'(Left.Re - Right, Left.Im);
 269    end "-";
 270
 271    function "-" (Left : Real'Base; Right : Complex) return Complex is
 272    begin
 273       return Complex'(Left - Right.Re, -Right.Im);
 274    end "-";
 275
 276    function "-" (Left : Complex; Right : Imaginary) return Complex is
 277    begin
 278       return Complex'(Left.Re, Left.Im - R (Right));
 279    end "-";
 280
 281    function "-" (Left : Imaginary; Right : Complex) return Complex is
 282    begin
 283       return Complex'(-Right.Re, R (Left) - Right.Im);
 284    end "-";
 285
 286    function "-" (Left : Imaginary; Right : Real'Base) return Complex is
 287    begin
 288       return Complex'(-Right, R (Left));
 289    end "-";
 290
 291    function "-" (Left : Real'Base; Right : Imaginary) return Complex is
 292    begin
 293       return Complex'(Left, -R (Right));
 294    end "-";
 295
 296    ---------
 297    -- "/" --
 298    ---------
 299
 300    function "/" (Left, Right : Complex) return Complex is
 301       a : constant R := Left.Re;
 302       b : constant R := Left.Im;
 303       c : constant R := Right.Re;
 304       d : constant R := Right.Im;
 305
 306    begin
 307       if c = 0.0 and then d = 0.0 then
 308          raise Constraint_Error;
 309       else
 310          return Complex'(Re => ((a * c) + (b * d)) / (c ** 2 + d ** 2),
 311                          Im => ((b * c) - (a * d)) / (c ** 2 + d ** 2));
 312       end if;
 313    end "/";
 314
 315    function "/" (Left, Right : Imaginary) return Real'Base is
 316    begin
 317       return R (Left) / R (Right);
 318    end "/";
 319
 320    function "/" (Left : Complex; Right : Real'Base) return Complex is
 321    begin
 322       return Complex'(Left.Re / Right, Left.Im / Right);
 323    end "/";
 324
 325    function "/" (Left : Real'Base; Right : Complex) return Complex is
 326       a : constant R := Left;
 327       c : constant R := Right.Re;
 328       d : constant R := Right.Im;
 329    begin
 330       return Complex'(Re =>   (a * c) / (c ** 2 + d ** 2),
 331                       Im => -((a * d) / (c ** 2 + d ** 2)));
 332    end "/";
 333
 334    function "/" (Left : Complex; Right : Imaginary) return Complex is
 335       a : constant R := Left.Re;
 336       b : constant R := Left.Im;
 337       d : constant R := R (Right);
 338
 339    begin
 340       return (b / d,  -(a / d));
 341    end "/";
 342
 343    function "/" (Left : Imaginary; Right : Complex) return Complex is
 344       b : constant R := R (Left);
 345       c : constant R := Right.Re;
 346       d : constant R := Right.Im;
 347
 348    begin
 349       return (Re => b * d / (c ** 2 + d ** 2),
 350               Im => b * c / (c ** 2 + d ** 2));
 351    end "/";
 352
 353    function "/" (Left : Imaginary; Right : Real'Base) return Imaginary is
 354    begin
 355       return Imaginary (R (Left) / Right);
 356    end "/";
 357
 358    function "/" (Left : Real'Base; Right : Imaginary) return Imaginary is
 359    begin
 360       return Imaginary (-(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    -- ">=" --
 392    ----------
 393
 394    function ">=" (Left, Right : Imaginary) return Boolean is
 395    begin
 396       return R (Left) >= R (Right);
 397    end ">=";
 398
 399    -----------
 400    -- "abs" --
 401    -----------
 402
 403    function "abs" (Right : Imaginary) return Real'Base is
 404    begin
 405       return abs R (Right);
 406    end "abs";
 407
 408    --------------
 409    -- Argument --
 410    --------------
 411
 412    function Argument (X : Complex) return Real'Base is
 413       a   : constant R := X.Re;
 414       b   : constant R := X.Im;
 415       arg : R;
 416
 417    begin
 418       if b = 0.0 then
 419
 420          if a >= 0.0 then
 421             return 0.0;
 422          else
 423             return R'Copy_Sign (Pi, b);
 424          end if;
 425
 426       elsif a = 0.0 then
 427
 428          if b >= 0.0 then
 429             return Half_Pi;
 430          else
 431             return -Half_Pi;
 432          end if;
 433
 434       else
 435          arg := R (Atan (Double (abs (b / a))));
 436
 437          if a > 0.0 then
 438             if b > 0.0 then
 439                return arg;
 440             else                  --  b < 0.0
 441                return -arg;
 442             end if;
 443
 444          else                     --  a < 0.0
 445             if b >= 0.0 then
 446                return Pi - arg;
 447             else                  --  b < 0.0
 448                return -(Pi - arg);
 449             end if;
 450          end if;
 451       end if;
 452
 453    exception
 454       when Constraint_Error =>
 455          if b > 0.0 then
 456             return Half_Pi;
 457          else
 458             return -Half_Pi;
 459          end if;
 460    end Argument;
 461
 462    function Argument (X : Complex; Cycle : Real'Base) return Real'Base is
 463    begin
 464       if Cycle > 0.0 then
 465          return Argument (X) * Cycle / Two_Pi;
 466       else
 467          raise Argument_Error;
 468       end if;
 469    end Argument;
 470
 471    ----------------------------
 472    -- Compose_From_Cartesian --
 473    ----------------------------
 474
 475    function Compose_From_Cartesian (Re, Im : Real'Base) return Complex is
 476    begin
 477       return (Re, Im);
 478    end Compose_From_Cartesian;
 479
 480    function Compose_From_Cartesian (Re : Real'Base) return Complex is
 481    begin
 482       return (Re, 0.0);
 483    end Compose_From_Cartesian;
 484
 485    function Compose_From_Cartesian (Im : Imaginary) return Complex is
 486    begin
 487       return (0.0, R (Im));
 488    end Compose_From_Cartesian;
 489
 490    ------------------------
 491    -- Compose_From_Polar --
 492    ------------------------
 493
 494    function Compose_From_Polar (
 495      Modulus, Argument : Real'Base)
 496      return Complex
 497    is
 498    begin
 499       if Modulus = 0.0 then
 500          return (0.0, 0.0);
 501       else
 502          return (Modulus * R (Cos (Double (Argument))),
 503                  Modulus * R (Sin (Double (Argument))));
 504       end if;
 505    end Compose_From_Polar;
 506
 507    function Compose_From_Polar (
 508      Modulus, Argument, Cycle : Real'Base)
 509      return Complex
 510    is
 511       Arg : Real'Base;
 512
 513    begin
 514       if Modulus = 0.0 then
 515          return (0.0, 0.0);
 516
 517       elsif Cycle > 0.0 then
 518          if Argument = 0.0 then
 519             return (Modulus, 0.0);
 520
 521          elsif Argument = Cycle / 4.0 then
 522             return (0.0, Modulus);
 523
 524          elsif Argument = Cycle / 2.0 then
 525             return (-Modulus, 0.0);
 526
 527          elsif Argument = 3.0 * Cycle / R (4.0) then
 528             return (0.0, -Modulus);
 529          else
 530             Arg := Two_Pi * Argument / Cycle;
 531             return (Modulus * R (Cos (Double (Arg))),
 532                     Modulus * R (Sin (Double (Arg))));
 533          end if;
 534       else
 535          raise Argument_Error;
 536       end if;
 537    end Compose_From_Polar;
 538
 539    ---------------
 540    -- Conjugate --
 541    ---------------
 542
 543    function Conjugate (X : Complex) return Complex is
 544    begin
 545       return Complex'(X.Re, -X.Im);
 546    end Conjugate;
 547
 548    --------
 549    -- Im --
 550    --------
 551
 552    function Im (X : Complex) return Real'Base is
 553    begin
 554       return X.Im;
 555    end Im;
 556
 557    function Im (X : Imaginary) return Real'Base is
 558    begin
 559       return R (X);
 560    end Im;
 561
 562    -------------
 563    -- Modulus --
 564    -------------
 565
 566    function Modulus (X : Complex) return Real'Base is
 567       Re2, Im2 : R;
 568
 569    begin
 570
 571       begin
 572          Re2 := X.Re ** 2;
 573
 574          --  To compute (a**2 + b**2) ** (0.5) when a**2 may be out of bounds,
 575          --  compute a * (1 + (b/a) **2) ** (0.5). On a machine where the
 576          --  squaring does not raise constraint_error but generates infinity,
 577          --  we can use an explicit comparison to determine whether to use
 578          --  the scaling expression.
 579
 580          --  The scaling expression is computed in double format throughout
 581          --  in order to prevent inaccuracies on machines where not all
 582          --  immediate expressions are rounded, such as PowerPC.
 583
 584          --  ??? same weird test, why not Re2 > R'Last ???
 585          if not (Re2 <= R'Last) then
 586             raise Constraint_Error;
 587          end if;
 588
 589       exception
 590          when Constraint_Error =>
 591             return R (Double (abs (X.Re))
 592               * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
 593       end;
 594
 595       begin
 596          Im2 := X.Im ** 2;
 597
 598          --  ??? same weird test
 599          if not (Im2 <= R'Last) then
 600             raise Constraint_Error;
 601          end if;
 602
 603       exception
 604          when Constraint_Error =>
 605             return R (Double (abs (X.Im))
 606               * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
 607       end;
 608
 609       --  Now deal with cases of underflow. If only one of the squares
 610       --  underflows, return the modulus of the other component. If both
 611       --  squares underflow, use scaling as above.
 612
 613       if Re2 = 0.0 then
 614
 615          if X.Re = 0.0 then
 616             return abs (X.Im);
 617
 618          elsif Im2 = 0.0 then
 619
 620             if X.Im = 0.0 then
 621                return abs (X.Re);
 622
 623             else
 624                if abs (X.Re) > abs (X.Im) then
 625                   return
 626                     R (Double (abs (X.Re))
 627                       * Sqrt (1.0 + (Double (X.Im) / Double (X.Re)) ** 2));
 628                else
 629                   return
 630                     R (Double (abs (X.Im))
 631                       * Sqrt (1.0 + (Double (X.Re) / Double (X.Im)) ** 2));
 632                end if;
 633             end if;
 634
 635          else
 636             return abs (X.Im);
 637          end if;
 638
 639       elsif Im2 = 0.0 then
 640          return abs (X.Re);
 641
 642       --  In all other cases, the naive computation will do
 643
 644       else
 645          return R (Sqrt (Double (Re2 + Im2)));
 646       end if;
 647    end Modulus;
 648
 649    --------
 650    -- Re --
 651    --------
 652
 653    function Re (X : Complex) return Real'Base is
 654    begin
 655       return X.Re;
 656    end Re;
 657
 658    ------------
 659    -- Set_Im --
 660    ------------
 661
 662    procedure Set_Im (X : in out Complex; Im : Real'Base) is
 663    begin
 664       X.Im := Im;
 665    end Set_Im;
 666
 667    procedure Set_Im (X : out Imaginary; Im : Real'Base) is
 668    begin
 669       X := Imaginary (Im);
 670    end Set_Im;
 671
 672    ------------
 673    -- Set_Re --
 674    ------------
 675
 676    procedure Set_Re (X : in out Complex; Re : Real'Base) is
 677    begin
 678       X.Re := Re;
 679    end Set_Re;
 680
 681 end Ada.Numerics.Generic_Complex_Types;