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