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 --                                                                          --
  10 --          Copyright (C) 1992-2001 Free Software Foundation, Inc.
  11 --                                                                          --
  12 -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  13 -- terms of the  GNU General Public License as published  by the Free Soft- --
  14 -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  15 -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  16 -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  17 -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  18 -- for  more details.  You should have  received  a copy of the GNU General --
  19 -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  20 -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  21 -- MA 02111-1307, USA.                                                      --
  22 --                                                                          --
  23 -- As a special exception,  if other files  instantiate  generics from this --
  24 -- unit, or you link  this unit with other files  to produce an executable, --
  25 -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  26 -- covered  by the  GNU  General  Public  License.  This exception does not --
  27 -- however invalidate  any other reasons why  the executable file  might be --
  28 -- covered by the  GNU Public License.                                      --
  29 --                                                                          --
  30 -- GNAT was originally developed  by the GNAT team at  New York University. --
  31 -- Extensive contributions were provided by Ada Core Technologies Inc.      --
  32 --                                                                          --
  33 ------------------------------------------------------------------------------
  34
  35 with Ada.Numerics.Generic_Elementary_Functions;
  36
  37 package body Ada.Numerics.Generic_Complex_Elementary_Functions is
  38
  39    package Elementary_Functions is new
  40       Ada.Numerics.Generic_Elementary_Functions (Real'Base);
  41    use Elementary_Functions;
  42
  43    PI      : constant := 3.14159_26535_89793_23846_26433_83279_50288_41971;
  44    PI_2    : constant := PI / 2.0;
  45    Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
  46    Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  47
  48    subtype T is Real'Base;
  49
  50    Epsilon                 : constant T := 2.0      ** (1 - T'Model_Mantissa);
  51    Square_Root_Epsilon     : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
  52    Inv_Square_Root_Epsilon : constant T := Sqrt_Two ** (T'Model_Mantissa - 1);
  53    Root_Root_Epsilon       : constant T := Sqrt_Two **
  54                                                  ((1 - T'Model_Mantissa) / 2);
  55    Log_Inverse_Epsilon_2   : constant T := T (T'Model_Mantissa - 1) / 2.0;
  56
  57    Complex_Zero : constant Complex := (0.0,  0.0);
  58    Complex_One  : constant Complex := (1.0,  0.0);
  59    Complex_I    : constant Complex := (0.0,  1.0);
  60    Half_Pi      : constant Complex := (PI_2, 0.0);
  61
  62    --------
  63    -- ** --
  64    --------
  65
  66    function "**" (Left : Complex; Right : Complex) return Complex is
  67    begin
  68       if Re (Right) = 0.0
  69         and then Im (Right) = 0.0
  70         and then Re (Left)  = 0.0
  71         and then Im (Left)  = 0.0
  72       then
  73          raise Argument_Error;
  74
  75       elsif Re (Left) = 0.0
  76         and then Im (Left) = 0.0
  77         and then Re (Right) < 0.0
  78       then
  79          raise Constraint_Error;
  80
  81       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
  82          return Left;
  83
  84       elsif Right = (0.0, 0.0)  then
  85          return Complex_One;
  86
  87       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
  88          return 1.0 + Right;
  89
  90       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
  91          return Left;
  92
  93       else
  94          return Exp (Right * Log (Left));
  95       end if;
  96    end "**";
  97
  98    function "**" (Left : Real'Base; Right : Complex) return Complex is
  99    begin
 100       if Re (Right) = 0.0 and then Im (Right) = 0.0 and then Left = 0.0 then
 101          raise Argument_Error;
 102
 103       elsif Left = 0.0 and then Re (Right) < 0.0 then
 104          raise Constraint_Error;
 105
 106       elsif Left = 0.0 then
 107          return Compose_From_Cartesian (Left, 0.0);
 108
 109       elsif Re (Right) = 0.0 and then Im (Right) = 0.0 then
 110          return Complex_One;
 111
 112       elsif Re (Right) = 1.0 and then Im (Right) = 0.0 then
 113          return Compose_From_Cartesian (Left, 0.0);
 114
 115       else
 116          return Exp (Log (Left) * Right);
 117       end if;
 118    end "**";
 119
 120    function "**" (Left : Complex; Right : Real'Base) return Complex is
 121    begin
 122       if Right = 0.0
 123         and then Re (Left) = 0.0
 124         and then Im (Left) = 0.0
 125       then
 126          raise Argument_Error;
 127
 128       elsif Re (Left) = 0.0
 129         and then Im (Left) = 0.0
 130         and then Right < 0.0
 131       then
 132          raise Constraint_Error;
 133
 134       elsif Re (Left) = 0.0 and then Im (Left) = 0.0 then
 135          return Left;
 136
 137       elsif Right = 0.0 then
 138          return Complex_One;
 139
 140       elsif Right = 1.0 then
 141          return Left;
 142
 143       else
 144          return Exp (Right * Log (Left));
 145       end if;
 146    end "**";
 147
 148    ------------
 149    -- Arccos --
 150    ------------
 151
 152    function Arccos (X : Complex) return Complex is
 153       Result : Complex;
 154
 155    begin
 156       if X = Complex_One then
 157          return Complex_Zero;
 158
 159       elsif abs Re (X) < Square_Root_Epsilon and then
 160          abs Im (X) < Square_Root_Epsilon
 161       then
 162          return Half_Pi - X;
 163
 164       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 165             abs Im (X) > Inv_Square_Root_Epsilon
 166       then
 167          return -2.0 * Complex_I * Log (Sqrt ((1.0 + X) / 2.0) +
 168                             Complex_I * Sqrt ((1.0 - X) / 2.0));
 169       end if;
 170
 171       Result := -Complex_I * Log (X + Complex_I * Sqrt (1.0 - X * X));
 172
 173       if Im (X) = 0.0
 174         and then abs Re (X) <= 1.00
 175       then
 176          Set_Im (Result, Im (X));
 177       end if;
 178
 179       return Result;
 180    end Arccos;
 181
 182    -------------
 183    -- Arccosh --
 184    -------------
 185
 186    function Arccosh (X : Complex) return Complex is
 187       Result : Complex;
 188
 189    begin
 190       if X = Complex_One then
 191          return Complex_Zero;
 192
 193       elsif abs Re (X) < Square_Root_Epsilon and then
 194          abs Im (X) < Square_Root_Epsilon
 195       then
 196          Result := Compose_From_Cartesian (-Im (X), -PI_2 + Re (X));
 197
 198       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 199             abs Im (X) > Inv_Square_Root_Epsilon
 200       then
 201          Result := Log_Two + Log (X);
 202
 203       else
 204          Result := 2.0 * Log (Sqrt ((1.0 + X) / 2.0) +
 205                               Sqrt ((X - 1.0) / 2.0));
 206       end if;
 207
 208       if Re (Result) <= 0.0 then
 209          Result := -Result;
 210       end if;
 211
 212       return Result;
 213    end Arccosh;
 214
 215    ------------
 216    -- Arccot --
 217    ------------
 218
 219    function Arccot (X : Complex) return Complex is
 220       Xt : Complex;
 221
 222    begin
 223       if abs Re (X) < Square_Root_Epsilon and then
 224          abs Im (X) < Square_Root_Epsilon
 225       then
 226          return Half_Pi - X;
 227
 228       elsif abs Re (X) > 1.0 / Epsilon or else
 229             abs Im (X) > 1.0 / Epsilon
 230       then
 231          Xt := Complex_One  /  X;
 232
 233          if Re (X) < 0.0 then
 234             Set_Re (Xt, PI - Re (Xt));
 235             return Xt;
 236          else
 237             return Xt;
 238          end if;
 239       end if;
 240
 241       Xt := Complex_I * Log ((X - Complex_I) / (X + Complex_I)) / 2.0;
 242
 243       if Re (Xt) < 0.0 then
 244          Xt := PI + Xt;
 245       end if;
 246
 247       return Xt;
 248    end Arccot;
 249
 250    --------------
 251    -- Arctcoth --
 252    --------------
 253
 254    function Arccoth (X : Complex) return Complex is
 255       R : Complex;
 256
 257    begin
 258       if X = (0.0, 0.0) then
 259          return Compose_From_Cartesian (0.0, PI_2);
 260
 261       elsif abs Re (X) < Square_Root_Epsilon
 262          and then abs Im (X) < Square_Root_Epsilon
 263       then
 264          return PI_2 * Complex_I + X;
 265
 266       elsif abs Re (X) > 1.0 / Epsilon or else
 267             abs Im (X) > 1.0 / Epsilon
 268       then
 269          if Im (X) > 0.0 then
 270             return (0.0, 0.0);
 271          else
 272             return PI * Complex_I;
 273          end if;
 274
 275       elsif Im (X) = 0.0 and then Re (X) = 1.0 then
 276          raise Constraint_Error;
 277
 278       elsif Im (X) = 0.0 and then Re (X) = -1.0 then
 279          raise Constraint_Error;
 280       end if;
 281
 282       begin
 283          R := Log ((1.0 + X) / (X - 1.0)) / 2.0;
 284
 285       exception
 286          when Constraint_Error =>
 287             R := (Log (1.0 + X) - Log (X - 1.0)) / 2.0;
 288       end;
 289
 290       if Im (R) < 0.0 then
 291          Set_Im (R, PI + Im (R));
 292       end if;
 293
 294       if Re (X) = 0.0 then
 295          Set_Re (R, Re (X));
 296       end if;
 297
 298       return R;
 299    end Arccoth;
 300
 301    ------------
 302    -- Arcsin --
 303    ------------
 304
 305    function Arcsin (X : Complex) return Complex is
 306       Result : Complex;
 307
 308    begin
 309       if abs Re (X) < Square_Root_Epsilon and then
 310          abs Im (X) < Square_Root_Epsilon
 311       then
 312          return X;
 313
 314       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 315             abs Im (X) > Inv_Square_Root_Epsilon
 316       then
 317          Result := -Complex_I * (Log (Complex_I * X) + Log (2.0 * Complex_I));
 318
 319          if Im (Result) > PI_2 then
 320             Set_Im (Result, PI - Im (X));
 321
 322          elsif Im (Result) < -PI_2 then
 323             Set_Im (Result, -(PI + Im (X)));
 324          end if;
 325       end if;
 326
 327       Result := -Complex_I * Log (Complex_I * X + Sqrt (1.0 - X * X));
 328
 329       if Re (X) = 0.0 then
 330          Set_Re (Result, Re (X));
 331
 332       elsif Im (X) = 0.0
 333         and then abs Re (X) <= 1.00
 334       then
 335          Set_Im (Result, Im (X));
 336       end if;
 337
 338       return Result;
 339    end Arcsin;
 340
 341    -------------
 342    -- Arcsinh --
 343    -------------
 344
 345    function Arcsinh (X : Complex) return Complex is
 346       Result : Complex;
 347
 348    begin
 349       if abs Re (X) < Square_Root_Epsilon and then
 350          abs Im (X) < Square_Root_Epsilon
 351       then
 352          return X;
 353
 354       elsif abs Re (X) > Inv_Square_Root_Epsilon or else
 355             abs Im (X) > Inv_Square_Root_Epsilon
 356       then
 357          Result := Log_Two + Log (X); -- may have wrong sign
 358
 359          if (Re (X) < 0.0 and Re (Result) > 0.0)
 360            or else (Re (X) > 0.0 and Re (Result) < 0.0)
 361          then
 362             Set_Re (Result, -Re (Result));
 363          end if;
 364
 365          return Result;
 366       end if;
 367
 368       Result := Log (X + Sqrt (1.0 + X * X));
 369
 370       if Re (X) = 0.0 then
 371          Set_Re (Result, Re (X));
 372       elsif Im  (X) = 0.0 then
 373          Set_Im (Result, Im  (X));
 374       end if;
 375
 376       return Result;
 377    end Arcsinh;
 378
 379    ------------
 380    -- Arctan --
 381    ------------
 382
 383    function Arctan (X : Complex) return Complex is
 384    begin
 385       if abs Re (X) < Square_Root_Epsilon and then
 386          abs Im (X) < Square_Root_Epsilon
 387       then
 388          return X;
 389
 390       else
 391          return -Complex_I * (Log (1.0 + Complex_I * X)
 392                             - Log (1.0 - Complex_I * X)) / 2.0;
 393       end if;
 394    end Arctan;
 395
 396    -------------
 397    -- Arctanh --
 398    -------------
 399
 400    function Arctanh (X : Complex) return Complex is
 401    begin
 402       if abs Re (X) < Square_Root_Epsilon and then
 403          abs Im (X) < Square_Root_Epsilon
 404       then
 405          return X;
 406       else
 407          return (Log (1.0 + X) - Log (1.0 - X)) / 2.0;
 408       end if;
 409    end Arctanh;
 410
 411    ---------
 412    -- Cos --
 413    ---------
 414
 415    function Cos (X : Complex) return Complex is
 416    begin
 417       return
 418         Compose_From_Cartesian
 419           (Cos (Re (X))  * Cosh (Im (X)),
 420            -Sin (Re (X)) * Sinh (Im (X)));
 421    end Cos;
 422
 423    ----------
 424    -- Cosh --
 425    ----------
 426
 427    function Cosh (X : Complex) return Complex is
 428    begin
 429       return
 430         Compose_From_Cartesian
 431           (Cosh (Re (X)) * Cos (Im (X)),
 432            Sinh (Re (X)) * Sin (Im (X)));
 433    end Cosh;
 434
 435    ---------
 436    -- Cot --
 437    ---------
 438
 439    function Cot (X : Complex) return Complex is
 440    begin
 441       if abs Re (X) < Square_Root_Epsilon and then
 442          abs Im (X) < Square_Root_Epsilon
 443       then
 444          return Complex_One  /  X;
 445
 446       elsif Im (X) > Log_Inverse_Epsilon_2 then
 447          return -Complex_I;
 448
 449       elsif Im (X) < -Log_Inverse_Epsilon_2 then
 450          return Complex_I;
 451       end if;
 452
 453       return Cos (X) / Sin (X);
 454    end Cot;
 455
 456    ----------
 457    -- Coth --
 458    ----------
 459
 460    function Coth (X : Complex) return Complex is
 461    begin
 462       if abs Re (X) < Square_Root_Epsilon and then
 463          abs Im (X) < Square_Root_Epsilon
 464       then
 465          return Complex_One  /  X;
 466
 467       elsif Re (X) > Log_Inverse_Epsilon_2 then
 468          return Complex_One;
 469
 470       elsif Re (X) < -Log_Inverse_Epsilon_2 then
 471          return -Complex_One;
 472
 473       else
 474          return Cosh (X) / Sinh (X);
 475       end if;
 476    end Coth;
 477
 478    ---------
 479    -- Exp --
 480    ---------
 481
 482    function Exp (X : Complex) return Complex is
 483       EXP_RE_X : Real'Base := Exp (Re (X));
 484
 485    begin
 486       return Compose_From_Cartesian (EXP_RE_X * Cos (Im (X)),
 487                                      EXP_RE_X * Sin (Im (X)));
 488    end Exp;
 489
 490
 491    function Exp (X : Imaginary) return Complex is
 492       ImX : 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;