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-2001 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       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       end if;
 325
 326       Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
 327
 328       if Re (X) = 0.0 then
 329          Set_Re (Result, Re (X));
 330
 331       elsif Im (X) = 0.0
 332         and then abs Re (X) <= 1.00
 333       then
 334          Set_Im (Result, Im (X));
 335       end if;
 336
 337       return Result;
 338    end Arcsin;
 339
 340    -------------
 341    -- Arcsinh --
 342    -------------
 343
 344    function Arcsinh (X : Complex) return Complex is
 345       Result : Complex;
 346
 347    begin
 348       if abs Re (X) < Square_Root_Epsilon and then
 349          abs Im (X) < Square_Root_Epsilon
 350       then
 351          return X;
 352
 353       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 354             abs Im (X) > Inv_Square_Root_Epsilon
 355       then
 356          Result := Log_Two + Log (X); -- may have wrong sign
 357
 358          if (Re (X) < 0.0 and Re (Result) > 0.0)
 359            or else (Re (X) > 0.0 and Re (Result) < 0.0)
 360          then
 361             Set_Re (Result, -Re (Result));
 362          end if;
 363
 364          return Result;
 365       end if;
 366
 367       Result := Log (X + Sqrt (1.0 + X * X));
 368
 369       if Re (X) = 0.0 then
 370          Set_Re (Result, Re (X));
 371       elsif Im  (X) = 0.0 then
 372          Set_Im (Result, Im  (X));
 373       end if;
 374
 375       return Result;
 376    end Arcsinh;
 377
 378    ------------
 379    -- Arctan --
 380    ------------
 381
 382    function Arctan (X : Complex) return Complex is
 383    begin
 384       if abs Re (X) < Square_Root_Epsilon and then
 385          abs Im (X) < Square_Root_Epsilon
 386       then
 387          return X;
 388
 389       else
 390          return -Complex_I * (Log (1.0 + Complex_I * X)
 391                             - Log (1.0 - Complex_I * X)) / 2.0;
 392       end if;
 393    end Arctan;
 394
 395    -------------
 396    -- Arctanh --
 397    -------------
 398
 399    function Arctanh (X : Complex) return Complex is
 400    begin
 401       if abs Re (X) < Square_Root_Epsilon and then
 402          abs Im (X) < Square_Root_Epsilon
 403       then
 404          return X;
 405       else
 406          return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
 407       end if;
 408    end Arctanh;
 409
 410    ---------
 411    -- Cos --
 412    ---------
 413
 414    function Cos (X : Complex) return Complex is
 415    begin
 416       return
 417         Compose_From_Cartesian
 418           (Cos (Re (X))  * Cosh (Im (X)),
 419            -Sin (Re (X)) * Sinh (Im (X)));
 420    end Cos;
 421
 422    ----------
 423    -- Cosh --
 424    ----------
 425
 426    function Cosh (X : Complex) return Complex is
 427    begin
 428       return
 429         Compose_From_Cartesian
 430           (Cosh (Re (X)) * Cos (Im (X)),
 431            Sinh (Re (X)) * Sin (Im (X)));
 432    end Cosh;
 433
 434    ---------
 435    -- Cot --
 436    ---------
 437
 438    function Cot (X : Complex) return Complex is
 439    begin
 440       if abs Re (X) < Square_Root_Epsilon and then
 441          abs Im (X) < Square_Root_Epsilon
 442       then
 443          return Complex_One  /  X;
 444
 445       elsif Im (X) > Log_Inverse_Epsilon_2 then
 446          return -Complex_I;
 447
 448       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 449          return Complex_I;
 450       end if;
 451
 452       return Cos (X) / Sin (X);
 453    end Cot;
 454
 455    ----------
 456    -- Coth --
 457    ----------
 458
 459    function Coth (X : Complex) return Complex is
 460    begin
 461       if abs Re (X) < Square_Root_Epsilon and then
 462          abs Im (X) < Square_Root_Epsilon
 463       then
 464          return Complex_One  /  X;
 465
 466       elsif Re (X) > Log_Inverse_Epsilon_2 then
 467          return Complex_One;
 468
 469       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 470          return -Complex_One;
 471
 472       else
 473          return Cosh (X) / Sinh (X);
 474       end if;
 475    end Coth;
 476
 477    ---------
 478    -- Exp --
 479    ---------
 480
 481    function Exp (X : Complex) return Complex is
 482       EXP_RE_X : Real'Base := Exp (Re (X));
 483
 484    begin
 485       return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
 486                                      EXP_RE_X * Sin (Im (X)));
 487    end Exp;
 488
 489
 490    function Exp (X : Imaginary) return Complex is
 491       ImX : Real'Base := Im (X);
 492
 493    begin
 494       return Compose_From_Cartesian (Cos (ImX), Sin (ImX));
 495    end Exp;
 496
 497    ---------
 498    -- Log --
 499    ---------
 500
 501    function Log (X : Complex) return Complex is
 502       ReX : Real'Base;
 503       ImX : Real'Base;
 504       Z   : Complex;
 505
 506    begin
 507       if Re (X) = 0.0 and then Im (X) = 0.0 then
 508          raise Constraint_Error;
 509
 510       elsif abs (1.0 - Re (X)) < Root_Root_Epsilon
 511         and then abs Im (X) < Root_Root_Epsilon
 512       then
 513          Z := X;
 514          Set_Re (Z, Re (Z) - 1.0);
 515
 516          return (1.0 - (1.0 / 2.0 -
 517                        (1.0 / 3.0 - (1.0 / 4.0) * Z) * Z) * Z) * Z;
 518       end if;
 519
 520       begin
 521          ReX := Log (Modulus (X));
 522
 523       exception
 524          when Constraint_Error =>
 525             ReX := Log (Modulus (X / 2.0)) - Log_Two;
 526       end;
 527
 528       ImX := Arctan (Im (X), Re (X));
 529
 530       if ImX > PI then
 531          ImX := ImX - 2.0 * PI;
 532       end if;
 533
 534       return Compose_From_Cartesian (ReX, ImX);
 535    end Log;
 536
 537    ---------
 538    -- Sin --
 539    ---------
 540
 541    function Sin (X : Complex) return Complex is
 542    begin
 543       if abs Re (X) < Square_Root_Epsilon and then
 544          abs Im (X) < Square_Root_Epsilon then
 545          return X;
 546       end if;
 547
 548       return
 549         Compose_From_Cartesian
 550           (Sin (Re (X)) * Cosh (Im (X)),
 551            Cos (Re (X)) * Sinh (Im (X)));
 552    end Sin;
 553
 554    ----------
 555    -- Sinh --
 556    ----------
 557
 558    function Sinh (X : Complex) return Complex is
 559    begin
 560       if abs Re (X) < Square_Root_Epsilon and then
 561          abs Im (X) < Square_Root_Epsilon
 562       then
 563          return X;
 564
 565       else
 566          return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
 567                                         Cosh (Re (X)) * Sin (Im (X)));
 568       end if;
 569    end Sinh;
 570
 571    ----------
 572    -- Sqrt --
 573    ----------
 574
 575    function Sqrt (X : Complex) return Complex is
 576       ReX : constant Real'Base := Re (X);
 577       ImX : constant Real'Base := Im (X);
 578       XR  : constant Real'Base := abs Re (X);
 579       YR  : constant Real'Base := abs Im (X);
 580       R   : Real'Base;
 581       R_X : Real'Base;
 582       R_Y : Real'Base;
 583
 584    begin
 585       --  Deal with pure real case, see (RM G.1.2(39))
 586
 587       if ImX = 0.0 then
 588          if ReX > 0.0 then
 589             return
 590               Compose_From_Cartesian
 591                 (Sqrt (ReX), 0.0);
 592
 593          elsif ReX = 0.0 then
 594             return X;
 595
 596          else
 597             return
 598               Compose_From_Cartesian
 599                 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
 600          end if;
 601
 602       elsif ReX = 0.0 then
 603          R_X := Sqrt (YR / 2.0);
 604
 605          if ImX > 0.0 then
 606             return Compose_From_Cartesian (R_X, R_X);
 607          else
 608             return Compose_From_Cartesian (R_X, -R_X);
 609          end if;
 610
 611       else
 612          R  := Sqrt (XR ** 2 + YR ** 2);
 613
 614          --  If the square of the modulus overflows, try rescaling the
 615          --  real and imaginary parts. We cannot depend on an exception
 616          --  being raised on all targets.
 617
 618          if R > Real'Base'Last then
 619             raise Constraint_Error;
 620          end if;
 621
 622          --  We are solving the system
 623
 624          --  XR = R_X ** 2 - Y_R ** 2      (1)
 625          --  YR = 2.0 * R_X * R_Y          (2)
 626          --
 627          --  The symmetric solution involves square roots for both R_X and
 628          --  R_Y, but it is more accurate to use the square root with the
 629          --  larger argument for either R_X or R_Y, and equation (2) for the
 630          --  other.
 631
 632          if ReX < 0.0 then
 633             R_Y := Sqrt (0.5 * (R - ReX));
 634             R_X := YR / (2.0 * R_Y);
 635
 636          else
 637             R_X := Sqrt (0.5 * (R + ReX));
 638             R_Y := YR / (2.0 * R_X);
 639          end if;
 640       end if;
 641
 642       if Im (X) < 0.0 then                 -- halve angle, Sqrt of magnitude
 643          R_Y := -R_Y;
 644       end if;
 645       return Compose_From_Cartesian (R_X, R_Y);
 646
 647    exception
 648       when Constraint_Error =>
 649
 650          --  Rescale and try again.
 651
 652          R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
 653          R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
 654          R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
 655
 656          if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 657             R_Y := -R_Y;
 658          end if;
 659
 660          return Compose_From_Cartesian (R_X, R_Y);
 661    end Sqrt;
 662
 663    ---------
 664    -- Tan --
 665    ---------
 666
 667    function Tan (X : Complex) return Complex is
 668    begin
 669       if abs Re (X) < Square_Root_Epsilon and then
 670          abs Im (X) < Square_Root_Epsilon
 671       then
 672          return X;
 673
 674       elsif Im (X) > Log_Inverse_Epsilon_2 then
 675          return Complex_I;
 676
 677       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 678          return -Complex_I;
 679
 680       else
 681          return Sin (X) / Cos (X);
 682       end if;
 683    end Tan;
 684
 685    ----------
 686    -- Tanh --
 687    ----------
 688
 689    function Tanh (X : Complex) return Complex is
 690    begin
 691       if abs Re (X) < Square_Root_Epsilon and then
 692          abs Im (X) < Square_Root_Epsilon
 693       then
 694          return X;
 695
 696       elsif Re (X) > Log_Inverse_Epsilon_2 then
 697          return Complex_One;
 698
 699       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 700          return -Complex_One;
 701
 702       else
 703          return Sinh (X) / Cosh (X);
 704       end if;
 705    end Tanh;
 706
 707 end Ada.Numerics.Generic_Complex_Elementary_Functions;