gcc/ada/a-ngcefu.adb

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