1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME COMPONENTS --
10 -- Copyright (C) 1992-2000 Free Software Foundation, Inc. --
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. --
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. --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- Extensive contributions were provided by Ada Core Technologies Inc. --
33 ------------------------------------------------------------------------------
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.
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.
44 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
45 -- sinh, cosh, tanh from C library via math.h
47 -- This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
48 -- provides a compatible body for the DEC Math_Lib package.
50 with Ada
.Numerics
.Aux
;
51 use type Ada
.Numerics
.Aux
.Double
;
52 with Ada
.Numerics
; use Ada
.Numerics
;
54 package body Math_Lib
is
56 Log_Two
: constant := 0.69314_71805_59945_30941_72321_21458_17656_80755
;
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;
64 subtype Double is Aux.Double;
66 DEpsilon : constant Double := Double (Epsilon);
67 DIEpsilon : constant Double := Double (IEpsilon);
69 -----------------------
70 -- Local Subprograms --
71 -----------------------
84 function Exact_Remainder
88 -- Computes exact remainder of A divided by Y
90 function Half_Log_Epsilon return Real;
91 -- Function to provide constant: 0.5 * Log (Epsilon)
97 -- Common code for arc tangent after cyele reduction
99 function Log_Inverse_Epsilon return Real;
100 -- Function to provide constant: Log (1.0 / Epsilon)
102 function Square_Root_Epsilon return Real;
103 -- Function to provide constant: Sqrt (Epsilon)
109 function "**" (A1, A2 : Real) return Real is
115 raise Argument_Error;
118 raise Argument_Error;
125 raise Constraint_Error;
142 Real (Aux.pow (Double (A1), Double (A2)));
147 raise Constraint_Error;
158 function Arccos (A : Real) return Real is
163 raise Argument_Error;
165 elsif abs A < Square_Root_Epsilon then
175 Temp := Real (Aux.acos (Double (A)));
186 function Arccos (A, Cycle : Real) return Real is
191 raise Argument_Error;
193 elsif abs A > 1.0 then
194 raise Argument_Error;
196 elsif abs A < Square_Root_Epsilon then
206 Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
209 Temp := Cycle / 2.0 + Temp;
219 function Arccosh (A : Real) return Real is
221 -- Return Log (A - Sqrt (A * A - 1.0)); double valued,
222 -- only positive value returned
223 -- What is this comment ???
226 raise Argument_Error;
228 elsif A < 1.0 + Square_Root_Epsilon then
231 elsif abs A > 1.0 / Square_Root_Epsilon then
232 return Log (A) + Log_Two;
235 return Log (A + Sqrt (A * A - 1.0));
251 -- Just reverse arguments
253 return Arctan (Y, A);
265 -- Just reverse arguments
267 return Arctan (Y, A, Cycle);
274 function Arccoth (A : Real) return Real is
277 raise Constraint_Error;
279 elsif abs A < 1.0 then
280 raise Argument_Error;
282 elsif abs A > 1.0 / Epsilon then
286 return 0.5 * Log ((1.0 + A) / (A - 1.0));
296 function Arcsin (A : Real) return Real is
299 raise Argument_Error;
301 elsif abs A < Square_Root_Epsilon then
311 return Real (Aux.asin (Double (A)));
316 function Arcsin (A, Cycle : Real) return Real is
319 raise Argument_Error;
321 elsif abs A > 1.0 then
322 raise Argument_Error;
334 return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
341 function Arcsinh (A : Real) return Real is
343 if abs A < Square_Root_Epsilon then
346 elsif A > 1.0 / Square_Root_Epsilon then
347 return Log (A) + Log_Two;
349 elsif A < -1.0 / Square_Root_Epsilon then
350 return -(Log (-A) + Log_Two);
353 return -Log (abs A + Sqrt (A * A + 1.0));
356 return Log (A + Sqrt (A * A + 1.0));
375 raise Argument_Error;
392 return Local_Atan (Y, A);
406 raise Argument_Error;
411 raise Argument_Error;
428 return Local_Atan (Y, A) * Cycle / Two_Pi;
436 function Arctanh (A : Real) return Real is
439 raise Constraint_Error;
441 elsif abs A > 1.0 then
442 raise Argument_Error;
444 elsif abs A < Square_Root_Epsilon then
448 return 0.5 * Log ((1.0 + A) / (1.0 - A));
458 function Cos (A : Real) return Real is
463 elsif abs A < Square_Root_Epsilon then
468 return Real (Aux.Cos (Double (A)));
473 function Cos (A, Cycle : Real) return Real is
478 raise Argument_Error;
484 T := Exact_Remainder (abs (A), Cycle) / Cycle;
493 elsif T = 0.5 or T = -0.5 then
497 return Real (Aux.Cos (Double (T * Two_Pi)));
504 function Cosh (A : Real) return Real is
506 if abs A < Square_Root_Epsilon then
509 elsif abs A > Log_Inverse_Epsilon then
510 return Exp ((abs A) - Log_Two);
513 return Real (Aux.cosh (Double (A)));
517 raise Constraint_Error;
526 function Cot (A : Real) return Real is
529 raise Constraint_Error;
531 elsif abs A < Square_Root_Epsilon then
535 return Real (1.0 / Real'Base (Aux.tan (Double (A))));
540 function Cot (A, Cycle : Real) return Real is
545 raise Argument_Error;
548 raise Constraint_Error;
550 elsif abs A < Square_Root_Epsilon then
554 T := Exact_Remainder (A, Cycle) / Cycle;
556 if T = 0.0 or T = 0.5 or T = -0.5 then
557 raise Constraint_Error;
559 return Cos (T * Two_Pi) / Sin (T * Two_Pi);
567 function Coth (A : Real) return Real is
570 raise Constraint_Error;
572 elsif A < Half_Log_Epsilon then
575 elsif A > -Half_Log_Epsilon then
578 elsif abs A < Square_Root_Epsilon then
582 return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
585 ---------------------
586 -- Exact_Remainder --
587 ---------------------
589 function Exact_Remainder
594 Denominator : Real'Base := abs A;
595 Divisor : Real'Base := abs Y;
597 Sign : Real'Base := 1.0;
601 raise Constraint_Error;
609 elsif Denominator < Divisor then
613 while Denominator >= Divisor loop
615 -- Put divisors mantissa with denominators exponent to make reducer
620 while Reducer * 1_048_576.0 < Denominator loop
621 Reducer := Reducer * 1_048_576.0;
629 while Reducer * 1_024.0 < Denominator loop
630 Reducer := Reducer * 1_024.0;
638 while Reducer * 2.0 < Denominator loop
639 Reducer := Reducer * 2.0;
646 Denominator := Denominator - Reducer;
660 function Exp (A : Real) return Real is
668 Result := Real (Aux.Exp (Double (A)));
670 -- The check here catches the case of Exp returning IEEE infinity
672 if Result > Real'Last then
673 raise Constraint_Error;
680 ----------------------
681 -- Half_Log_Epsilon --
682 ----------------------
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!
687 function Half_Log_Epsilon return Real is
689 return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
690 end Half_Log_Epsilon;
702 Raw_Atan : Real'Base;
705 if abs Y > abs A then
711 if Z < Square_Root_Epsilon then
715 Raw_Atan := Pi / 4.0;
717 elsif Z < Square_Root_Epsilon then
721 Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
724 if abs Y > abs A then
725 Raw_Atan := Half_Pi - Raw_Atan;
737 return Pi - Raw_Atan;
739 return -(Pi - Raw_Atan);
750 function Log (A : Real) return Real is
753 raise Argument_Error;
756 raise Constraint_Error;
762 return Real (Aux.Log (Double (A)));
767 function Log (A, Base : Real) return Real is
770 raise Argument_Error;
772 elsif Base <= 0.0 or else Base = 1.0 then
773 raise Argument_Error;
776 raise Constraint_Error;
782 return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
785 -------------------------
786 -- Log_Inverse_Epsilon --
787 -------------------------
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!
792 function Log_Inverse_Epsilon return Real is
794 return Real (Aux.Log (DIEpsilon));
795 end Log_Inverse_Epsilon;
803 function Sin (A : Real) return Real is
805 if abs A < Square_Root_Epsilon then
809 return Real (Aux.Sin (Double (A)));
814 function Sin (A, Cycle : Real) return Real is
819 raise Argument_Error;
825 T := Exact_Remainder (A, Cycle) / Cycle;
827 if T = 0.0 or T = 0.5 or T = -0.5 then
830 elsif T = 0.25 or T = -0.75 then
833 elsif T = -0.25 or T = 0.75 then
838 return Real (Aux.Sin (Double (T * Two_Pi)));
845 function Sinh (A : Real) return Real is
847 if abs A < Square_Root_Epsilon then
850 elsif A > Log_Inverse_Epsilon then
851 return Exp (A - Log_Two);
853 elsif A < -Log_Inverse_Epsilon then
854 return -Exp ((-A) - Log_Two);
857 return Real (Aux.Sinh (Double (A)));
861 raise Constraint_Error;
864 -------------------------
865 -- Square_Root_Epsilon --
866 -------------------------
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!
871 function Square_Root_Epsilon return Real is
873 return Real (Aux.Sqrt (DEpsilon));
874 end Square_Root_Epsilon;
880 function Sqrt (A : Real) return Real is
883 raise Argument_Error;
885 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
890 -- Sqrt (1.0) must be exact for good complex accuracy
897 return Real (Aux.Sqrt (Double (A)));
906 function Tan (A : Real) return Real is
908 if abs A < Square_Root_Epsilon then
911 elsif abs A = Pi / 2.0 then
912 raise Constraint_Error;
915 return Real (Aux.tan (Double (A)));
920 function Tan (A, Cycle : Real) return Real is
925 raise Argument_Error;
931 T := Exact_Remainder (A, Cycle) / Cycle;
938 raise Constraint_Error;
941 return Sin (T * Two_Pi) / Cos (T * Two_Pi);
949 function Tanh (A : Real) return Real is
951 if A < Half_Log_Epsilon then
954 elsif A > -Half_Log_Epsilon then
957 elsif abs A < Square_Root_Epsilon then
961 return Real (Aux.tanh (Double (A)));
964 ----------------------------
965 -- DEC-Specific functions --
966 ----------------------------
968 function LOG10 (A : REAL) return REAL is
970 return Log (A, 10.0);
973 function LOG2 (A : REAL) return REAL is
978 function ASIN (A : REAL) return REAL renames Arcsin;
979 function ACOS (A : REAL) return REAL renames Arccos;
981 function ATAN (A : REAL) return REAL is
983 return Arctan (A, 1.0);
986 function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
988 function SIND (A : REAL) return REAL is
990 return Sin (A, 360.0);
993 function COSD (A : REAL) return REAL is
995 return Cos (A, 360.0);
998 function TAND (A : REAL) return REAL is
1000 return Tan (A, 360.0);
1003 function ASIND (A : REAL) return REAL is
1005 return Arcsin (A, 360.0);
1008 function ACOSD (A : REAL) return REAL is
1010 return Arccos (A, 360.0);
1013 function Arctan (A : REAL) return REAL is
1015 return Arctan (A, 1.0, 360.0);
1018 function ATAND (A : REAL) return REAL is
1020 return Arctan (A, 1.0, 360.0);
1023 function ATAN2D (A1, A2 : REAL) return REAL is
1025 return Arctan (A1, A2, 360.0);