1 ------------------------------------------------------------------------------
3 -- GNAT RUNTIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS --
10 -- Copyright (C) 1992-2001, 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 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
42 -- sinh, cosh, tanh from C library via math.h
44 with Ada
.Numerics
.Aux
;
46 package body Ada
.Numerics
.Generic_Elementary_Functions
is
48 use type Ada
.Numerics
.Aux
.Double
;
50 Sqrt_Two
: constant := 1.41421_35623_73095_04880_16887_24209_69807_85696
;
51 Log_Two
: constant := 0.69314_71805_59945_30941_72321_21458_17656_80755
;
52 Half_Log_Two
: constant := Log_Two
/ 2;
54 subtype T
is Float_Type
'Base;
55 subtype Double
is Aux
.Double
;
57 Two_Pi
: constant T
:= 2.0 * Pi
;
58 Half_Pi
: constant T
:= Pi
/ 2.0;
59 Fourth_Pi
: constant T
:= Pi
/ 4.0;
61 Epsilon
: constant T
:= 2.0 ** (1 - T
'Model_Mantissa);
62 IEpsilon
: constant T
:= 2.0 ** (T
'Model_Mantissa - 1);
63 Log_Epsilon
: constant T
:= T
(1 - T
'Model_Mantissa) * Log_Two
;
64 Half_Log_Epsilon
: constant T
:= T
(1 - T
'Model_Mantissa) * Half_Log_Two
;
65 Log_Inverse_Epsilon
: constant T
:= T
(T
'Model_Mantissa - 1) * Log_Two
;
66 Sqrt_Epsilon
: constant T
:= Sqrt_Two
** (1 - T
'Model_Mantissa);
68 DEpsilon
: constant Double
:= Double
(Epsilon
);
69 DIEpsilon
: constant Double
:= Double
(IEpsilon
);
71 -----------------------
72 -- Local Subprograms --
73 -----------------------
75 function Exp_Strict
(X
: Float_Type
'Base) return Float_Type
'Base;
76 -- Cody/Waite routine, supposedly more precise than the library
77 -- version. Currently only needed for Sinh/Cosh on X86 with the largest
82 X
: Float_Type
'Base := 1.0)
83 return Float_Type
'Base;
84 -- Common code for arc tangent after cyele reduction
90 function "**" (Left
, Right
: Float_Type
'Base) return Float_Type
'Base is
91 A_Right
: Float_Type
'Base;
93 Result
: Float_Type
'Base;
95 Rest
: Float_Type
'Base;
101 raise Argument_Error
;
103 elsif Left
< 0.0 then
104 raise Argument_Error
;
106 elsif Right
= 0.0 then
109 elsif Left
= 0.0 then
111 raise Constraint_Error
;
116 elsif Left
= 1.0 then
119 elsif Right
= 1.0 then
127 elsif Right
= 0.5 then
131 A_Right
:= abs (Right
);
133 -- If exponent is larger than one, compute integer exponen-
134 -- tiation if possible, and evaluate fractional part with
135 -- more precision. The relative error is now proportional
136 -- to the fractional part of the exponent only.
139 and then A_Right
< Float_Type
'Base (Integer'Last)
141 Int_Part
:= Integer (Float_Type
'Base'Truncation (A_Right));
142 Result := Left ** Int_Part;
143 Rest := A_Right - Float_Type'Base (Int_Part);
145 -- Compute with two leading bits of the mantissa using
146 -- square roots. Bound to be better than logarithms, and
147 -- easily extended to greater precision.
151 Result := Result * R1;
155 Result := Result * Sqrt (R1);
159 elsif Rest >= 0.25 then
160 Result := Result * Sqrt (Sqrt (Left));
165 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));
170 return (1.0 / Result);
174 Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
180 raise Constraint_Error;
191 function Arccos (X : Float_Type'Base) return Float_Type'Base is
192 Temp : Float_Type'Base;
196 raise Argument_Error;
198 elsif abs X < Sqrt_Epsilon then
208 Temp := Float_Type'Base (Aux.Acos (Double (X)));
219 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
220 Temp : Float_Type'Base;
224 raise Argument_Error;
226 elsif abs X > 1.0 then
227 raise Argument_Error;
229 elsif abs X < Sqrt_Epsilon then
239 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);
242 Temp := Cycle / 2.0 + Temp;
252 function Arccosh (X : Float_Type'Base) return Float_Type'Base is
254 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or
255 -- the proper approximation for X close to 1 or >> 1.
258 raise Argument_Error;
260 elsif X < 1.0 + Sqrt_Epsilon then
261 return Sqrt (2.0 * (X - 1.0));
263 elsif X > 1.0 / Sqrt_Epsilon then
264 return Log (X) + Log_Two;
267 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
278 (X : Float_Type'Base;
279 Y : Float_Type'Base := 1.0)
280 return Float_Type'Base
283 -- Just reverse arguments
285 return Arctan (Y, X);
291 (X : Float_Type'Base;
292 Y : Float_Type'Base := 1.0;
293 Cycle : Float_Type'Base)
294 return Float_Type'Base
297 -- Just reverse arguments
299 return Arctan (Y, X, Cycle);
306 function Arccoth (X : Float_Type'Base) return Float_Type'Base is
309 return Arctanh (1.0 / X);
311 elsif abs X = 1.0 then
312 raise Constraint_Error;
314 elsif abs X < 1.0 then
315 raise Argument_Error;
318 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the
319 -- other has error 0 or Epsilon.
321 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
331 function Arcsin (X : Float_Type'Base) return Float_Type'Base is
334 raise Argument_Error;
336 elsif abs X < Sqrt_Epsilon then
346 return Float_Type'Base (Aux.Asin (Double (X)));
351 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
354 raise Argument_Error;
356 elsif abs X > 1.0 then
357 raise Argument_Error;
369 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
376 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
378 if abs X < Sqrt_Epsilon then
381 elsif X > 1.0 / Sqrt_Epsilon then
382 return Log (X) + Log_Two;
384 elsif X < -1.0 / Sqrt_Epsilon then
385 return -(Log (-X) + Log_Two);
388 return -Log (abs X + Sqrt (X * X + 1.0));
391 return Log (X + Sqrt (X * X + 1.0));
402 (Y : Float_Type'Base;
403 X : Float_Type'Base := 1.0)
404 return Float_Type'Base
410 raise Argument_Error;
416 return Pi * Float_Type'Copy_Sign (1.0, Y);
427 return Local_Atan (Y, X);
434 (Y : Float_Type'Base;
435 X : Float_Type'Base := 1.0;
436 Cycle : Float_Type'Base)
437 return Float_Type'Base
441 raise Argument_Error;
446 raise Argument_Error;
452 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
463 return Local_Atan (Y, X) * Cycle / Two_Pi;
471 function Arctanh (X : Float_Type'Base) return Float_Type'Base is
472 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;
473 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa
;
476 -- The naive formula:
478 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X)
480 -- is not well-behaved numerically when X < 0.5 and when X is close
481 -- to one. The following is accurate but probably not optimal.
484 raise Constraint_Error
;
486 elsif abs X
>= 1.0 - 2.0 ** (-Mantissa
) then
489 raise Argument_Error
;
492 -- The one case that overflows if put through the method below:
493 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is
494 -- accurate. This simplifies to:
496 return Float_Type
'Copy_Sign (
497 Half_Log_Two
* Float_Type
'Base (Mantissa
+ 1), X
);
500 -- elsif abs X <= 0.5 then
501 -- why is above line commented out ???
504 -- Use several piecewise linear approximations.
505 -- A is close to X, chosen so 1.0 + A, 1.0 - A, and X - A are exact.
506 -- The two scalings remove the low-order bits of X.
508 A
:= Float_Type
'Base'Scaling (
509 Float_Type'Base (Long_Long_Integer
510 (Float_Type'Base'Scaling
(X
, Mantissa
- 1))), 1 - Mantissa
);
512 B
:= X
- A
; -- This is exact; abs B <= 2**(-Mantissa).
513 A_Plus_1
:= 1.0 + A
; -- This is exact.
514 A_From_1
:= 1.0 - A
; -- Ditto.
515 D
:= A_Plus_1
* A_From_1
; -- 1 - A*A.
517 -- use one term of the series expansion:
518 -- f (x + e) = f(x) + e * f'(x) + ..
520 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is
521 -- A*(B/D)**2 (if a quadratic approximation is ever needed).
523 return 0.5 * (Log
(A_Plus_1
) - Log
(A_From_1
)) + B
/ D
;
526 -- return 0.5 * Log ((X + 1.0) / (1.0 - X));
527 -- why are above lines commented out ???
537 function Cos
(X
: Float_Type
'Base) return Float_Type
'Base is
542 elsif abs X
< Sqrt_Epsilon
then
547 return Float_Type
'Base (Aux
.Cos
(Double
(X
)));
552 function Cos
(X
, Cycle
: Float_Type
'Base) return Float_Type
'Base is
554 -- Just reuse the code for Sin. The potential small
555 -- loss of speed is negligible with proper (front-end) inlining.
557 return -Sin
(abs X
- Cycle
* 0.25, Cycle
);
564 function Cosh
(X
: Float_Type
'Base) return Float_Type
'Base is
565 Lnv
: constant Float_Type
'Base := 8#
0.542714#
;
566 V2minus1
: constant Float_Type
'Base := 0.13830_27787_96019_02638E
-4
;
567 Y
: Float_Type
'Base := abs X
;
571 if Y
< Sqrt_Epsilon
then
574 elsif Y
> Log_Inverse_Epsilon
then
575 Z
:= Exp_Strict
(Y
- Lnv
);
576 return (Z
+ V2minus1
* Z
);
580 return 0.5 * (Z
+ 1.0 / Z
);
591 function Cot
(X
: Float_Type
'Base) return Float_Type
'Base is
594 raise Constraint_Error
;
596 elsif abs X
< Sqrt_Epsilon
then
600 return 1.0 / Float_Type
'Base (Aux
.Tan
(Double
(X
)));
605 function Cot
(X
, Cycle
: Float_Type
'Base) return Float_Type
'Base is
610 raise Argument_Error
;
613 T
:= Float_Type
'Base'Remainder (X, Cycle);
615 if T = 0.0 or abs T = 0.5 * Cycle then
616 raise Constraint_Error;
618 elsif abs T < Sqrt_Epsilon then
621 elsif abs T = 0.25 * Cycle then
625 T := T / Cycle * Two_Pi;
626 return Cos (T) / Sin (T);
634 function Coth (X : Float_Type'Base) return Float_Type'Base is
637 raise Constraint_Error;
639 elsif X < Half_Log_Epsilon then
642 elsif X > -Half_Log_Epsilon then
645 elsif abs X < Sqrt_Epsilon then
649 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
656 function Exp (X : Float_Type'Base) return Float_Type'Base is
657 Result : Float_Type'Base;
664 Result := Float_Type'Base (Aux.Exp (Double (X)));
666 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
667 -- is False, then we can just leave it as an infinity (and indeed we
668 -- prefer to do so). But if Machine_Overflows is True, then we have
669 -- to raise a Constraint_Error exception as required by the RM.
671 if Float_Type'Machine_Overflows and then not Result'Valid then
672 raise Constraint_Error;
682 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
686 P0 : constant := 0.25000_00000_00000_00000;
687 P1 : constant := 0.75753_18015_94227_76666E-2;
688 P2 : constant := 0.31555_19276_56846_46356E-4;
690 Q0 : constant := 0.5;
691 Q1 : constant := 0.56817_30269_85512_21787E-1;
692 Q2 : constant := 0.63121_89437_43985_02557E-3;
693 Q3 : constant := 0.75104_02839_98700_46114E-6;
695 C1 : constant := 8#0.543#;
696 C2 : constant := -2.1219_44400_54690_58277E-4;
697 Le : constant := 1.4426_95040_88896_34074;
699 XN : Float_Type'Base;
700 P, Q, R : Float_Type'Base;
707 XN := Float_Type'Base'Rounding
(X
* Le
);
708 G
:= (X
- XN
* C1
) - XN
* C2
;
710 P
:= G
* ((P2
* Z
+ P1
) * Z
+ P0
);
711 Q
:= ((Q3
* Z
+ Q2
) * Z
+ Q1
) * Z
+ Q0
;
712 R
:= 0.5 + P
/ (Q
- P
);
714 R
:= Float_Type
'Base'Scaling (R, Integer (XN) + 1);
716 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
717 -- is False, then we can just leave it as an infinity (and indeed we
718 -- prefer to do so). But if Machine_Overflows is True, then we have
719 -- to raise a Constraint_Error exception as required by the RM.
721 if Float_Type'Machine_Overflows and then not R'Valid then
722 raise Constraint_Error;
734 (Y : Float_Type'Base;
735 X : Float_Type'Base := 1.0)
736 return Float_Type'Base
739 Raw_Atan : Float_Type'Base;
742 if abs Y > abs X then
748 if Z < Sqrt_Epsilon then
752 Raw_Atan := Pi / 4.0;
755 Raw_Atan := Float_Type'Base (Aux.Atan (Double (Z)));
758 if abs Y > abs X then
759 Raw_Atan := Half_Pi - Raw_Atan;
771 return Pi - Raw_Atan;
773 return -(Pi - Raw_Atan);
784 function Log (X : Float_Type'Base) return Float_Type'Base is
787 raise Argument_Error;
790 raise Constraint_Error;
796 return Float_Type'Base (Aux.Log (Double (X)));
801 function Log (X, Base : Float_Type'Base) return Float_Type'Base is
804 raise Argument_Error;
806 elsif Base <= 0.0 or else Base = 1.0 then
807 raise Argument_Error;
810 raise Constraint_Error;
816 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
825 function Sin (X : Float_Type'Base) return Float_Type'Base is
827 if abs X < Sqrt_Epsilon then
831 return Float_Type'Base (Aux.Sin (Double (X)));
836 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
841 raise Argument_Error;
844 -- Is this test really needed on any machine ???
848 T := Float_Type'Base'Remainder
(X
, Cycle
);
850 -- The following two reductions reduce the argument
851 -- to the interval [-0.25 * Cycle, 0.25 * Cycle].
852 -- This reduction is exact and is needed to prevent
853 -- inaccuracy that may result if the sinus function
854 -- a different (more accurate) value of Pi in its
855 -- reduction than is used in the multiplication with Two_Pi.
857 if abs T
> 0.25 * Cycle
then
858 T
:= 0.5 * Float_Type
'Copy_Sign (Cycle
, T
) - T
;
861 -- Could test for 12.0 * abs T = Cycle, and return
862 -- an exact value in those cases. It is not clear that
863 -- this is worth the extra test though.
865 return Float_Type
'Base (Aux
.Sin
(Double
(T
/ Cycle
* Two_Pi
)));
872 function Sinh
(X
: Float_Type
'Base) return Float_Type
'Base is
873 Lnv
: constant Float_Type
'Base := 8#
0.542714#
;
874 V2minus1
: constant Float_Type
'Base := 0.13830_27787_96019_02638E
-4
;
875 Y
: Float_Type
'Base := abs X
;
876 F
: constant Float_Type
'Base := Y
* Y
;
879 Float_Digits_1_6
: constant Boolean := Float_Type
'Digits < 7;
882 if Y
< Sqrt_Epsilon
then
885 elsif Y
> Log_Inverse_Epsilon
then
886 Z
:= Exp_Strict
(Y
- Lnv
);
887 Z
:= Z
+ V2minus1
* Z
;
891 if Float_Digits_1_6
then
893 -- Use expansion provided by Cody and Waite, p. 226. Note that
894 -- leading term of the polynomial in Q is exactly 1.0.
897 P0
: constant := -0.71379_3159E
+1
;
898 P1
: constant := -0.19033_3399E
+0
;
899 Q0
: constant := -0.42827_7109E
+2
;
902 Z
:= Y
+ Y
* F
* (P1
* F
+ P0
) / (F
+ Q0
);
907 P0
: constant := -0.35181_28343_01771_17881E
+6
;
908 P1
: constant := -0.11563_52119_68517_68270E
+5
;
909 P2
: constant := -0.16375_79820_26307_51372E
+3
;
910 P3
: constant := -0.78966_12741_73570_99479E
+0
;
911 Q0
: constant := -0.21108_77005_81062_71242E
+7
;
912 Q1
: constant := 0.36162_72310_94218_36460E
+5
;
913 Q2
: constant := -0.27773_52311_96507_01667E
+3
;
916 Z
:= Y
+ Y
* F
* (((P3
* F
+ P2
) * F
+ P1
) * F
+ P0
)
917 / (((F
+ Q2
) * F
+ Q1
) * F
+ Q0
);
923 Z
:= 0.5 * (Z
- 1.0 / Z
);
937 function Sqrt
(X
: Float_Type
'Base) return Float_Type
'Base is
940 raise Argument_Error
;
942 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
949 return Float_Type
'Base (Aux
.Sqrt
(Double
(X
)));
958 function Tan
(X
: Float_Type
'Base) return Float_Type
'Base is
960 if abs X
< Sqrt_Epsilon
then
963 elsif abs X
= Pi
/ 2.0 then
964 raise Constraint_Error
;
967 return Float_Type
'Base (Aux
.Tan
(Double
(X
)));
972 function Tan
(X
, Cycle
: Float_Type
'Base) return Float_Type
'Base is
977 raise Argument_Error
;
983 T
:= Float_Type
'Base'Remainder (X, Cycle);
985 if abs T = 0.25 * Cycle then
986 raise Constraint_Error;
988 elsif abs T = 0.5 * Cycle then
992 T := T / Cycle * Two_Pi;
993 return Sin (T) / Cos (T);
1002 function Tanh (X : Float_Type'Base) return Float_Type'Base is
1003 P0 : constant Float_Type'Base := -0.16134_11902E4;
1004 P1 : constant Float_Type'Base := -0.99225_92967E2;
1005 P2 : constant Float_Type'Base := -0.96437_49299E0;
1007 Q0 : constant Float_Type'Base := 0.48402_35707E4;
1008 Q1 : constant Float_Type'Base := 0.22337_72071E4;
1009 Q2 : constant Float_Type'Base := 0.11274_47438E3;
1010 Q3 : constant Float_Type'Base := 0.10000000000E1;
1012 Half_Ln3 : constant Float_Type'Base := 0.54930_61443;
1014 P, Q, R : Float_Type'Base;
1015 Y : Float_Type'Base := abs X;
1016 G : Float_Type'Base := Y * Y;
1018 Float_Type_Digits_15_Or_More : constant Boolean :=
1019 Float_Type'Digits > 14;
1022 if X < Half_Log_Epsilon then
1025 elsif X > -Half_Log_Epsilon then
1028 elsif Y < Sqrt_Epsilon then
1032 and then Float_Type_Digits_15_Or_More
1034 P := (P2 * G + P1) * G + P0;
1035 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
1040 return Float_Type'Base (Aux.Tanh (Double (X)));
1044 end Ada.Numerics.Generic_Elementary_Functions;