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