gcc/ada/libgnat/a-ngcefu.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUN-TIME COMPONENTS                         --
   4 --                                                                          --
   5 --            ADA.NUMERICS.GENERIC_COMPLEX_ELEMENTARY_FUNCTIONS             --
   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.Generic_Elementary_Functions;
  33
  34 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
  35
  36    package Elementary_Functions is new
  37       Ada.Numerics.Generic_Elementary_Functions (Real'Base);
  38    use Elementary_Functions;
  39
  40    PI       : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
  41    PI_2     : constant := PI / 2.0;
  42    Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
  43    Log_Two  : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  44
  45    subtype T is Real'Base;
  46
  47    Epsilon                 : constant T := 2.0      ** (1 - T'Model_Mantissa);
  48    Square_Root_Epsilon     : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
  49    Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
  50    Root_Root_Epsilon       : constant T := Sqrt_Two **
  51                                                  ((1 - T'Model_Mantissa) / 2);
  52    Log_Inverse_Epsilon_2   : constant T := T (T'Model_Mantissa - 1) / 2.0;
  53
  54    Complex_Zero : constant Complex := (0.0,  0.0);
  55    Complex_One  : constant Complex := (1.0,  0.0);
  56    Complex_I    : constant Complex := (0.0,  1.0);
  57    Half_Pi      : constant Complex := (PI_2, 0.0);
  58
  59    --------
  60    -- ** --
  61    --------
  62
  63    function "**" (Left : Complex; Right : Complex) return Complex is
  64    begin
  65       if Re (Right) = 0.0
  66         and then Im (Right) = 0.0
  67         and then Re (Left)  = 0.0
  68         and then Im (Left)  = 0.0
  69       then
  70          raise Argument_Error;
  71
  72       elsif Re (Left) = 0.0
  73         and then Im (Left) = 0.0
  74         and then Re (Right) < 0.0
  75       then
  76          raise Constraint_Error;
  77
  78       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
  79          return Left;
  80
  81       elsif Right = (0.0, 0.0) then
  82          return Complex_One;
  83
  84       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
  85          return 1.0 + Right;
  86
  87       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
  88          return Left;
  89
  90       else
  91          return Exp (Right * Log (Left));
  92       end if;
  93    end "**";
  94
  95    function "**" (Left : Real'Base; Right : Complex) return Complex is
  96    begin
  97       if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
  98          raise Argument_Error;
  99
 100       elsif Left = 0.0 and then Re (Right) < 0.0 then
 101          raise Constraint_Error;
 102
 103       elsif Left = 0.0 then
 104          return Compose_From_Cartesian (Left, 0.0);
 105
 106       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
 107          return Complex_One;
 108
 109       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
 110          return Compose_From_Cartesian (Left, 0.0);
 111
 112       else
 113          return Exp (Log (Left) * Right);
 114       end if;
 115    end "**";
 116
 117    function "**" (Left : Complex; Right : Real'Base) return Complex is
 118    begin
 119       if Right = 0.0
 120         and then Re (Left) = 0.0
 121         and then Im (Left) = 0.0
 122       then
 123          raise Argument_Error;
 124
 125       elsif Re (Left) = 0.0
 126         and then Im (Left) = 0.0
 127         and then Right < 0.0
 128       then
 129          raise Constraint_Error;
 130
 131       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
 132          return Left;
 133
 134       elsif Right = 0.0 then
 135          return Complex_One;
 136
 137       elsif Right = 1.0 then
 138          return Left;
 139
 140       else
 141          return Exp (Right * Log (Left));
 142       end if;
 143    end "**";
 144
 145    ------------
 146    -- Arccos --
 147    ------------
 148
 149    function Arccos (X : Complex) return Complex is
 150       Result : Complex;
 151
 152    begin
 153       if X = Complex_One then
 154          return Complex_Zero;
 155
 156       elsif abs Re (X) < Square_Root_Epsilon and then
 157          abs Im (X) < Square_Root_Epsilon
 158       then
 159          return Half_Pi - X;
 160
 161       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 162             abs Im (X) > Inv_Square_Root_Epsilon
 163       then
 164          return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
 165                             Complex_I * Sqrt ((1.0 - X) / 2.0));
 166       end if;
 167
 168       Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
 169
 170       if Im (X) = 0.0
 171         and then abs Re (X) <= 1.00
 172       then
 173          Set_Im (Result, Im (X));
 174       end if;
 175
 176       return Result;
 177    end Arccos;
 178
 179    -------------
 180    -- Arccosh --
 181    -------------
 182
 183    function Arccosh (X : Complex) return Complex is
 184       Result : Complex;
 185
 186    begin
 187       if X = Complex_One then
 188          return Complex_Zero;
 189
 190       elsif abs Re (X) < Square_Root_Epsilon and then
 191          abs Im (X) < Square_Root_Epsilon
 192       then
 193          Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
 194
 195       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 196             abs Im (X) > Inv_Square_Root_Epsilon
 197       then
 198          Result := Log_Two + Log (X);
 199
 200       else
 201          Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
 202                               Sqrt ((X - 1.0) / 2.0));
 203       end if;
 204
 205       if Re (Result) <= 0.0 then
 206          Result := -Result;
 207       end if;
 208
 209       return Result;
 210    end Arccosh;
 211
 212    ------------
 213    -- Arccot --
 214    ------------
 215
 216    function Arccot (X : Complex) return Complex is
 217       Xt : Complex;
 218
 219    begin
 220       if abs Re (X) < Square_Root_Epsilon and then
 221          abs Im (X) < Square_Root_Epsilon
 222       then
 223          return Half_Pi - X;
 224
 225       elsif abs Re (X) > 1.0 / Epsilon or else
 226             abs Im (X) > 1.0 / Epsilon
 227       then
 228          Xt := Complex_One / X;
 229
 230          if Re (X) < 0.0 then
 231             Set_Re (Xt, PI - Re (Xt));
 232             return Xt;
 233          else
 234             return Xt;
 235          end if;
 236       end if;
 237
 238       Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
 239
 240       if Re (Xt) < 0.0 then
 241          Xt := PI + Xt;
 242       end if;
 243
 244       return Xt;
 245    end Arccot;
 246
 247    --------------
 248    -- Arccoth --
 249    --------------
 250
 251    function Arccoth (X : Complex) return Complex is
 252       R : Complex;
 253
 254    begin
 255       if X = (0.0, 0.0) then
 256          return Compose_From_Cartesian (0.0, PI_2);
 257
 258       elsif abs Re (X) < Square_Root_Epsilon
 259          and then abs Im (X) < Square_Root_Epsilon
 260       then
 261          return PI_2 * Complex_I + X;
 262
 263       elsif abs Re (X) > 1.0 / Epsilon or else
 264             abs Im (X) > 1.0 / Epsilon
 265       then
 266          if Im (X) > 0.0 then
 267             return (0.0, 0.0);
 268          else
 269             return PI * Complex_I;
 270          end if;
 271
 272       elsif Im (X) = 0.0 and then Re (X) = 1.0 then
 273          raise Constraint_Error;
 274
 275       elsif Im (X) = 0.0 and then Re (X) = -1.0 then
 276          raise Constraint_Error;
 277       end if;
 278
 279       begin
 280          R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
 281
 282       exception
 283          when Constraint_Error =>
 284             R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
 285       end;
 286
 287       if Im (R) < 0.0 then
 288          Set_Im (R, PI + Im (R));
 289       end if;
 290
 291       if Re (X) = 0.0 then
 292          Set_Re (R, Re (X));
 293       end if;
 294
 295       return R;
 296    end Arccoth;
 297
 298    ------------
 299    -- Arcsin --
 300    ------------
 301
 302    function Arcsin (X : Complex) return Complex is
 303       Result : Complex;
 304
 305    begin
 306       --  For very small argument, sin (x) = x
 307
 308       if abs Re (X) < Square_Root_Epsilon and then
 309          abs Im (X) < Square_Root_Epsilon
 310       then
 311          return X;
 312
 313       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 314             abs Im (X) > Inv_Square_Root_Epsilon
 315       then
 316          Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
 317
 318          if Im (Result) > PI_2 then
 319             Set_Im (Result, PI - Im (X));
 320
 321          elsif Im (Result) < -PI_2 then
 322             Set_Im (Result, -(PI + Im (X)));
 323          end if;
 324
 325          return Result;
 326       end if;
 327
 328       Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
 329
 330       if Re (X) = 0.0 then
 331          Set_Re (Result, Re (X));
 332
 333       elsif Im (X) = 0.0
 334         and then abs Re (X) <= 1.00
 335       then
 336          Set_Im (Result, Im (X));
 337       end if;
 338
 339       return Result;
 340    end Arcsin;
 341
 342    -------------
 343    -- Arcsinh --
 344    -------------
 345
 346    function Arcsinh (X : Complex) return Complex is
 347       Result : Complex;
 348
 349    begin
 350       if abs Re (X) < Square_Root_Epsilon and then
 351          abs Im (X) < Square_Root_Epsilon
 352       then
 353          return X;
 354
 355       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 356             abs Im (X) > Inv_Square_Root_Epsilon
 357       then
 358          Result := Log_Two + Log (X); -- may have wrong sign
 359
 360          if (Re (X) < 0.0 and then Re (Result) > 0.0)
 361            or else (Re (X) > 0.0 and then Re (Result) < 0.0)
 362          then
 363             Set_Re (Result, -Re (Result));
 364          end if;
 365
 366          return Result;
 367       end if;
 368
 369       Result := Log (X + Sqrt (1.0 + X * X));
 370
 371       if Re (X) = 0.0 then
 372          Set_Re (Result, Re (X));
 373       elsif Im  (X) = 0.0 then
 374          Set_Im (Result, Im  (X));
 375       end if;
 376
 377       return Result;
 378    end Arcsinh;
 379
 380    ------------
 381    -- Arctan --
 382    ------------
 383
 384    function Arctan (X : Complex) return Complex is
 385    begin
 386       if abs Re (X) < Square_Root_Epsilon and then
 387          abs Im (X) < Square_Root_Epsilon
 388       then
 389          return X;
 390
 391       else
 392          return -Complex_I * (Log (1.0 + Complex_I * X)
 393                             - Log (1.0 - Complex_I * X)) / 2.0;
 394       end if;
 395    end Arctan;
 396
 397    -------------
 398    -- Arctanh --
 399    -------------
 400
 401    function Arctanh (X : Complex) return Complex is
 402    begin
 403       if abs Re (X) < Square_Root_Epsilon and then
 404          abs Im (X) < Square_Root_Epsilon
 405       then
 406          return X;
 407       else
 408          return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
 409       end if;
 410    end Arctanh;
 411
 412    ---------
 413    -- Cos --
 414    ---------
 415
 416    function Cos (X : Complex) return Complex is
 417    begin
 418       return
 419         Compose_From_Cartesian
 420           (Cos (Re (X)) * Cosh (Im (X)),
 421            -(Sin (Re (X)) * Sinh (Im (X))));
 422    end Cos;
 423
 424    ----------
 425    -- Cosh --
 426    ----------
 427
 428    function Cosh (X : Complex) return Complex is
 429    begin
 430       return
 431         Compose_From_Cartesian
 432           (Cosh (Re (X)) * Cos (Im (X)),
 433            Sinh (Re (X)) * Sin (Im (X)));
 434    end Cosh;
 435
 436    ---------
 437    -- Cot --
 438    ---------
 439
 440    function Cot (X : Complex) return Complex is
 441    begin
 442       if abs Re (X) < Square_Root_Epsilon and then
 443          abs Im (X) < Square_Root_Epsilon
 444       then
 445          return Complex_One / X;
 446
 447       elsif Im (X) > Log_Inverse_Epsilon_2 then
 448          return -Complex_I;
 449
 450       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 451          return Complex_I;
 452       end if;
 453
 454       return Cos (X) / Sin (X);
 455    end Cot;
 456
 457    ----------
 458    -- Coth --
 459    ----------
 460
 461    function Coth (X : Complex) return Complex is
 462    begin
 463       if abs Re (X) < Square_Root_Epsilon and then
 464          abs Im (X) < Square_Root_Epsilon
 465       then
 466          return Complex_One / X;
 467
 468       elsif Re (X) > Log_Inverse_Epsilon_2 then
 469          return Complex_One;
 470
 471       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 472          return -Complex_One;
 473
 474       else
 475          return Cosh (X) / Sinh (X);
 476       end if;
 477    end Coth;
 478
 479    ---------
 480    -- Exp --
 481    ---------
 482
 483    function Exp (X : Complex) return Complex is
 484       ImX : constant Real'Base := Im (X);
 485       EXP_RE_X : constant Real'Base := Exp (Re (X));
 486
 487    begin
 488       return Compose_From_Cartesian (EXP_RE_X * Cos (ImX),
 489                                      EXP_RE_X * Sin (ImX));
 490    end Exp;
 491
 492    function Exp (X : Imaginary) return Complex is
 493       ImX : constant Real'Base := Im (X);
 494
 495    begin
 496       return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
 497    end Exp;
 498
 499    ---------
 500    -- Log --
 501    ---------
 502
 503    function Log (X : Complex) return Complex is
 504       ReX : Real'Base;
 505       ImX : Real'Base;
 506       Z   : Complex;
 507
 508    begin
 509       if Re (X) = 0.0 and then Im (X) = 0.0 then
 510          raise Constraint_Error;
 511
 512       elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
 513         and then abs Im (X) < Root_Root_Epsilon
 514       then
 515          Z := X;
 516          Set_Re (Z, Re (Z) - 1.0);
 517
 518          return (1.0 - (1.0 / 2.0 -
 519                        (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
 520       end if;
 521
 522       begin
 523          ReX := Log (Modulus (X));
 524
 525       exception
 526          when Constraint_Error =>
 527             ReX := Log (Modulus (X / 2.0)) - Log_Two;
 528       end;
 529
 530       ImX := Arctan (Im (X), Re (X));
 531
 532       if ImX > PI then
 533          ImX := ImX - 2.0 * PI;
 534       end if;
 535
 536       return Compose_From_Cartesian (ReX, ImX);
 537    end Log;
 538
 539    ---------
 540    -- Sin --
 541    ---------
 542
 543    function Sin (X : Complex) return Complex is
 544    begin
 545       if abs Re (X) < Square_Root_Epsilon
 546            and then
 547          abs Im (X) < Square_Root_Epsilon
 548       then
 549          return X;
 550       end if;
 551
 552       return
 553         Compose_From_Cartesian
 554           (Sin (Re (X)) * Cosh (Im (X)),
 555            Cos (Re (X)) * Sinh (Im (X)));
 556    end Sin;
 557
 558    ----------
 559    -- Sinh --
 560    ----------
 561
 562    function Sinh (X : Complex) return Complex is
 563    begin
 564       if abs Re (X) < Square_Root_Epsilon and then
 565          abs Im (X) < Square_Root_Epsilon
 566       then
 567          return X;
 568
 569       else
 570          return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
 571                                         Cosh (Re (X)) * Sin (Im (X)));
 572       end if;
 573    end Sinh;
 574
 575    ----------
 576    -- Sqrt --
 577    ----------
 578
 579    function Sqrt (X : Complex) return Complex is
 580       ReX : constant Real'Base := Re (X);
 581       ImX : constant Real'Base := Im (X);
 582       XR  : constant Real'Base := abs Re (X);
 583       YR  : constant Real'Base := abs Im (X);
 584       R   : Real'Base;
 585       R_X : Real'Base;
 586       R_Y : Real'Base;
 587
 588    begin
 589       --  Deal with pure real case, see (RM G.1.2(39))
 590
 591       if ImX = 0.0 then
 592          if ReX > 0.0 then
 593             return
 594               Compose_From_Cartesian
 595                 (Sqrt (ReX), 0.0);
 596
 597          elsif ReX = 0.0 then
 598             return X;
 599
 600          else
 601             return
 602               Compose_From_Cartesian
 603                 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
 604          end if;
 605
 606       elsif ReX = 0.0 then
 607          R_X := Sqrt (YR / 2.0);
 608
 609          if ImX > 0.0 then
 610             return Compose_From_Cartesian (R_X, R_X);
 611          else
 612             return Compose_From_Cartesian (R_X, -R_X);
 613          end if;
 614
 615       else
 616          R  := Sqrt (XR ** 2 + YR ** 2);
 617
 618          --  If the square of the modulus overflows, try rescaling the
 619          --  real and imaginary parts. We cannot depend on an exception
 620          --  being raised on all targets.
 621
 622          if R > Real'Base'Last then
 623             raise Constraint_Error;
 624          end if;
 625
 626          --  We are solving the system
 627
 628          --  XR = R_X ** 2 - Y_R ** 2      (1)
 629          --  YR = 2.0 * R_X * R_Y          (2)
 630          --
 631          --  The symmetric solution involves square roots for both R_X and
 632          --  R_Y, but it is more accurate to use the square root with the
 633          --  larger argument for either R_X or R_Y, and equation (2) for the
 634          --  other.
 635
 636          if ReX < 0.0 then
 637             R_Y := Sqrt (0.5 * (R - ReX));
 638             R_X := YR / (2.0 * R_Y);
 639
 640          else
 641             R_X := Sqrt (0.5 * (R + ReX));
 642             R_Y := YR / (2.0 * R_X);
 643          end if;
 644       end if;
 645
 646       if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 647          R_Y := -R_Y;
 648       end if;
 649       return Compose_From_Cartesian (R_X, R_Y);
 650
 651    exception
 652       when Constraint_Error =>
 653
 654          --  Rescale and try again
 655
 656          R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
 657          R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
 658          R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
 659
 660          if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 661             R_Y := -R_Y;
 662          end if;
 663
 664          return Compose_From_Cartesian (R_X, R_Y);
 665    end Sqrt;
 666
 667    ---------
 668    -- Tan --
 669    ---------
 670
 671    function Tan (X : Complex) return Complex is
 672    begin
 673       if abs Re (X) < Square_Root_Epsilon and then
 674          abs Im (X) < Square_Root_Epsilon
 675       then
 676          return X;
 677
 678       elsif Im (X) > Log_Inverse_Epsilon_2 then
 679          return Complex_I;
 680
 681       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 682          return -Complex_I;
 683
 684       else
 685          return Sin (X) / Cos (X);
 686       end if;
 687    end Tan;
 688
 689    ----------
 690    -- Tanh --
 691    ----------
 692
 693    function Tanh (X : Complex) return Complex is
 694    begin
 695       if abs Re (X) < Square_Root_Epsilon and then
 696          abs Im (X) < Square_Root_Epsilon
 697       then
 698          return X;
 699
 700       elsif Re (X) > Log_Inverse_Epsilon_2 then
 701          return Complex_One;
 702
 703       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 704          return -Complex_One;
 705
 706       else
 707          return Sinh (X) / Cosh (X);
 708       end if;
 709    end Tanh;
 710
 711 end Ada.Numerics.Generic_Complex_Elementary_Functions;