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 --                            $Revision: 1.13 $
  10 --                                                                          --
  11 --          Copyright (C) 1992-2001 Free Software Foundation, Inc.
  12 --                                                                          --
  13 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  14 -- terms of the  GNU General Public License as published  by the Free Soft- --
  15 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  16 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  17 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  18 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  19 -- for  more details.  You should have  received  a copy of the GNU General --
  20 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  21 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  22 -- MA 02111-1307, USA.                                                      --
  23 --                                                                          --
  24 -- As a special exception,  if other files  instantiate  generics from this --
  25 -- unit, or you link  this unit with other files  to produce an executable, --
  26 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  27 -- covered  by the  GNU  General  Public  License.  This exception does not --
  28 -- however invalidate  any other reasons why  the executable file  might be --
  29 -- covered by the  GNU Public License.                                      --
  30 --                                                                          --
  31 -- GNAT was originally developed  by the GNAT team at  New York University. --
  32 -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
  33 --                                                                          --
  34 ------------------------------------------------------------------------------
  35
  36 with Ada.Numerics.Generic_Elementary_Functions;
  37
  38 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
  39
  40    package Elementary_Functions is new
  41       Ada.Numerics.Generic_Elementary_Functions (Real'Base);
  42    use Elementary_Functions;
  43
  44    PI      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
  45    PI_2    : constant := PI / 2.0;
  46    Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
  47    Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  48
  49    subtype T is Real'Base;
  50
  51    Epsilon                 : constant T := 2.0      ** (1 - T'Model_Mantissa);
  52    Square_Root_Epsilon     : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
  53    Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
  54    Root_Root_Epsilon       : constant T := Sqrt_Two **
  55                                                  ((1 - T'Model_Mantissa) / 2);
  56    Log_Inverse_Epsilon_2   : constant T := T (T'Model_Mantissa - 1) / 2.0;
  57
  58    Complex_Zero : constant Complex := (0.0,  0.0);
  59    Complex_One  : constant Complex := (1.0,  0.0);
  60    Complex_I    : constant Complex := (0.0,  1.0);
  61    Half_Pi      : constant Complex := (PI_2, 0.0);
  62
  63    --------
  64    -- ** --
  65    --------
  66
  67    function "**" (Left : Complex; Right : Complex) return Complex is
  68    begin
  69       if Re (Right) = 0.0
  70         and then Im (Right) = 0.0
  71         and then Re (Left)  = 0.0
  72         and then Im (Left)  = 0.0
  73       then
  74          raise Argument_Error;
  75
  76       elsif Re (Left) = 0.0
  77         and then Im (Left) = 0.0
  78         and then Re (Right) < 0.0
  79       then
  80          raise Constraint_Error;
  81
  82       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
  83          return Left;
  84
  85       elsif Right = (0.0, 0.0)  then
  86          return Complex_One;
  87
  88       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
  89          return 1.0 + Right;
  90
  91       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
  92          return Left;
  93
  94       else
  95          return Exp (Right * Log (Left));
  96       end if;
  97    end "**";
  98
  99    function "**" (Left : Real'Base; Right : Complex) return Complex is
 100    begin
 101       if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
 102          raise Argument_Error;
 103
 104       elsif Left = 0.0 and then Re (Right) < 0.0 then
 105          raise Constraint_Error;
 106
 107       elsif Left = 0.0 then
 108          return Compose_From_Cartesian (Left, 0.0);
 109
 110       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
 111          return Complex_One;
 112
 113       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
 114          return Compose_From_Cartesian (Left, 0.0);
 115
 116       else
 117          return Exp (Log (Left) * Right);
 118       end if;
 119    end "**";
 120
 121    function "**" (Left : Complex; Right : Real'Base) return Complex is
 122    begin
 123       if Right = 0.0
 124         and then Re (Left) = 0.0
 125         and then Im (Left) = 0.0
 126       then
 127          raise Argument_Error;
 128
 129       elsif Re (Left) = 0.0
 130         and then Im (Left) = 0.0
 131         and then Right < 0.0
 132       then
 133          raise Constraint_Error;
 134
 135       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
 136          return Left;
 137
 138       elsif Right = 0.0 then
 139          return Complex_One;
 140
 141       elsif Right = 1.0 then
 142          return Left;
 143
 144       else
 145          return Exp (Right * Log (Left));
 146       end if;
 147    end "**";
 148
 149    ------------
 150    -- Arccos --
 151    ------------
 152
 153    function Arccos (X : Complex) return Complex is
 154       Result : Complex;
 155
 156    begin
 157       if X = Complex_One then
 158          return Complex_Zero;
 159
 160       elsif abs Re (X) < Square_Root_Epsilon and then
 161          abs Im (X) < Square_Root_Epsilon
 162       then
 163          return Half_Pi - X;
 164
 165       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 166             abs Im (X) > Inv_Square_Root_Epsilon
 167       then
 168          return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
 169                             Complex_I * Sqrt ((1.0 - X) / 2.0));
 170       end if;
 171
 172       Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
 173
 174       if Im (X) = 0.0
 175         and then abs Re (X) <= 1.00
 176       then
 177          Set_Im (Result, Im (X));
 178       end if;
 179
 180       return Result;
 181    end Arccos;
 182
 183    -------------
 184    -- Arccosh --
 185    -------------
 186
 187    function Arccosh (X : Complex) return Complex is
 188       Result : Complex;
 189
 190    begin
 191       if X = Complex_One then
 192          return Complex_Zero;
 193
 194       elsif abs Re (X) < Square_Root_Epsilon and then
 195          abs Im (X) < Square_Root_Epsilon
 196       then
 197          Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
 198
 199       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 200             abs Im (X) > Inv_Square_Root_Epsilon
 201       then
 202          Result := Log_Two + Log (X);
 203
 204       else
 205          Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
 206                               Sqrt ((X - 1.0) / 2.0));
 207       end if;
 208
 209       if Re (Result) <= 0.0 then
 210          Result := -Result;
 211       end if;
 212
 213       return Result;
 214    end Arccosh;
 215
 216    ------------
 217    -- Arccot --
 218    ------------
 219
 220    function Arccot (X : Complex) return Complex is
 221       Xt : Complex;
 222
 223    begin
 224       if abs Re (X) < Square_Root_Epsilon and then
 225          abs Im (X) < Square_Root_Epsilon
 226       then
 227          return Half_Pi - X;
 228
 229       elsif abs Re (X) > 1.0 / Epsilon or else
 230             abs Im (X) > 1.0 / Epsilon
 231       then
 232          Xt := Complex_One  /  X;
 233
 234          if Re (X) < 0.0 then
 235             Set_Re (Xt, PI - Re (Xt));
 236             return Xt;
 237          else
 238             return Xt;
 239          end if;
 240       end if;
 241
 242       Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
 243
 244       if Re (Xt) < 0.0 then
 245          Xt := PI + Xt;
 246       end if;
 247
 248       return Xt;
 249    end Arccot;
 250
 251    --------------
 252    -- Arctcoth --
 253    --------------
 254
 255    function Arccoth (X : Complex) return Complex is
 256       R : Complex;
 257
 258    begin
 259       if X = (0.0, 0.0) then
 260          return Compose_From_Cartesian (0.0, PI_2);
 261
 262       elsif abs Re (X) < Square_Root_Epsilon
 263          and then abs Im (X) < Square_Root_Epsilon
 264       then
 265          return PI_2 * Complex_I + X;
 266
 267       elsif abs Re (X) > 1.0 / Epsilon or else
 268             abs Im (X) > 1.0 / Epsilon
 269       then
 270          if Im (X) > 0.0 then
 271             return (0.0, 0.0);
 272          else
 273             return PI * Complex_I;
 274          end if;
 275
 276       elsif Im (X) = 0.0 and then Re (X) = 1.0 then
 277          raise Constraint_Error;
 278
 279       elsif Im (X) = 0.0 and then Re (X) = -1.0 then
 280          raise Constraint_Error;
 281       end if;
 282
 283       begin
 284          R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
 285
 286       exception
 287          when Constraint_Error =>
 288             R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
 289       end;
 290
 291       if Im (R) < 0.0 then
 292          Set_Im (R, PI + Im (R));
 293       end if;
 294
 295       if Re (X) = 0.0 then
 296          Set_Re (R, Re (X));
 297       end if;
 298
 299       return R;
 300    end Arccoth;
 301
 302    ------------
 303    -- Arcsin --
 304    ------------
 305
 306    function Arcsin (X : Complex) return Complex is
 307       Result : Complex;
 308
 309    begin
 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       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 Re (Result) > 0.0)
 361            or else (Re (X) > 0.0 and 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 : 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
 492    function Exp (X : Imaginary) return Complex is
 493       ImX : 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 and then
 546          abs Im (X) < Square_Root_Epsilon then
 547          return X;
 548       end if;
 549
 550       return
 551         Compose_From_Cartesian
 552           (Sin (Re (X)) * Cosh (Im (X)),
 553            Cos (Re (X)) * Sinh (Im (X)));
 554    end Sin;
 555
 556    ----------
 557    -- Sinh --
 558    ----------
 559
 560    function Sinh (X : Complex) return Complex is
 561    begin
 562       if abs Re (X) < Square_Root_Epsilon and then
 563          abs Im (X) < Square_Root_Epsilon
 564       then
 565          return X;
 566
 567       else
 568          return Compose_From_Cartesian (Sinh (Re (X)) * Cos (Im (X)),
 569                                         Cosh (Re (X)) * Sin (Im (X)));
 570       end if;
 571    end Sinh;
 572
 573    ----------
 574    -- Sqrt --
 575    ----------
 576
 577    function Sqrt (X : Complex) return Complex is
 578       ReX : constant Real'Base := Re (X);
 579       ImX : constant Real'Base := Im (X);
 580       XR  : constant Real'Base := abs Re (X);
 581       YR  : constant Real'Base := abs Im (X);
 582       R   : Real'Base;
 583       R_X : Real'Base;
 584       R_Y : Real'Base;
 585
 586    begin
 587       --  Deal with pure real case, see (RM G.1.2(39))
 588
 589       if ImX = 0.0 then
 590          if ReX > 0.0 then
 591             return
 592               Compose_From_Cartesian
 593                 (Sqrt (ReX), 0.0);
 594
 595          elsif ReX = 0.0 then
 596             return X;
 597
 598          else
 599             return
 600               Compose_From_Cartesian
 601                 (0.0, Real'Copy_Sign (Sqrt (-ReX), ImX));
 602          end if;
 603
 604       elsif ReX = 0.0 then
 605          R_X := Sqrt (YR / 2.0);
 606
 607          if ImX > 0.0 then
 608             return Compose_From_Cartesian (R_X, R_X);
 609          else
 610             return Compose_From_Cartesian (R_X, -R_X);
 611          end if;
 612
 613       else
 614          R  := Sqrt (XR ** 2 + YR ** 2);
 615
 616          --  If the square of the modulus overflows, try rescaling the
 617          --  real and imaginary parts. We cannot depend on an exception
 618          --  being raised on all targets.
 619
 620          if R > Real'Base'Last then
 621             raise Constraint_Error;
 622          end if;
 623
 624          --  We are solving the system
 625
 626          --  XR = R_X ** 2 - Y_R ** 2      (1)
 627          --  YR = 2.0 * R_X * R_Y          (2)
 628          --
 629          --  The symmetric solution involves square roots for both R_X and
 630          --  R_Y, but it is more accurate to use the square root with the
 631          --  larger argument for either R_X or R_Y, and equation (2) for the
 632          --  other.
 633
 634          if ReX < 0.0 then
 635             R_Y := Sqrt (0.5 * (R - ReX));
 636             R_X := YR / (2.0 * R_Y);
 637
 638          else
 639             R_X := Sqrt (0.5 * (R + ReX));
 640             R_Y := YR / (2.0 * R_X);
 641          end if;
 642       end if;
 643
 644       if Im (X) < 0.0 then                 -- halve angle, Sqrt of magnitude
 645          R_Y := -R_Y;
 646       end if;
 647       return Compose_From_Cartesian (R_X, R_Y);
 648
 649    exception
 650       when Constraint_Error =>
 651
 652          --  Rescale and try again.
 653
 654          R := Modulus (Compose_From_Cartesian (Re (X / 4.0), Im (X / 4.0)));
 655          R_X := 2.0 * Sqrt (0.5 * R + 0.5 * Re (X / 4.0));
 656          R_Y := 2.0 * Sqrt (0.5 * R - 0.5 * Re (X / 4.0));
 657
 658          if Im (X) < 0.0 then -- halve angle, Sqrt of magnitude
 659             R_Y := -R_Y;
 660          end if;
 661
 662          return Compose_From_Cartesian (R_X, R_Y);
 663    end Sqrt;
 664
 665    ---------
 666    -- Tan --
 667    ---------
 668
 669    function Tan (X : Complex) return Complex is
 670    begin
 671       if abs Re (X) < Square_Root_Epsilon and then
 672          abs Im (X) < Square_Root_Epsilon
 673       then
 674          return X;
 675
 676       elsif Im (X) > Log_Inverse_Epsilon_2 then
 677          return Complex_I;
 678
 679       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 680          return -Complex_I;
 681
 682       else
 683          return Sin (X) / Cos (X);
 684       end if;
 685    end Tan;
 686
 687    ----------
 688    -- Tanh --
 689    ----------
 690
 691    function Tanh (X : Complex) return Complex is
 692    begin
 693       if abs Re (X) < Square_Root_Epsilon and then
 694          abs Im (X) < Square_Root_Epsilon
 695       then
 696          return X;
 697
 698       elsif Re (X) > Log_Inverse_Epsilon_2 then
 699          return Complex_One;
 700
 701       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 702          return -Complex_One;
 703
 704       else
 705          return Sinh (X) / Cosh (X);
 706       end if;
 707    end Tanh;
 708
 709 end Ada.Numerics.Generic_Complex_Elementary_Functions;