gcc/ada/math_lib.adb

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