| blob | f6d01d8e017b77a9c1d1b94b6a607df886c2f43a |
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUNTIME COMPONENTS --
4 -- --
5 -- M A T H _ L I B --
6 -- --
7 -- B o d y --
8 -- --
9 -- --
10 -- Copyright (C) 1992-2000 Free Software Foundation, Inc. --
11 -- --
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. --
22 -- --
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. --
29 -- --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- Extensive contributions were provided by Ada Core Technologies Inc. --
32 -- --
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.
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 -----------------------
73 function Arctan
74 (Y : Real;
75 A : Real := 1.0)
76 return Real;
78 function Arctan
79 (Y : Real;
80 A : Real := 1.0;
81 Cycle : Real)
82 return Real;
84 function Exact_Remainder
85 (A : Real;
86 Y : Real)
87 return Real;
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)
93 function Local_Atan
94 (Y : Real;
95 A : Real := 1.0)
96 return Real;
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)
105 ----------
106 -- "**" --
107 ----------
109 function "**" (A1, A2 : Real) return Real is
111 begin
112 if A1 = 0.0
113 and then A2 = 0.0
114 then
115 raise Argument_Error;
117 elsif A1 < 0.0 then
118 raise Argument_Error;
120 elsif A2 = 0.0 then
121 return 1.0;
123 elsif A1 = 0.0 then
124 if A2 < 0.0 then
125 raise Constraint_Error;
126 else
127 return 0.0;
128 end if;
130 elsif A1 = 1.0 then
131 return 1.0;
133 elsif A2 = 1.0 then
134 return A1;
136 else
137 begin
138 if A2 = 2.0 then
139 return A1 * A1;
140 else
141 return
142 Real (Aux.pow (Double (A1), Double (A2)));
143 end if;
145 exception
146 when others =>
147 raise Constraint_Error;
148 end;
149 end if;
150 end "**";
152 ------------
153 -- Arccos --
154 ------------
156 -- Natural cycle
158 function Arccos (A : Real) return Real is
159 Temp : Real'Base;
161 begin
162 if abs A > 1.0 then
163 raise Argument_Error;
165 elsif abs A < Square_Root_Epsilon then
166 return Pi / 2.0 - A;
168 elsif A = 1.0 then
169 return 0.0;
171 elsif A = -1.0 then
172 return Pi;
173 end if;
175 Temp := Real (Aux.acos (Double (A)));
177 if Temp < 0.0 then
178 Temp := Pi + Temp;
179 end if;
181 return Temp;
182 end Arccos;
184 -- Arbitrary cycle
186 function Arccos (A, Cycle : Real) return Real is
187 Temp : Real'Base;
189 begin
190 if Cycle <= 0.0 then
191 raise Argument_Error;
193 elsif abs A > 1.0 then
194 raise Argument_Error;
196 elsif abs A < Square_Root_Epsilon then
197 return Cycle / 4.0;
199 elsif A = 1.0 then
200 return 0.0;
202 elsif A = -1.0 then
203 return Cycle / 2.0;
204 end if;
206 Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
208 if Temp < 0.0 then
209 Temp := Cycle / 2.0 + Temp;
210 end if;
212 return Temp;
213 end Arccos;
215 -------------
216 -- Arccosh --
217 -------------
219 function Arccosh (A : Real) return Real is
220 begin
221 -- Return Log (A - Sqrt (A * A - 1.0)); double valued,
222 -- only positive value returned
223 -- What is this comment ???
225 if A < 1.0 then
226 raise Argument_Error;
228 elsif A < 1.0 + Square_Root_Epsilon then
229 return A - 1.0;
231 elsif abs A > 1.0 / Square_Root_Epsilon then
232 return Log (A) + Log_Two;
234 else
235 return Log (A + Sqrt (A * A - 1.0));
236 end if;
237 end Arccosh;
239 ------------
240 -- Arccot --
241 ------------
243 -- Natural cycle
245 function Arccot
246 (A : Real;
247 Y : Real := 1.0)
248 return Real
249 is
250 begin
251 -- Just reverse arguments
253 return Arctan (Y, A);
254 end Arccot;
256 -- Arbitrary cycle
258 function Arccot
259 (A : Real;
260 Y : Real := 1.0;
261 Cycle : Real)
262 return Real
263 is
264 begin
265 -- Just reverse arguments
267 return Arctan (Y, A, Cycle);
268 end Arccot;
270 -------------
271 -- Arccoth --
272 -------------
274 function Arccoth (A : Real) return Real is
275 begin
276 if abs A = 1.0 then
277 raise Constraint_Error;
279 elsif abs A < 1.0 then
280 raise Argument_Error;
282 elsif abs A > 1.0 / Epsilon then
283 return 0.0;
285 else
286 return 0.5 * Log ((1.0 + A) / (A - 1.0));
287 end if;
288 end Arccoth;
290 ------------
291 -- Arcsin --
292 ------------
294 -- Natural cycle
296 function Arcsin (A : Real) return Real is
297 begin
298 if abs A > 1.0 then
299 raise Argument_Error;
301 elsif abs A < Square_Root_Epsilon then
302 return A;
304 elsif A = 1.0 then
305 return Pi / 2.0;
307 elsif A = -1.0 then
308 return -Pi / 2.0;
309 end if;
311 return Real (Aux.asin (Double (A)));
312 end Arcsin;
314 -- Arbitrary cycle
316 function Arcsin (A, Cycle : Real) return Real is
317 begin
318 if Cycle <= 0.0 then
319 raise Argument_Error;
321 elsif abs A > 1.0 then
322 raise Argument_Error;
324 elsif A = 0.0 then
325 return A;
327 elsif A = 1.0 then
328 return Cycle / 4.0;
330 elsif A = -1.0 then
331 return -Cycle / 4.0;
332 end if;
334 return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
335 end Arcsin;
337 -------------
338 -- Arcsinh --
339 -------------
341 function Arcsinh (A : Real) return Real is
342 begin
343 if abs A < Square_Root_Epsilon then
344 return A;
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);
352 elsif A < 0.0 then
353 return -Log (abs A + Sqrt (A * A + 1.0));
355 else
356 return Log (A + Sqrt (A * A + 1.0));
357 end if;
358 end Arcsinh;
360 ------------
361 -- Arctan --
362 ------------
364 -- Natural cycle
366 function Arctan
367 (Y : Real;
368 A : Real := 1.0)
369 return Real
370 is
371 begin
372 if A = 0.0
373 and then Y = 0.0
374 then
375 raise Argument_Error;
377 elsif Y = 0.0 then
378 if A > 0.0 then
379 return 0.0;
380 else -- A < 0.0
381 return Pi;
382 end if;
384 elsif A = 0.0 then
385 if Y > 0.0 then
386 return Half_Pi;
387 else -- Y < 0.0
388 return -Half_Pi;
389 end if;
391 else
392 return Local_Atan (Y, A);
393 end if;
394 end Arctan;
396 -- Arbitrary cycle
398 function Arctan
399 (Y : Real;
400 A : Real := 1.0;
401 Cycle : Real)
402 return Real
403 is
404 begin
405 if Cycle <= 0.0 then
406 raise Argument_Error;
408 elsif A = 0.0
409 and then Y = 0.0
410 then
411 raise Argument_Error;
413 elsif Y = 0.0 then
414 if A > 0.0 then
415 return 0.0;
416 else -- A < 0.0
417 return Cycle / 2.0;
418 end if;
420 elsif A = 0.0 then
421 if Y > 0.0 then
422 return Cycle / 4.0;
423 else -- Y < 0.0
424 return -Cycle / 4.0;
425 end if;
427 else
428 return Local_Atan (Y, A) * Cycle / Two_Pi;
429 end if;
430 end Arctan;
432 -------------
433 -- Arctanh --
434 -------------
436 function Arctanh (A : Real) return Real is
437 begin
438 if abs A = 1.0 then
439 raise Constraint_Error;
441 elsif abs A > 1.0 then
442 raise Argument_Error;
444 elsif abs A < Square_Root_Epsilon then
445 return A;
447 else
448 return 0.5 * Log ((1.0 + A) / (1.0 - A));
449 end if;
450 end Arctanh;
452 ---------
453 -- Cos --
454 ---------
456 -- Natural cycle
458 function Cos (A : Real) return Real is
459 begin
460 if A = 0.0 then
461 return 1.0;
463 elsif abs A < Square_Root_Epsilon then
464 return 1.0;
466 end if;
468 return Real (Aux.Cos (Double (A)));
469 end Cos;
471 -- Arbitrary cycle
473 function Cos (A, Cycle : Real) return Real is
474 T : Real'Base;
476 begin
477 if Cycle <= 0.0 then
478 raise Argument_Error;
480 elsif A = 0.0 then
481 return 1.0;
482 end if;
484 T := Exact_Remainder (abs (A), Cycle) / Cycle;
486 if T = 0.25
487 or else T = 0.75
488 or else T = -0.25
489 or else T = -0.75
490 then
491 return 0.0;
493 elsif T = 0.5 or T = -0.5 then
494 return -1.0;
495 end if;
497 return Real (Aux.Cos (Double (T * Two_Pi)));
498 end Cos;
500 ----------
501 -- Cosh --
502 ----------
504 function Cosh (A : Real) return Real is
505 begin
506 if abs A < Square_Root_Epsilon then
507 return 1.0;
509 elsif abs A > Log_Inverse_Epsilon then
510 return Exp ((abs A) - Log_Two);
511 end if;
513 return Real (Aux.cosh (Double (A)));
515 exception
516 when others =>
517 raise Constraint_Error;
518 end Cosh;
520 ---------
521 -- Cot --
522 ---------
524 -- Natural cycle
526 function Cot (A : Real) return Real is
527 begin
528 if A = 0.0 then
529 raise Constraint_Error;
531 elsif abs A < Square_Root_Epsilon then
532 return 1.0 / A;
533 end if;
535 return Real (1.0 / Real'Base (Aux.tan (Double (A))));
536 end Cot;
538 -- Arbitrary cycle
540 function Cot (A, Cycle : Real) return Real is
541 T : Real'Base;
543 begin
544 if Cycle <= 0.0 then
545 raise Argument_Error;
547 elsif A = 0.0 then
548 raise Constraint_Error;
550 elsif abs A < Square_Root_Epsilon then
551 return 1.0 / A;
552 end if;
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;
558 else
559 return Cos (T * Two_Pi) / Sin (T * Two_Pi);
560 end if;
561 end Cot;
563 ----------
564 -- Coth --
565 ----------
567 function Coth (A : Real) return Real is
568 begin
569 if A = 0.0 then
570 raise Constraint_Error;
572 elsif A < Half_Log_Epsilon then
573 return -1.0;
575 elsif A > -Half_Log_Epsilon then
576 return 1.0;
578 elsif abs A < Square_Root_Epsilon then
579 return 1.0 / A;
580 end if;
582 return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
583 end Coth;
585 ---------------------
586 -- Exact_Remainder --
587 ---------------------
589 function Exact_Remainder
590 (A : Real;
591 Y : Real)
592 return Real
593 is
594 Denominator : Real'Base := abs A;
595 Divisor : Real'Base := abs Y;
596 Reducer : Real'Base;
597 Sign : Real'Base := 1.0;
599 begin
600 if Y = 0.0 then
601 raise Constraint_Error;
603 elsif A = 0.0 then
604 return 0.0;
606 elsif A = Y then
607 return 0.0;
609 elsif Denominator < Divisor then
610 return A;
611 end if;
613 while Denominator >= Divisor loop
615 -- Put divisors mantissa with denominators exponent to make reducer
617 Reducer := Divisor;
619 begin
620 while Reducer * 1_048_576.0 < Denominator loop
621 Reducer := Reducer * 1_048_576.0;
622 end loop;
624 exception
625 when others => null;
626 end;
628 begin
629 while Reducer * 1_024.0 < Denominator loop
630 Reducer := Reducer * 1_024.0;
631 end loop;
633 exception
634 when others => null;
635 end;
637 begin
638 while Reducer * 2.0 < Denominator loop
639 Reducer := Reducer * 2.0;
640 end loop;
642 exception
643 when others => null;
644 end;
646 Denominator := Denominator - Reducer;
647 end loop;
649 if A < 0.0 then
650 return -Denominator;
651 else
652 return Denominator;
653 end if;
654 end Exact_Remainder;
656 ---------
657 -- Exp --
658 ---------
660 function Exp (A : Real) return Real is
661 Result : Real'Base;
663 begin
664 if A = 0.0 then
665 return 1.0;
667 else
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;
674 else
675 return Result;
676 end if;
677 end if;
678 end Exp;
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
688 begin
689 return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
690 end Half_Log_Epsilon;
692 ----------------
693 -- Local_Atan --
694 ----------------
696 function Local_Atan
697 (Y : Real;
698 A : Real := 1.0)
699 return Real
700 is
701 Z : Real'Base;
702 Raw_Atan : Real'Base;
704 begin
705 if abs Y > abs A then
706 Z := abs (A / Y);
707 else
708 Z := abs (Y / A);
709 end if;
711 if Z < Square_Root_Epsilon then
712 Raw_Atan := Z;
714 elsif Z = 1.0 then
715 Raw_Atan := Pi / 4.0;
717 elsif Z < Square_Root_Epsilon then
718 Raw_Atan := Z;
720 else
721 Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
722 end if;
724 if abs Y > abs A then
725 Raw_Atan := Half_Pi - Raw_Atan;
726 end if;
728 if A > 0.0 then
729 if Y > 0.0 then
730 return Raw_Atan;
731 else -- Y < 0.0
732 return -Raw_Atan;
733 end if;
735 else -- A < 0.0
736 if Y > 0.0 then
737 return Pi - Raw_Atan;
738 else -- Y < 0.0
739 return -(Pi - Raw_Atan);
740 end if;
741 end if;
742 end Local_Atan;
744 ---------
745 -- Log --
746 ---------
748 -- Natural base
750 function Log (A : Real) return Real is
751 begin
752 if A < 0.0 then
753 raise Argument_Error;
755 elsif A = 0.0 then
756 raise Constraint_Error;
758 elsif A = 1.0 then
759 return 0.0;
760 end if;
762 return Real (Aux.Log (Double (A)));
763 end Log;
765 -- Arbitrary base
767 function Log (A, Base : Real) return Real is
768 begin
769 if A < 0.0 then
770 raise Argument_Error;
772 elsif Base <= 0.0 or else Base = 1.0 then
773 raise Argument_Error;
775 elsif A = 0.0 then
776 raise Constraint_Error;
778 elsif A = 1.0 then
779 return 0.0;
780 end if;
782 return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
783 end Log;
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
793 begin
794 return Real (Aux.Log (DIEpsilon));
795 end Log_Inverse_Epsilon;
797 ---------
798 -- Sin --
799 ---------
801 -- Natural cycle
803 function Sin (A : Real) return Real is
804 begin
805 if abs A < Square_Root_Epsilon then
806 return A;
807 end if;
809 return Real (Aux.Sin (Double (A)));
810 end Sin;
812 -- Arbitrary cycle
814 function Sin (A, Cycle : Real) return Real is
815 T : Real'Base;
817 begin
818 if Cycle <= 0.0 then
819 raise Argument_Error;
821 elsif A = 0.0 then
822 return A;
823 end if;
825 T := Exact_Remainder (A, Cycle) / Cycle;
827 if T = 0.0 or T = 0.5 or T = -0.5 then
828 return 0.0;
830 elsif T = 0.25 or T = -0.75 then
831 return 1.0;
833 elsif T = -0.25 or T = 0.75 then
834 return -1.0;
836 end if;
838 return Real (Aux.Sin (Double (T * Two_Pi)));
839 end Sin;
841 ----------
842 -- Sinh --
843 ----------
845 function Sinh (A : Real) return Real is
846 begin
847 if abs A < Square_Root_Epsilon then
848 return A;
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);
855 end if;
857 return Real (Aux.Sinh (Double (A)));
859 exception
860 when others =>
861 raise Constraint_Error;
862 end Sinh;
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
872 begin
873 return Real (Aux.Sqrt (DEpsilon));
874 end Square_Root_Epsilon;
876 ----------
877 -- Sqrt --
878 ----------
880 function Sqrt (A : Real) return Real is
881 begin
882 if A < 0.0 then
883 raise Argument_Error;
885 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
887 elsif A = 0.0 then
888 return A;
890 -- Sqrt (1.0) must be exact for good complex accuracy
892 elsif A = 1.0 then
893 return 1.0;
895 end if;
897 return Real (Aux.Sqrt (Double (A)));
898 end Sqrt;
900 ---------
901 -- Tan --
902 ---------
904 -- Natural cycle
906 function Tan (A : Real) return Real is
907 begin
908 if abs A < Square_Root_Epsilon then
909 return A;
911 elsif abs A = Pi / 2.0 then
912 raise Constraint_Error;
913 end if;
915 return Real (Aux.tan (Double (A)));
916 end Tan;
918 -- Arbitrary cycle
920 function Tan (A, Cycle : Real) return Real is
921 T : Real'Base;
923 begin
924 if Cycle <= 0.0 then
925 raise Argument_Error;
927 elsif A = 0.0 then
928 return A;
929 end if;
931 T := Exact_Remainder (A, Cycle) / Cycle;
933 if T = 0.25
934 or else T = 0.75
935 or else T = -0.25
936 or else T = -0.75
937 then
938 raise Constraint_Error;
940 else
941 return Sin (T * Two_Pi) / Cos (T * Two_Pi);
942 end if;
943 end Tan;
945 ----------
946 -- Tanh --
947 ----------
949 function Tanh (A : Real) return Real is
950 begin
951 if A < Half_Log_Epsilon then
952 return -1.0;
954 elsif A > -Half_Log_Epsilon then
955 return 1.0;
957 elsif abs A < Square_Root_Epsilon then
958 return A;
959 end if;
961 return Real (Aux.tanh (Double (A)));
962 end Tanh;
964 ----------------------------
965 -- DEC-Specific functions --
966 ----------------------------
968 function LOG10 (A : REAL) return REAL is
969 begin
970 return Log (A, 10.0);
971 end LOG10;
973 function LOG2 (A : REAL) return REAL is
974 begin
975 return Log (A, 2.0);
976 end LOG2;
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
982 begin
983 return Arctan (A, 1.0);
984 end ATAN;
986 function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
988 function SIND (A : REAL) return REAL is
989 begin
990 return Sin (A, 360.0);
991 end SIND;
993 function COSD (A : REAL) return REAL is
994 begin
995 return Cos (A, 360.0);
996 end COSD;
998 function TAND (A : REAL) return REAL is
999 begin
1000 return Tan (A, 360.0);
1001 end TAND;
1003 function ASIND (A : REAL) return REAL is
1004 begin
1005 return Arcsin (A, 360.0);
1006 end ASIND;
1008 function ACOSD (A : REAL) return REAL is
1009 begin
1010 return Arccos (A, 360.0);
1011 end ACOSD;
1013 function Arctan (A : REAL) return REAL is
1014 begin
1015 return Arctan (A, 1.0, 360.0);
1016 end Arctan;
1018 function ATAND (A : REAL) return REAL is
1019 begin
1020 return Arctan (A, 1.0, 360.0);
1021 end ATAND;
1023 function ATAN2D (A1, A2 : REAL) return REAL is
1024 begin
1025 return Arctan (A1, A2, 360.0);
1026 end ATAN2D;
1028 end Math_Lib;