1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
9 -- Copyright (C) 1992-2005 Free Software Foundation, Inc. --
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. --
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. --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 ------------------------------------------------------------------------------
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.
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.
43 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
44 -- sinh, cosh, tanh from C library via math.h
46 -- This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
47 -- provides a compatible body for the DEC Math_Lib package.
49 with Ada
.Numerics
.Aux
;
50 use type Ada
.Numerics
.Aux
.Double
;
51 with Ada
.Numerics
; use Ada
.Numerics
;
53 package body Math_Lib
is
55 Log_Two
: constant := 0.69314_71805_59945_30941_72321_21458_17656_80755
;
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;
63 subtype Double is Aux.Double;
65 DEpsilon : constant Double := Double (Epsilon);
66 DIEpsilon : constant Double := Double (IEpsilon);
68 -----------------------
69 -- Local Subprograms --
70 -----------------------
83 function Exact_Remainder
87 -- Computes exact remainder of A divided by Y
89 function Half_Log_Epsilon return Real;
90 -- Function to provide constant: 0.5 * Log (Epsilon)
96 -- Common code for arc tangent after cyele reduction
98 function Log_Inverse_Epsilon return Real;
99 -- Function to provide constant: Log (1.0 / Epsilon)
101 function Square_Root_Epsilon return Real;
102 -- Function to provide constant: Sqrt (Epsilon)
108 function "**" (A1, A2 : Real) return Real is
114 raise Argument_Error;
117 raise Argument_Error;
124 raise Constraint_Error;
141 Real (Aux.pow (Double (A1), Double (A2)));
146 raise Constraint_Error;
157 function Arccos (A : Real) return Real is
162 raise Argument_Error;
164 elsif abs A < Square_Root_Epsilon then
174 Temp := Real (Aux.acos (Double (A)));
185 function Arccos (A, Cycle : Real) return Real is
190 raise Argument_Error;
192 elsif abs A > 1.0 then
193 raise Argument_Error;
195 elsif abs A < Square_Root_Epsilon then
205 Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
208 Temp := Cycle / 2.0 + Temp;
218 function Arccosh (A : Real) return Real is
220 -- Return Log (A - Sqrt (A * A - 1.0)); double valued,
221 -- only positive value returned
222 -- What is this comment ???
225 raise Argument_Error;
227 elsif A < 1.0 + Square_Root_Epsilon then
230 elsif abs A > 1.0 / Square_Root_Epsilon then
231 return Log (A) + Log_Two;
234 return Log (A + Sqrt (A * A - 1.0));
250 -- Just reverse arguments
252 return Arctan (Y, A);
264 -- Just reverse arguments
266 return Arctan (Y, A, Cycle);
273 function Arccoth (A : Real) return Real is
276 raise Constraint_Error;
278 elsif abs A < 1.0 then
279 raise Argument_Error;
281 elsif abs A > 1.0 / Epsilon then
285 return 0.5 * Log ((1.0 + A) / (A - 1.0));
295 function Arcsin (A : Real) return Real is
298 raise Argument_Error;
300 elsif abs A < Square_Root_Epsilon then
310 return Real (Aux.asin (Double (A)));
315 function Arcsin (A, Cycle : Real) return Real is
318 raise Argument_Error;
320 elsif abs A > 1.0 then
321 raise Argument_Error;
333 return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
340 function Arcsinh (A : Real) return Real is
342 if abs A < Square_Root_Epsilon then
345 elsif A > 1.0 / Square_Root_Epsilon then
346 return Log (A) + Log_Two;
348 elsif A < -1.0 / Square_Root_Epsilon then
349 return -(Log (-A) + Log_Two);
352 return -Log (abs A + Sqrt (A * A + 1.0));
355 return Log (A + Sqrt (A * A + 1.0));
374 raise Argument_Error;
391 return Local_Atan (Y, A);
405 raise Argument_Error;
410 raise Argument_Error;
427 return Local_Atan (Y, A) * Cycle / Two_Pi;
435 function Arctanh (A : Real) return Real is
438 raise Constraint_Error;
440 elsif abs A > 1.0 then
441 raise Argument_Error;
443 elsif abs A < Square_Root_Epsilon then
447 return 0.5 * Log ((1.0 + A) / (1.0 - A));
457 function Cos (A : Real) return Real is
462 elsif abs A < Square_Root_Epsilon then
467 return Real (Aux.Cos (Double (A)));
472 function Cos (A, Cycle : Real) return Real is
477 raise Argument_Error;
483 T := Exact_Remainder (abs (A), Cycle) / Cycle;
492 elsif T = 0.5 or T = -0.5 then
496 return Real (Aux.Cos (Double (T * Two_Pi)));
503 function Cosh (A : Real) return Real is
505 if abs A < Square_Root_Epsilon then
508 elsif abs A > Log_Inverse_Epsilon then
509 return Exp ((abs A) - Log_Two);
512 return Real (Aux.cosh (Double (A)));
516 raise Constraint_Error;
525 function Cot (A : Real) return Real is
528 raise Constraint_Error;
530 elsif abs A < Square_Root_Epsilon then
534 return Real (1.0 / Real'Base (Aux.tan (Double (A))));
539 function Cot (A, Cycle : Real) return Real is
544 raise Argument_Error;
547 raise Constraint_Error;
549 elsif abs A < Square_Root_Epsilon then
553 T := Exact_Remainder (A, Cycle) / Cycle;
555 if T = 0.0 or T = 0.5 or T = -0.5 then
556 raise Constraint_Error;
558 return Cos (T * Two_Pi) / Sin (T * Two_Pi);
566 function Coth (A : Real) return Real is
569 raise Constraint_Error;
571 elsif A < Half_Log_Epsilon then
574 elsif A > -Half_Log_Epsilon then
577 elsif abs A < Square_Root_Epsilon then
581 return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
584 ---------------------
585 -- Exact_Remainder --
586 ---------------------
588 function Exact_Remainder
593 Denominator : Real'Base := abs A;
594 Divisor : Real'Base := abs Y;
596 Sign : Real'Base := 1.0;
600 raise Constraint_Error;
608 elsif Denominator < Divisor then
612 while Denominator >= Divisor loop
614 -- Put divisors mantissa with denominators exponent to make reducer
619 while Reducer * 1_048_576.0 < Denominator loop
620 Reducer := Reducer * 1_048_576.0;
628 while Reducer * 1_024.0 < Denominator loop
629 Reducer := Reducer * 1_024.0;
637 while Reducer * 2.0 < Denominator loop
638 Reducer := Reducer * 2.0;
645 Denominator := Denominator - Reducer;
659 function Exp (A : Real) return Real is
667 Result := Real (Aux.Exp (Double (A)));
669 -- The check here catches the case of Exp returning IEEE infinity
671 if Result > Real'Last then
672 raise Constraint_Error;
679 ----------------------
680 -- Half_Log_Epsilon --
681 ----------------------
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!
686 function Half_Log_Epsilon return Real is
688 return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
689 end Half_Log_Epsilon;
701 Raw_Atan : Real'Base;
704 if abs Y > abs A then
710 if Z < Square_Root_Epsilon then
714 Raw_Atan := Pi / 4.0;
716 elsif Z < Square_Root_Epsilon then
720 Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
723 if abs Y > abs A then
724 Raw_Atan := Half_Pi - Raw_Atan;
736 return Pi - Raw_Atan;
738 return -(Pi - Raw_Atan);
749 function Log (A : Real) return Real is
752 raise Argument_Error;
755 raise Constraint_Error;
761 return Real (Aux.Log (Double (A)));
766 function Log (A, Base : Real) return Real is
769 raise Argument_Error;
771 elsif Base <= 0.0 or else Base = 1.0 then
772 raise Argument_Error;
775 raise Constraint_Error;
781 return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
784 -------------------------
785 -- Log_Inverse_Epsilon --
786 -------------------------
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!
791 function Log_Inverse_Epsilon return Real is
793 return Real (Aux.Log (DIEpsilon));
794 end Log_Inverse_Epsilon;
802 function Sin (A : Real) return Real is
804 if abs A < Square_Root_Epsilon then
808 return Real (Aux.Sin (Double (A)));
813 function Sin (A, Cycle : Real) return Real is
818 raise Argument_Error;
824 T := Exact_Remainder (A, Cycle) / Cycle;
826 if T = 0.0 or T = 0.5 or T = -0.5 then
829 elsif T = 0.25 or T = -0.75 then
832 elsif T = -0.25 or T = 0.75 then
837 return Real (Aux.Sin (Double (T * Two_Pi)));
844 function Sinh (A : Real) return Real is
846 if abs A < Square_Root_Epsilon then
849 elsif A > Log_Inverse_Epsilon then
850 return Exp (A - Log_Two);
852 elsif A < -Log_Inverse_Epsilon then
853 return -Exp ((-A) - Log_Two);
856 return Real (Aux.Sinh (Double (A)));
860 raise Constraint_Error;
863 -------------------------
864 -- Square_Root_Epsilon --
865 -------------------------
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!
870 function Square_Root_Epsilon return Real is
872 return Real (Aux.Sqrt (DEpsilon));
873 end Square_Root_Epsilon;
879 function Sqrt (A : Real) return Real is
882 raise Argument_Error;
884 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
889 -- Sqrt (1.0) must be exact for good complex accuracy
896 return Real (Aux.Sqrt (Double (A)));
905 function Tan (A : Real) return Real is
907 if abs A < Square_Root_Epsilon then
910 elsif abs A = Pi / 2.0 then
911 raise Constraint_Error;
914 return Real (Aux.tan (Double (A)));
919 function Tan (A, Cycle : Real) return Real is
924 raise Argument_Error;
930 T := Exact_Remainder (A, Cycle) / Cycle;
937 raise Constraint_Error;
940 return Sin (T * Two_Pi) / Cos (T * Two_Pi);
948 function Tanh (A : Real) return Real is
950 if A < Half_Log_Epsilon then
953 elsif A > -Half_Log_Epsilon then
956 elsif abs A < Square_Root_Epsilon then
960 return Real (Aux.tanh (Double (A)));
963 ----------------------------
964 -- DEC-Specific functions --
965 ----------------------------
967 function LOG10 (A : REAL) return REAL is
969 return Log (A, 10.0);
972 function LOG2 (A : REAL) return REAL is
977 function ASIN (A : REAL) return REAL renames Arcsin;
978 function ACOS (A : REAL) return REAL renames Arccos;
980 function ATAN (A : REAL) return REAL is
982 return Arctan (A, 1.0);
985 function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
987 function SIND (A : REAL) return REAL is
989 return Sin (A, 360.0);
992 function COSD (A : REAL) return REAL is
994 return Cos (A, 360.0);
997 function TAND (A : REAL) return REAL is
999 return Tan (A, 360.0);
1002 function ASIND (A : REAL) return REAL is
1004 return Arcsin (A, 360.0);
1007 function ACOSD (A : REAL) return REAL is
1009 return Arccos (A, 360.0);
1012 function Arctan (A : REAL) return REAL is
1014 return Arctan (A, 1.0, 360.0);
1017 function ATAND (A : REAL) return REAL is
1019 return Arctan (A, 1.0, 360.0);
1022 function ATAN2D (A1, A2 : REAL) return REAL is
1024 return Arctan (A1, A2, 360.0);