1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS --
9 -- Copyright (C) 1992-2009, 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 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. --
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. --
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/>. --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 ------------------------------------------------------------------------------
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.
38 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, sinh,
39 -- cosh, tanh from C library via math.h
41 with Ada.Numerics.Aux;
43 package body Ada.Numerics.Generic_Elementary_Functions is
45 use type Ada.Numerics.Aux.Double;
47 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
48 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
50 Half_Log_Two : constant := Log_Two / 2;
52 subtype T is Float_Type'Base;
53 subtype Double is Aux.Double;
55 Two_Pi : constant T := 2.0 * Pi;
56 Half_Pi : constant T := Pi / 2.0;
58 Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two;
59 Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two;
60 Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
62 -----------------------
63 -- Local Subprograms --
64 -----------------------
66 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base;
67 -- Cody/Waite routine, supposedly more precise than the library version.
68 -- Currently only needed for Sinh/Cosh on X86 with the largest FP type.
72 X : Float_Type'Base := 1.0) return Float_Type'Base;
73 -- Common code for arc tangent after cycle reduction
79 function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
80 A_Right : Float_Type'Base;
82 Result : Float_Type'Base;
84 Rest : Float_Type'Base;
95 elsif Right = 0.0 then
100 raise Constraint_Error;
105 elsif Left = 1.0 then
108 elsif Right = 1.0 then
116 elsif Right = 0.5 then
120 A_Right := abs (Right);
122 -- If exponent is larger than one, compute integer exponen-
123 -- tiation if possible, and evaluate fractional part with more
124 -- precision. The relative error is now proportional to the
125 -- fractional part of the exponent only.
128 and then A_Right < Float_Type'Base (Integer'Last)
130 Int_Part := Integer (Float_Type'Base'Truncation (A_Right));
131 Result := Left ** Int_Part;
132 Rest := A_Right - Float_Type'Base (Int_Part);
134 -- Compute with two leading bits of the mantissa using
135 -- square roots. Bound to be better than logarithms, and
136 -- easily extended to greater precision.
140 Result := Result * R1;
144 Result := Result * Sqrt (R1);
148 elsif Rest >= 0.25 then
149 Result := Result * Sqrt (Sqrt (Left));
154 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));
159 return (1.0 / Result);
163 Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
169 raise Constraint_Error;
180 function Arccos (X : Float_Type'Base) return Float_Type'Base is
181 Temp : Float_Type'Base;
185 raise Argument_Error;
187 elsif abs X < Sqrt_Epsilon then
197 Temp := Float_Type'Base (Aux.Acos (Double (X)));
208 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
209 Temp : Float_Type'Base;
213 raise Argument_Error;
215 elsif abs X > 1.0 then
216 raise Argument_Error;
218 elsif abs X < Sqrt_Epsilon then
228 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);
231 Temp := Cycle / 2.0 + Temp;
241 function Arccosh (X : Float_Type'Base) return Float_Type'Base is
243 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or the proper
244 -- approximation for X close to 1 or >> 1.
247 raise Argument_Error;
249 elsif X < 1.0 + Sqrt_Epsilon then
250 return Sqrt (2.0 * (X - 1.0));
252 elsif X > 1.0 / Sqrt_Epsilon then
253 return Log (X) + Log_Two;
256 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
267 (X : Float_Type'Base;
268 Y : Float_Type'Base := 1.0)
269 return Float_Type'Base
272 -- Just reverse arguments
274 return Arctan (Y, X);
280 (X : Float_Type'Base;
281 Y : Float_Type'Base := 1.0;
282 Cycle : Float_Type'Base)
283 return Float_Type'Base
286 -- Just reverse arguments
288 return Arctan (Y, X, Cycle);
295 function Arccoth (X : Float_Type'Base) return Float_Type'Base is
298 return Arctanh (1.0 / X);
300 elsif abs X = 1.0 then
301 raise Constraint_Error;
303 elsif abs X < 1.0 then
304 raise Argument_Error;
307 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the other
308 -- has error 0 or Epsilon.
310 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
320 function Arcsin (X : Float_Type'Base) return Float_Type'Base is
323 raise Argument_Error;
325 elsif abs X < Sqrt_Epsilon then
335 return Float_Type'Base (Aux.Asin (Double (X)));
340 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
343 raise Argument_Error;
345 elsif abs X > 1.0 then
346 raise Argument_Error;
355 return -(Cycle / 4.0);
358 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
365 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
367 if abs X < Sqrt_Epsilon then
370 elsif X > 1.0 / Sqrt_Epsilon then
371 return Log (X) + Log_Two;
373 elsif X < -(1.0 / Sqrt_Epsilon) then
374 return -(Log (-X) + Log_Two);
377 return -Log (abs X + Sqrt (X * X + 1.0));
380 return Log (X + Sqrt (X * X + 1.0));
391 (Y : Float_Type'Base;
392 X : Float_Type'Base := 1.0)
393 return Float_Type'Base
396 if X = 0.0 and then Y = 0.0 then
397 raise Argument_Error;
403 return Pi * Float_Type'Copy_Sign (1.0, Y);
407 return Float_Type'Copy_Sign (Half_Pi, Y);
410 return Local_Atan (Y, X);
417 (Y : Float_Type'Base;
418 X : Float_Type'Base := 1.0;
419 Cycle : Float_Type'Base)
420 return Float_Type'Base
424 raise Argument_Error;
426 elsif X = 0.0 and then Y = 0.0 then
427 raise Argument_Error;
433 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
437 return Float_Type'Copy_Sign (Cycle / 4.0, Y);
440 return Local_Atan (Y, X) * Cycle / Two_Pi;
448 function Arctanh (X : Float_Type'Base) return Float_Type'Base is
449 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;
451 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa;
454 -- The naive formula:
456 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X)
458 -- is not well-behaved numerically when X < 0.5 and when X is close
459 -- to one. The following is accurate but probably not optimal.
462 raise Constraint_Error;
464 elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then
467 raise Argument_Error;
470 -- The one case that overflows if put through the method below:
471 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is
472 -- accurate. This simplifies to:
474 return Float_Type'Copy_Sign (
475 Half_Log_Two * Float_Type'Base (Mantissa + 1), X);
478 -- elsif abs X <= 0.5 then
479 -- why is above line commented out ???
482 -- Use several piecewise linear approximations. A is close to X,
483 -- chosen so 1.0 + A, 1.0 - A, and X - A are exact. The two scalings
484 -- remove the low-order bits of X.
486 A := Float_Type'Base'Scaling (
487 Float_Type'Base (Long_Long_Integer
488 (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa);
490 B := X - A; -- This is exact; abs B <= 2**(-Mantissa).
491 A_Plus_1 := 1.0 + A; -- This is exact.
492 A_From_1 := 1.0 - A; -- Ditto.
493 D := A_Plus_1 * A_From_1; -- 1 - A*A.
495 -- use one term of the series expansion:
497 -- f (x + e) = f(x) + e * f'(x) + ..
499 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is
500 -- A*(B/D)**2 (if a quadratic approximation is ever needed).
502 return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D;
512 function Cos (X : Float_Type'Base) return Float_Type'Base is
517 elsif abs X < Sqrt_Epsilon then
522 return Float_Type'Base (Aux.Cos (Double (X)));
527 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
529 -- Just reuse the code for Sin. The potential small loss of speed is
530 -- negligible with proper (front-end) inlining.
532 return -Sin (abs X - Cycle * 0.25, Cycle);
539 function Cosh (X : Float_Type'Base) return Float_Type'Base is
540 Lnv : constant Float_Type'Base := 8#0.542714#;
541 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
542 Y : constant Float_Type'Base := abs X;
546 if Y < Sqrt_Epsilon then
549 elsif Y > Log_Inverse_Epsilon then
550 Z := Exp_Strict (Y - Lnv);
551 return (Z + V2minus1 * Z);
555 return 0.5 * (Z + 1.0 / Z);
566 function Cot (X : Float_Type'Base) return Float_Type'Base is
569 raise Constraint_Error;
571 elsif abs X < Sqrt_Epsilon then
575 return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
580 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
585 raise Argument_Error;
588 T := Float_Type'Base'Remainder (X, Cycle);
590 if T = 0.0 or else abs T = 0.5 * Cycle then
591 raise Constraint_Error;
593 elsif abs T < Sqrt_Epsilon then
596 elsif abs T = 0.25 * Cycle then
600 T := T / Cycle * Two_Pi;
601 return Cos (T) / Sin (T);
609 function Coth (X : Float_Type'Base) return Float_Type'Base is
612 raise Constraint_Error;
614 elsif X < Half_Log_Epsilon then
617 elsif X > -Half_Log_Epsilon then
620 elsif abs X < Sqrt_Epsilon then
624 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
631 function Exp (X : Float_Type'Base) return Float_Type'Base is
632 Result : Float_Type'Base;
639 Result := Float_Type'Base (Aux.Exp (Double (X)));
641 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
642 -- is False, then we can just leave it as an infinity (and indeed we
643 -- prefer to do so). But if Machine_Overflows is True, then we have
644 -- to raise a Constraint_Error exception as required by the RM.
646 if Float_Type'Machine_Overflows and then not Result'Valid then
647 raise Constraint_Error;
657 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
661 P0 : constant := 0.25000_00000_00000_00000;
662 P1 : constant := 0.75753_18015_94227_76666E-2;
663 P2 : constant := 0.31555_19276_56846_46356E-4;
665 Q0 : constant := 0.5;
666 Q1 : constant := 0.56817_30269_85512_21787E-1;
667 Q2 : constant := 0.63121_89437_43985_02557E-3;
668 Q3 : constant := 0.75104_02839_98700_46114E-6;
670 C1 : constant := 8#0.543#;
671 C2 : constant := -2.1219_44400_54690_58277E-4;
672 Le : constant := 1.4426_95040_88896_34074;
674 XN : Float_Type'Base;
675 P, Q, R : Float_Type'Base;
682 XN := Float_Type'Base'Rounding (X * Le);
683 G := (X - XN * C1) - XN * C2;
685 P := G * ((P2 * Z + P1) * Z + P0);
686 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
687 R := 0.5 + P / (Q - P);
689 R := Float_Type'Base'Scaling (R, Integer (XN) + 1);
691 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
692 -- is False, then we can just leave it as an infinity (and indeed we
693 -- prefer to do so). But if Machine_Overflows is True, then we have to
694 -- raise a Constraint_Error exception as required by the RM.
696 if Float_Type'Machine_Overflows and then not R'Valid then
697 raise Constraint_Error;
709 (Y : Float_Type'Base;
710 X : Float_Type'Base := 1.0) return Float_Type'Base
713 Raw_Atan : Float_Type'Base;
716 Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));
719 (if Z < Sqrt_Epsilon then Z
720 elsif Z = 1.0 then Pi / 4.0
721 else Float_Type'Base (Aux.Atan (Double (Z))));
723 if abs Y > abs X then
724 Raw_Atan := Half_Pi - Raw_Atan;
728 return Float_Type'Copy_Sign (Raw_Atan, Y);
730 return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
740 function Log (X : Float_Type'Base) return Float_Type'Base is
743 raise Argument_Error;
746 raise Constraint_Error;
752 return Float_Type'Base (Aux.Log (Double (X)));
757 function Log (X, Base : Float_Type'Base) return Float_Type'Base is
760 raise Argument_Error;
762 elsif Base <= 0.0 or else Base = 1.0 then
763 raise Argument_Error;
766 raise Constraint_Error;
772 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
781 function Sin (X : Float_Type'Base) return Float_Type'Base is
783 if abs X < Sqrt_Epsilon then
787 return Float_Type'Base (Aux.Sin (Double (X)));
792 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
797 raise Argument_Error;
799 -- If X is zero, return it as the result, preserving the argument sign.
800 -- Is this test really needed on any machine ???
806 T := Float_Type'Base'Remainder (X, Cycle);
808 -- The following two reductions reduce the argument to the interval
809 -- [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed
810 -- to prevent inaccuracy that may result if the sinus function uses a
811 -- different (more accurate) value of Pi in its reduction than is used
812 -- in the multiplication with Two_Pi.
814 if abs T > 0.25 * Cycle then
815 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
818 -- Could test for 12.0 * abs T = Cycle, and return an exact value in
819 -- those cases. It is not clear this is worth the extra test though.
821 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
828 function Sinh (X : Float_Type'Base) return Float_Type'Base is
829 Lnv : constant Float_Type'Base := 8#0.542714#;
830 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
831 Y : constant Float_Type'Base := abs X;
832 F : constant Float_Type'Base := Y * Y;
835 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;
838 if Y < Sqrt_Epsilon then
841 elsif Y > Log_Inverse_Epsilon then
842 Z := Exp_Strict (Y - Lnv);
843 Z := Z + V2minus1 * Z;
847 if Float_Digits_1_6 then
849 -- Use expansion provided by Cody and Waite, p. 226. Note that
850 -- leading term of the polynomial in Q is exactly 1.0.
853 P0 : constant := -0.71379_3159E+1;
854 P1 : constant := -0.19033_3399E+0;
855 Q0 : constant := -0.42827_7109E+2;
858 Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
863 P0 : constant := -0.35181_28343_01771_17881E+6;
864 P1 : constant := -0.11563_52119_68517_68270E+5;
865 P2 : constant := -0.16375_79820_26307_51372E+3;
866 P3 : constant := -0.78966_12741_73570_99479E+0;
867 Q0 : constant := -0.21108_77005_81062_71242E+7;
868 Q1 : constant := 0.36162_72310_94218_36460E+5;
869 Q2 : constant := -0.27773_52311_96507_01667E+3;
872 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
873 / (((F + Q2) * F + Q1) * F + Q0);
879 Z := 0.5 * (Z - 1.0 / Z);
893 function Sqrt (X : Float_Type'Base) return Float_Type'Base is
896 raise Argument_Error;
898 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
904 return Float_Type'Base (Aux.Sqrt (Double (X)));
913 function Tan (X : Float_Type'Base) return Float_Type'Base is
915 if abs X < Sqrt_Epsilon then
918 elsif abs X = Pi / 2.0 then
919 raise Constraint_Error;
922 return Float_Type'Base (Aux.Tan (Double (X)));
927 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
932 raise Argument_Error;
938 T := Float_Type'Base'Remainder (X, Cycle);
940 if abs T = 0.25 * Cycle then
941 raise Constraint_Error;
943 elsif abs T = 0.5 * Cycle then
947 T := T / Cycle * Two_Pi;
948 return Sin (T) / Cos (T);
957 function Tanh (X : Float_Type'Base) return Float_Type'Base is
958 P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
959 P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
960 P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;
962 Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4;
963 Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4;
964 Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3;
965 Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1;
967 Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
969 P, Q, R : Float_Type'Base;
970 Y : constant Float_Type'Base := abs X;
971 G : constant Float_Type'Base := Y * Y;
973 Float_Type_Digits_15_Or_More : constant Boolean :=
974 Float_Type'Digits > 14;
977 if X < Half_Log_Epsilon then
980 elsif X > -Half_Log_Epsilon then
983 elsif Y < Sqrt_Epsilon then
987 and then Float_Type_Digits_15_Or_More
989 P := (P2 * G + P1) * G + P0;
990 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
995 return Float_Type'Base (Aux.Tanh (Double (X)));
999 end Ada.Numerics.Generic_Elementary_Functions;