gcc/ada/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-2009, 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       EXP_RE_X : constant Real'Base := Exp (Re (X));
 485
 486    begin
 487       return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
 488                                      EXP_RE_X * Sin (Im (X)));
 489    end Exp;
 490
 491    function Exp (X : Imaginary) return Complex is
 492       ImX : constant Real'Base := Im (X);
 493
 494    begin
 495       return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
 496    end Exp;
 497
 498    ---------
 499    -- Log --
 500    ---------
 501
 502    function Log (X : Complex) return Complex is
 503       ReX : Real'Base;
 504       ImX : Real'Base;
 505       Z   : Complex;
 506
 507    begin
 508       if Re (X) = 0.0 and then Im (X) = 0.0 then
 509          raise Constraint_Error;
 510
 511       elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
 512         and then abs Im (X) < Root_Root_Epsilon
 513       then
 514          Z := X;
 515          Set_Re (Z, Re (Z) - 1.0);
 516
 517          return (1.0 - (1.0 / 2.0 -
 518                        (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
 519       end if;
 520
 521       begin
 522          ReX := Log (Modulus (X));
 523
 524       exception
 525          when Constraint_Error =>
 526             ReX := Log (Modulus (X / 2.0)) - Log_Two;
 527       end;
 528
 529       ImX := Arctan (Im (X), Re (X));
 530
 531       if ImX > PI then
 532          ImX := ImX - 2.0 * PI;
 533       end if;
 534
 535       return Compose_From_Cartesian (ReX, ImX);
 536    end Log;
 537
 538    ---------
 539    -- Sin --
 540    ---------
 541
 542    function Sin (X : Complex) return Complex is
 543    begin
 544       if abs Re (X) < Square_Root_Epsilon and then
 545          abs Im (X) < Square_Root_Epsilon then
 546          return X;
 547       end if;
 548
 549       return
 550         Compose_From_Cartesian
 551           (Sin (Re (X)) * Cosh (Im (X)),
 552            Cos (Re (X)) * Sinh (Im (X)));
 553    end Sin;
 554
 555    ----------
 556    -- Sinh --
 557    ----------
 558
 559    function Sinh (X : Complex) return Complex is
 560    begin
 561       if abs Re (X) < Square_Root_Epsilon and then
 562          abs Im (X) < Square_Root_Epsilon
 563       then
 564          return X;
 565
 566       else
 567          return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
 568                                         Cosh (Re (X)) * Sin (Im (X)));
 569       end if;
 570    end Sinh;
 571
 572    ----------
 573    -- Sqrt --
 574    ----------
 575
 576    function Sqrt (X : Complex) return Complex is
 577       ReX : constant Real'Base := Re (X);
 578       ImX : constant Real'Base := Im (X);
 579       XR  : constant Real'Base := abs Re (X);
 580       YR  : constant Real'Base := abs Im (X);
 581       R   : Real'Base;
 582       R_X : Real'Base;
 583       R_Y : Real'Base;
 584
 585    begin
 586       --  Deal with pure real case, see (RM G.1.2(39))
 587
 588       if ImX = 0.0 then
 589          if ReX > 0.0 then
 590             return
 591               Compose_From_Cartesian
 592                 (Sqrt (ReX), 0.0);
 593
 594          elsif ReX = 0.0 then
 595             return X;
 596
 597          else
 598             return
 599               Compose_From_Cartesian
 600                 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
 601          end if;
 602
 603       elsif ReX = 0.0 then
 604          R_X := Sqrt (YR / 2.0);
 605
 606          if ImX > 0.0 then
 607             return Compose_From_Cartesian (R_X, R_X);
 608          else
 609             return Compose_From_Cartesian (R_X, -R_X);
 610          end if;
 611
 612       else
 613          R  := Sqrt (XR ** 2 + YR ** 2);
 614
 615          --  If the square of the modulus overflows, try rescaling the
 616          --  real and imaginary parts. We cannot depend on an exception
 617          --  being raised on all targets.
 618
 619          if R > Real'Base'Last then
 620             raise Constraint_Error;
 621          end if;
 622
 623          --  We are solving the system
 624
 625          --  XR = R_X ** 2 - Y_R ** 2      (1)
 626          --  YR = 2.0 * R_X * R_Y          (2)
 627          --
 628          --  The symmetric solution involves square roots for both R_X and
 629          --  R_Y, but it is more accurate to use the square root with the
 630          --  larger argument for either R_X or R_Y, and equation (2) for the
 631          --  other.
 632
 633          if ReX < 0.0 then
 634             R_Y := Sqrt (0.5 * (R - ReX));
 635             R_X := YR / (2.0 * R_Y);
 636
 637          else
 638             R_X := Sqrt (0.5 * (R + ReX));
 639             R_Y := YR / (2.0 * R_X);
 640          end if;
 641       end if;
 642
 643       if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 644          R_Y := -R_Y;
 645       end if;
 646       return Compose_From_Cartesian (R_X, R_Y);
 647
 648    exception
 649       when Constraint_Error =>
 650
 651          --  Rescale and try again
 652
 653          R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
 654          R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
 655          R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
 656
 657          if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 658             R_Y := -R_Y;
 659          end if;
 660
 661          return Compose_From_Cartesian (R_X, R_Y);
 662    end Sqrt;
 663
 664    ---------
 665    -- Tan --
 666    ---------
 667
 668    function Tan (X : Complex) return Complex is
 669    begin
 670       if abs Re (X) < Square_Root_Epsilon and then
 671          abs Im (X) < Square_Root_Epsilon
 672       then
 673          return X;
 674
 675       elsif Im (X) > Log_Inverse_Epsilon_2 then
 676          return Complex_I;
 677
 678       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 679          return -Complex_I;
 680
 681       else
 682          return Sin (X) / Cos (X);
 683       end if;
 684    end Tan;
 685
 686    ----------
 687    -- Tanh --
 688    ----------
 689
 690    function Tanh (X : Complex) return Complex is
 691    begin
 692       if abs Re (X) < Square_Root_Epsilon and then
 693          abs Im (X) < Square_Root_Epsilon
 694       then
 695          return X;
 696
 697       elsif Re (X) > Log_Inverse_Epsilon_2 then
 698          return Complex_One;
 699
 700       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 701          return -Complex_One;
 702
 703       else
 704          return Sinh (X) / Cosh (X);
 705       end if;
 706    end Tanh;
 707
 708 end Ada.Numerics.Generic_Complex_Elementary_Functions;