gcc/ada/math_lib.adb

   1 ------------------------------------------------------------------------------
   2 --                                                                          --
   3 --                         GNAT RUNTIME COMPONENTS                          --
   4 --                                                                          --
   5 --                             M A T H _ L I B                              --
   6 --                                                                          --
   7 --                                 B o d y                                  --
   8 --                                                                          --
   9 --                                                                          --
  10 --          Copyright (C) 1992-2000 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 --  This body is specifically for using an Ada interface to C math.h to get
  36 --  the computation engine. Many special cases are handled locally to avoid
  37 --  unnecessary calls. This is not a "strict" implementation, but takes full
  38 --  advantage of the C functions, e.g. in providing interface to hardware
  39 --  provided versions of the elementary functions.
  40
  41 --  A known weakness is that on the x86, all computation is done in Double,
  42 --  which means that a lot of accuracy is lost for the Long_Long_Float case.
  43
  44 --  Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
  45 --  sinh, cosh, tanh from C library via math.h
  46
  47 --  This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
  48 --  provides a compatible body for the DEC Math_Lib package.
  49
  50 with Ada.Numerics.Aux;
  51 use type Ada.Numerics.Aux.Double;
  52 with Ada.Numerics; use Ada.Numerics;
  53
  54 package body Math_Lib is
  55
  56    Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
  57
  58    Two_Pi     : constant Real'Base := 2.0 * Pi;
  59    Half_Pi    : constant Real'Base := Pi / 2.0;
  60    Fourth_Pi  : constant Real'Base := Pi / 4.0;
  61    Epsilon    : constant Real'Base := Real'Base'Epsilon;
  62    IEpsilon   : constant Real'Base := 1.0 / Epsilon;
  63
  64    subtype Double is Aux.Double;
  65
  66    DEpsilon    : constant Double := Double (Epsilon);
  67    DIEpsilon   : constant Double := Double (IEpsilon);
  68
  69    -----------------------
  70    -- Local Subprograms --
  71    -----------------------
  72
  73    function Arctan
  74      (Y    : Real;
  75       A    : Real := 1.0)
  76       return Real;
  77
  78    function Arctan
  79      (Y     : Real;
  80       A     : Real := 1.0;
  81       Cycle : Real)
  82       return  Real;
  83
  84    function Exact_Remainder
  85      (A    : Real;
  86       Y    : Real)
  87       return Real;
  88    --  Computes exact remainder of A divided by Y
  89
  90    function Half_Log_Epsilon return Real;
  91    --  Function to provide constant: 0.5 * Log (Epsilon)
  92
  93    function Local_Atan
  94      (Y    : Real;
  95       A    : Real := 1.0)
  96       return Real;
  97    --  Common code for arc tangent after cyele reduction
  98
  99    function Log_Inverse_Epsilon return Real;
 100    --  Function to provide constant: Log (1.0 / Epsilon)
 101
 102    function Square_Root_Epsilon return Real;
 103    --  Function to provide constant: Sqrt (Epsilon)
 104
 105    ----------
 106    -- "**" --
 107    ----------
 108
 109    function "**" (A1, A2 : Real) return Real is
 110
 111    begin
 112       if A1 = 0.0
 113         and then A2 = 0.0
 114       then
 115          raise Argument_Error;
 116
 117       elsif A1 < 0.0 then
 118          raise Argument_Error;
 119
 120       elsif A2 = 0.0 then
 121          return 1.0;
 122
 123       elsif A1 = 0.0 then
 124          if A2 < 0.0 then
 125             raise Constraint_Error;
 126          else
 127             return 0.0;
 128          end if;
 129
 130       elsif A1 = 1.0 then
 131          return 1.0;
 132
 133       elsif A2 = 1.0 then
 134          return A1;
 135
 136       else
 137          begin
 138             if A2 = 2.0 then
 139                return A1 * A1;
 140             else
 141                return
 142                  Real (Aux.pow (Double (A1), Double (A2)));
 143             end if;
 144
 145          exception
 146             when others =>
 147                raise Constraint_Error;
 148          end;
 149       end if;
 150    end "**";
 151
 152    ------------
 153    -- Arccos --
 154    ------------
 155
 156    --  Natural cycle
 157
 158    function Arccos (A : Real) return Real is
 159       Temp : Real'Base;
 160
 161    begin
 162       if abs A > 1.0 then
 163          raise Argument_Error;
 164
 165       elsif abs A < Square_Root_Epsilon then
 166          return Pi / 2.0 - A;
 167
 168       elsif A = 1.0 then
 169          return 0.0;
 170
 171       elsif A = -1.0 then
 172          return Pi;
 173       end if;
 174
 175       Temp := Real (Aux.acos (Double (A)));
 176
 177       if Temp < 0.0 then
 178          Temp := Pi + Temp;
 179       end if;
 180
 181       return Temp;
 182    end Arccos;
 183
 184    --  Arbitrary cycle
 185
 186    function Arccos (A, Cycle : Real) return Real is
 187       Temp : Real'Base;
 188
 189    begin
 190       if Cycle <= 0.0 then
 191          raise Argument_Error;
 192
 193       elsif abs A > 1.0 then
 194          raise Argument_Error;
 195
 196       elsif abs A < Square_Root_Epsilon then
 197          return Cycle / 4.0;
 198
 199       elsif A = 1.0 then
 200          return 0.0;
 201
 202       elsif A = -1.0 then
 203          return Cycle / 2.0;
 204       end if;
 205
 206       Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
 207
 208       if Temp < 0.0 then
 209          Temp := Cycle / 2.0 + Temp;
 210       end if;
 211
 212       return Temp;
 213    end Arccos;
 214
 215    -------------
 216    -- Arccosh --
 217    -------------
 218
 219    function Arccosh (A : Real) return Real is
 220    begin
 221       --  Return Log (A - Sqrt (A * A - 1.0));  double valued,
 222       --    only positive value returned
 223       --  What is this comment ???
 224
 225       if A < 1.0 then
 226          raise Argument_Error;
 227
 228       elsif A < 1.0 + Square_Root_Epsilon then
 229          return A - 1.0;
 230
 231       elsif abs A > 1.0 / Square_Root_Epsilon then
 232          return Log (A) + Log_Two;
 233
 234       else
 235          return Log (A + Sqrt (A * A - 1.0));
 236       end if;
 237    end Arccosh;
 238
 239    ------------
 240    -- Arccot --
 241    ------------
 242
 243    --  Natural cycle
 244
 245    function Arccot
 246      (A    : Real;
 247       Y    : Real := 1.0)
 248       return Real
 249    is
 250    begin
 251       --  Just reverse arguments
 252
 253       return Arctan (Y, A);
 254    end Arccot;
 255
 256    --  Arbitrary cycle
 257
 258    function Arccot
 259      (A     : Real;
 260       Y     : Real := 1.0;
 261       Cycle : Real)
 262       return  Real
 263    is
 264    begin
 265       --  Just reverse arguments
 266
 267       return Arctan (Y, A, Cycle);
 268    end Arccot;
 269
 270    -------------
 271    -- Arccoth --
 272    -------------
 273
 274    function Arccoth (A : Real) return Real is
 275    begin
 276       if abs A = 1.0 then
 277          raise Constraint_Error;
 278
 279       elsif abs A < 1.0 then
 280          raise Argument_Error;
 281
 282       elsif abs A > 1.0 / Epsilon then
 283          return 0.0;
 284
 285       else
 286          return 0.5 * Log ((1.0 + A) / (A - 1.0));
 287       end if;
 288    end Arccoth;
 289
 290    ------------
 291    -- Arcsin --
 292    ------------
 293
 294    --  Natural cycle
 295
 296    function Arcsin (A : Real) return Real is
 297    begin
 298       if abs A > 1.0 then
 299          raise Argument_Error;
 300
 301       elsif abs A < Square_Root_Epsilon then
 302          return A;
 303
 304       elsif A = 1.0 then
 305          return Pi / 2.0;
 306
 307       elsif A = -1.0 then
 308          return -Pi / 2.0;
 309       end if;
 310
 311       return Real (Aux.asin (Double (A)));
 312    end Arcsin;
 313
 314    --  Arbitrary cycle
 315
 316    function Arcsin (A, Cycle : Real) return Real is
 317    begin
 318       if Cycle <= 0.0 then
 319          raise Argument_Error;
 320
 321       elsif abs A > 1.0 then
 322          raise Argument_Error;
 323
 324       elsif A = 0.0 then
 325          return A;
 326
 327       elsif A = 1.0 then
 328          return Cycle / 4.0;
 329
 330       elsif A = -1.0 then
 331          return -Cycle / 4.0;
 332       end if;
 333
 334       return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
 335    end Arcsin;
 336
 337    -------------
 338    -- Arcsinh --
 339    -------------
 340
 341    function Arcsinh (A : Real) return Real is
 342    begin
 343       if abs A < Square_Root_Epsilon then
 344          return A;
 345
 346       elsif A > 1.0 / Square_Root_Epsilon then
 347          return Log (A) + Log_Two;
 348
 349       elsif A < -1.0 / Square_Root_Epsilon then
 350          return -(Log (-A) + Log_Two);
 351
 352       elsif A < 0.0 then
 353          return -Log (abs A + Sqrt (A * A + 1.0));
 354
 355       else
 356          return Log (A + Sqrt (A * A + 1.0));
 357       end if;
 358    end Arcsinh;
 359
 360    ------------
 361    -- Arctan --
 362    ------------
 363
 364    --  Natural cycle
 365
 366    function Arctan
 367      (Y    : Real;
 368       A    : Real := 1.0)
 369       return Real
 370    is
 371    begin
 372       if A = 0.0
 373         and then Y = 0.0
 374       then
 375          raise Argument_Error;
 376
 377       elsif Y = 0.0 then
 378          if A > 0.0 then
 379             return 0.0;
 380          else -- A < 0.0
 381             return Pi;
 382          end if;
 383
 384       elsif A = 0.0 then
 385          if Y > 0.0 then
 386             return Half_Pi;
 387          else -- Y < 0.0
 388             return -Half_Pi;
 389          end if;
 390
 391       else
 392          return Local_Atan (Y, A);
 393       end if;
 394    end Arctan;
 395
 396    --  Arbitrary cycle
 397
 398    function Arctan
 399      (Y     : Real;
 400       A     : Real := 1.0;
 401       Cycle : Real)
 402       return  Real
 403    is
 404    begin
 405       if Cycle <= 0.0 then
 406          raise Argument_Error;
 407
 408       elsif A = 0.0
 409         and then Y = 0.0
 410       then
 411          raise Argument_Error;
 412
 413       elsif Y = 0.0 then
 414          if A > 0.0 then
 415             return 0.0;
 416          else -- A < 0.0
 417             return Cycle / 2.0;
 418          end if;
 419
 420       elsif A = 0.0 then
 421          if Y > 0.0 then
 422             return Cycle / 4.0;
 423          else -- Y < 0.0
 424             return -Cycle / 4.0;
 425          end if;
 426
 427       else
 428          return Local_Atan (Y, A) *  Cycle / Two_Pi;
 429       end if;
 430    end Arctan;
 431
 432    -------------
 433    -- Arctanh --
 434    -------------
 435
 436    function Arctanh (A : Real) return Real is
 437    begin
 438       if abs A = 1.0 then
 439          raise Constraint_Error;
 440
 441       elsif abs A > 1.0 then
 442          raise Argument_Error;
 443
 444       elsif abs A < Square_Root_Epsilon then
 445          return A;
 446
 447       else
 448          return 0.5 * Log ((1.0 + A) / (1.0 - A));
 449       end if;
 450    end Arctanh;
 451
 452    ---------
 453    -- Cos --
 454    ---------
 455
 456    --  Natural cycle
 457
 458    function Cos (A : Real) return Real is
 459    begin
 460       if A = 0.0 then
 461          return 1.0;
 462
 463       elsif abs A < Square_Root_Epsilon then
 464          return 1.0;
 465
 466       end if;
 467
 468       return Real (Aux.Cos (Double (A)));
 469    end Cos;
 470
 471    --  Arbitrary cycle
 472
 473    function Cos (A, Cycle : Real) return Real is
 474       T : Real'Base;
 475
 476    begin
 477       if Cycle <= 0.0 then
 478          raise Argument_Error;
 479
 480       elsif A = 0.0 then
 481          return 1.0;
 482       end if;
 483
 484       T := Exact_Remainder (abs (A), Cycle) / Cycle;
 485
 486       if T = 0.25
 487         or else T = 0.75
 488         or else T = -0.25
 489         or else T = -0.75
 490       then
 491          return 0.0;
 492
 493       elsif T = 0.5 or T = -0.5 then
 494          return -1.0;
 495       end if;
 496
 497       return Real (Aux.Cos (Double (T * Two_Pi)));
 498    end Cos;
 499
 500    ----------
 501    -- Cosh --
 502    ----------
 503
 504    function Cosh (A : Real) return Real is
 505    begin
 506       if abs A < Square_Root_Epsilon then
 507          return 1.0;
 508
 509       elsif abs A > Log_Inverse_Epsilon then
 510          return Exp ((abs A) - Log_Two);
 511       end if;
 512
 513       return Real (Aux.cosh (Double (A)));
 514
 515    exception
 516       when others =>
 517          raise Constraint_Error;
 518    end Cosh;
 519
 520    ---------
 521    -- Cot --
 522    ---------
 523
 524    --  Natural cycle
 525
 526    function Cot (A : Real) return Real is
 527    begin
 528       if A = 0.0 then
 529          raise Constraint_Error;
 530
 531       elsif abs A < Square_Root_Epsilon then
 532          return 1.0 / A;
 533       end if;
 534
 535       return Real (1.0 / Real'Base (Aux.tan (Double (A))));
 536    end Cot;
 537
 538    --  Arbitrary cycle
 539
 540    function Cot (A, Cycle : Real) return Real is
 541       T : Real'Base;
 542
 543    begin
 544       if Cycle <= 0.0 then
 545          raise Argument_Error;
 546
 547       elsif A = 0.0 then
 548          raise Constraint_Error;
 549
 550       elsif abs A < Square_Root_Epsilon then
 551          return 1.0 / A;
 552       end if;
 553
 554       T := Exact_Remainder (A, Cycle) / Cycle;
 555
 556       if T = 0.0 or T = 0.5 or T = -0.5 then
 557          raise Constraint_Error;
 558       else
 559          return  Cos (T * Two_Pi) / Sin (T * Two_Pi);
 560       end if;
 561    end Cot;
 562
 563    ----------
 564    -- Coth --
 565    ----------
 566
 567    function Coth (A : Real) return Real is
 568    begin
 569       if A = 0.0 then
 570          raise Constraint_Error;
 571
 572       elsif A < Half_Log_Epsilon then
 573          return -1.0;
 574
 575       elsif A > -Half_Log_Epsilon then
 576          return 1.0;
 577
 578       elsif abs A < Square_Root_Epsilon then
 579          return 1.0 / A;
 580       end if;
 581
 582       return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
 583    end Coth;
 584
 585    ---------------------
 586    -- Exact_Remainder --
 587    ---------------------
 588
 589    function Exact_Remainder
 590      (A    : Real;
 591       Y    : Real)
 592       return Real
 593    is
 594       Denominator : Real'Base := abs A;
 595       Divisor     : Real'Base := abs Y;
 596       Reducer     : Real'Base;
 597       Sign        : Real'Base := 1.0;
 598
 599    begin
 600       if Y = 0.0 then
 601          raise Constraint_Error;
 602
 603       elsif A = 0.0 then
 604          return 0.0;
 605
 606       elsif A = Y then
 607          return 0.0;
 608
 609       elsif Denominator < Divisor then
 610          return A;
 611       end if;
 612
 613       while Denominator >= Divisor loop
 614
 615          --  Put divisors mantissa with denominators exponent to make reducer
 616
 617          Reducer := Divisor;
 618
 619          begin
 620             while Reducer * 1_048_576.0 < Denominator loop
 621                Reducer := Reducer * 1_048_576.0;
 622             end loop;
 623
 624          exception
 625             when others => null;
 626          end;
 627
 628          begin
 629             while Reducer * 1_024.0 < Denominator loop
 630                Reducer := Reducer * 1_024.0;
 631             end loop;
 632
 633          exception
 634             when others => null;
 635          end;
 636
 637          begin
 638             while Reducer * 2.0 < Denominator loop
 639                Reducer := Reducer * 2.0;
 640             end loop;
 641
 642          exception
 643             when others => null;
 644          end;
 645
 646          Denominator := Denominator - Reducer;
 647       end loop;
 648
 649       if A < 0.0 then
 650          return -Denominator;
 651       else
 652          return Denominator;
 653       end if;
 654    end Exact_Remainder;
 655
 656    ---------
 657    -- Exp --
 658    ---------
 659
 660    function Exp (A : Real) return Real is
 661       Result : Real'Base;
 662
 663    begin
 664       if A = 0.0 then
 665          return 1.0;
 666
 667       else
 668          Result := Real (Aux.Exp (Double (A)));
 669
 670          --  The check here catches the case of Exp returning IEEE infinity
 671
 672          if Result > Real'Last then
 673             raise Constraint_Error;
 674          else
 675             return Result;
 676          end if;
 677       end if;
 678    end Exp;
 679
 680    ----------------------
 681    -- Half_Log_Epsilon --
 682    ----------------------
 683
 684    --  Cannot precompute this constant, because this is required to be a
 685    --  pure package, which allows no state. A pity, but no way around it!
 686
 687    function Half_Log_Epsilon return Real is
 688    begin
 689       return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
 690    end Half_Log_Epsilon;
 691
 692    ----------------
 693    -- Local_Atan --
 694    ----------------
 695
 696    function Local_Atan
 697      (Y    : Real;
 698       A    : Real := 1.0)
 699       return Real
 700    is
 701       Z        : Real'Base;
 702       Raw_Atan : Real'Base;
 703
 704    begin
 705       if abs Y > abs A then
 706          Z := abs (A / Y);
 707       else
 708          Z := abs (Y / A);
 709       end if;
 710
 711       if Z < Square_Root_Epsilon then
 712          Raw_Atan := Z;
 713
 714       elsif Z = 1.0 then
 715          Raw_Atan := Pi / 4.0;
 716
 717       elsif Z < Square_Root_Epsilon then
 718          Raw_Atan := Z;
 719
 720       else
 721          Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
 722       end if;
 723
 724       if abs Y > abs A then
 725          Raw_Atan := Half_Pi - Raw_Atan;
 726       end if;
 727
 728       if A > 0.0 then
 729          if Y > 0.0 then
 730             return Raw_Atan;
 731          else                 --  Y < 0.0
 732             return -Raw_Atan;
 733          end if;
 734
 735       else                    --  A < 0.0
 736          if Y > 0.0 then
 737             return Pi - Raw_Atan;
 738          else                  --  Y < 0.0
 739             return -(Pi - Raw_Atan);
 740          end if;
 741       end if;
 742    end Local_Atan;
 743
 744    ---------
 745    -- Log --
 746    ---------
 747
 748    --  Natural base
 749
 750    function Log (A : Real) return Real is
 751    begin
 752       if A < 0.0 then
 753          raise Argument_Error;
 754
 755       elsif A = 0.0 then
 756          raise Constraint_Error;
 757
 758       elsif A = 1.0 then
 759          return 0.0;
 760       end if;
 761
 762       return Real (Aux.Log (Double (A)));
 763    end Log;
 764
 765    --  Arbitrary base
 766
 767    function Log (A, Base : Real) return Real is
 768    begin
 769       if A < 0.0 then
 770          raise Argument_Error;
 771
 772       elsif Base <= 0.0 or else Base = 1.0 then
 773          raise Argument_Error;
 774
 775       elsif A = 0.0 then
 776          raise Constraint_Error;
 777
 778       elsif A = 1.0 then
 779          return 0.0;
 780       end if;
 781
 782       return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
 783    end Log;
 784
 785    -------------------------
 786    -- Log_Inverse_Epsilon --
 787    -------------------------
 788
 789    --  Cannot precompute this constant, because this is required to be a
 790    --  pure package, which allows no state. A pity, but no way around it!
 791
 792    function Log_Inverse_Epsilon return Real is
 793    begin
 794       return Real (Aux.Log (DIEpsilon));
 795    end Log_Inverse_Epsilon;
 796
 797    ---------
 798    -- Sin --
 799    ---------
 800
 801    --  Natural cycle
 802
 803    function Sin (A : Real) return Real is
 804    begin
 805       if abs A < Square_Root_Epsilon then
 806          return A;
 807       end if;
 808
 809       return Real (Aux.Sin (Double (A)));
 810    end Sin;
 811
 812    --  Arbitrary cycle
 813
 814    function Sin (A, Cycle : Real) return Real is
 815       T : Real'Base;
 816
 817    begin
 818       if Cycle <= 0.0 then
 819          raise Argument_Error;
 820
 821       elsif A = 0.0 then
 822          return A;
 823       end if;
 824
 825       T := Exact_Remainder (A, Cycle) / Cycle;
 826
 827       if T = 0.0 or T = 0.5 or T = -0.5 then
 828          return 0.0;
 829
 830       elsif T = 0.25 or T = -0.75 then
 831          return 1.0;
 832
 833       elsif T = -0.25 or T = 0.75 then
 834          return -1.0;
 835
 836       end if;
 837
 838       return Real (Aux.Sin (Double (T * Two_Pi)));
 839    end Sin;
 840
 841    ----------
 842    -- Sinh --
 843    ----------
 844
 845    function Sinh (A : Real) return Real is
 846    begin
 847       if abs A < Square_Root_Epsilon then
 848          return A;
 849
 850       elsif  A > Log_Inverse_Epsilon then
 851          return Exp (A - Log_Two);
 852
 853       elsif A < -Log_Inverse_Epsilon then
 854          return -Exp ((-A) - Log_Two);
 855       end if;
 856
 857       return Real (Aux.Sinh (Double (A)));
 858
 859    exception
 860       when others =>
 861          raise Constraint_Error;
 862    end Sinh;
 863
 864    -------------------------
 865    -- Square_Root_Epsilon --
 866    -------------------------
 867
 868    --  Cannot precompute this constant, because this is required to be a
 869    --  pure package, which allows no state. A pity, but no way around it!
 870
 871    function Square_Root_Epsilon return Real is
 872    begin
 873       return Real (Aux.Sqrt (DEpsilon));
 874    end Square_Root_Epsilon;
 875
 876    ----------
 877    -- Sqrt --
 878    ----------
 879
 880    function Sqrt (A : Real) return Real is
 881    begin
 882       if A < 0.0 then
 883          raise Argument_Error;
 884
 885       --  Special case Sqrt (0.0) to preserve possible minus sign per IEEE
 886
 887       elsif A = 0.0 then
 888          return A;
 889
 890       --  Sqrt (1.0) must be exact for good complex accuracy
 891
 892       elsif A = 1.0 then
 893          return 1.0;
 894
 895       end if;
 896
 897       return Real (Aux.Sqrt (Double (A)));
 898    end Sqrt;
 899
 900    ---------
 901    -- Tan --
 902    ---------
 903
 904    --  Natural cycle
 905
 906    function Tan (A : Real) return Real is
 907    begin
 908       if abs A < Square_Root_Epsilon then
 909          return A;
 910
 911       elsif abs A = Pi / 2.0 then
 912          raise Constraint_Error;
 913       end if;
 914
 915       return Real (Aux.tan (Double (A)));
 916    end Tan;
 917
 918    --  Arbitrary cycle
 919
 920    function Tan (A, Cycle : Real) return Real is
 921       T : Real'Base;
 922
 923    begin
 924       if Cycle <= 0.0 then
 925          raise Argument_Error;
 926
 927       elsif A = 0.0 then
 928          return A;
 929       end if;
 930
 931       T := Exact_Remainder (A, Cycle) / Cycle;
 932
 933       if T = 0.25
 934         or else T = 0.75
 935         or else T = -0.25
 936         or else T = -0.75
 937       then
 938          raise Constraint_Error;
 939
 940       else
 941          return  Sin (T * Two_Pi) / Cos (T * Two_Pi);
 942       end if;
 943    end Tan;
 944
 945    ----------
 946    -- Tanh --
 947    ----------
 948
 949    function Tanh (A : Real) return Real is
 950    begin
 951       if A < Half_Log_Epsilon then
 952          return -1.0;
 953
 954       elsif A > -Half_Log_Epsilon then
 955          return 1.0;
 956
 957       elsif abs A < Square_Root_Epsilon then
 958          return A;
 959       end if;
 960
 961       return Real (Aux.tanh (Double (A)));
 962    end Tanh;
 963
 964    ----------------------------
 965    -- DEC-Specific functions --
 966    ----------------------------
 967
 968    function LOG10  (A : REAL) return REAL is
 969    begin
 970       return Log (A, 10.0);
 971    end LOG10;
 972
 973    function LOG2   (A : REAL) return REAL is
 974    begin
 975       return Log (A, 2.0);
 976    end LOG2;
 977
 978    function ASIN (A : REAL) return REAL renames Arcsin;
 979    function ACOS (A : REAL) return REAL renames Arccos;
 980
 981    function ATAN (A : REAL) return REAL is
 982    begin
 983       return Arctan (A, 1.0);
 984    end ATAN;
 985
 986    function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
 987
 988    function SIND   (A : REAL) return REAL is
 989    begin
 990       return Sin (A, 360.0);
 991    end SIND;
 992
 993    function COSD   (A : REAL) return REAL is
 994    begin
 995       return  Cos (A, 360.0);
 996    end COSD;
 997
 998    function TAND   (A : REAL) return REAL is
 999    begin
1000       return  Tan (A, 360.0);
1001    end TAND;
1002
1003    function ASIND  (A : REAL) return REAL is
1004    begin
1005       return  Arcsin (A, 360.0);
1006    end ASIND;
1007
1008    function ACOSD  (A : REAL) return REAL is
1009    begin
1010       return  Arccos (A, 360.0);
1011    end ACOSD;
1012
1013    function Arctan  (A : REAL) return REAL is
1014    begin
1015       return  Arctan (A, 1.0, 360.0);
1016    end Arctan;
1017
1018    function ATAND (A : REAL) return REAL is
1019    begin
1020       return Arctan (A, 1.0, 360.0);
1021    end ATAND;
1022
1023    function ATAN2D (A1, A2 : REAL) return REAL is
1024    begin
1025       return Arctan (A1, A2, 360.0);
1026    end ATAN2D;
1027
1028 end Math_Lib;