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