| blob | af978419ffa36b46e9101fc2bfa42513a4d56d4e |
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUNTIME COMPONENTS --
4 -- --
5 -- M A T H _ L I B --
6 -- --
7 -- B o d y --
8 -- --
9 -- $Revision: 1.1.16.1 $
10 -- --
11 -- Copyright (C) 1992-2000 Free Software Foundation, Inc. --
12 -- --
13 -- GNAT is free software; you can redistribute it and/or modify it under --
14 -- terms of the GNU General Public License as published by the Free Soft- --
15 -- ware Foundation; either version 2, or (at your option) any later ver- --
16 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
17 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
18 -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License --
19 -- for more details. You should have received a copy of the GNU General --
20 -- Public License distributed with GNAT; see file COPYING. If not, write --
21 -- to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, --
22 -- MA 02111-1307, USA. --
23 -- --
24 -- As a special exception, if other files instantiate generics from this --
25 -- unit, or you link this unit with other files to produce an executable, --
26 -- this unit does not by itself cause the resulting executable to be --
27 -- covered by the GNU General Public License. This exception does not --
28 -- however invalidate any other reasons why the executable file might be --
29 -- covered by the GNU Public License. --
30 -- --
31 -- GNAT was originally developed by the GNAT team at New York University. --
32 -- Extensive contributions were provided by Ada Core Technologies Inc. --
33 -- --
34 ------------------------------------------------------------------------------
36 -- This body is specifically for using an Ada interface to C math.h to get
37 -- the computation engine. Many special cases are handled locally to avoid
38 -- unnecessary calls. This is not a "strict" implementation, but takes full
39 -- advantage of the C functions, e.g. in providing interface to hardware
40 -- provided versions of the elementary functions.
42 -- A known weakness is that on the x86, all computation is done in Double,
43 -- which means that a lot of accuracy is lost for the Long_Long_Float case.
45 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
46 -- sinh, cosh, tanh from C library via math.h
48 -- This is an adaptation of Ada.Numerics.Generic_Elementary_Functions that
49 -- provides a compatible body for the DEC Math_Lib package.
63 IEpsilon : constant Real'Base := 1.0 / Epsilon;
65 subtype Double is Aux.Double;
67 DEpsilon : constant Double := Double (Epsilon);
68 DIEpsilon : constant Double := Double (IEpsilon);
70 -----------------------
71 -- Local Subprograms --
72 -----------------------
74 function Arctan
75 (Y : Real;
76 A : Real := 1.0)
77 return Real;
79 function Arctan
80 (Y : Real;
81 A : Real := 1.0;
82 Cycle : Real)
83 return Real;
85 function Exact_Remainder
86 (A : Real;
87 Y : Real)
88 return Real;
89 -- Computes exact remainder of A divided by Y
91 function Half_Log_Epsilon return Real;
92 -- Function to provide constant: 0.5 * Log (Epsilon)
94 function Local_Atan
95 (Y : Real;
96 A : Real := 1.0)
97 return Real;
98 -- Common code for arc tangent after cyele reduction
100 function Log_Inverse_Epsilon return Real;
101 -- Function to provide constant: Log (1.0 / Epsilon)
103 function Square_Root_Epsilon return Real;
104 -- Function to provide constant: Sqrt (Epsilon)
106 ----------
107 -- "**" --
108 ----------
110 function "**" (A1, A2 : Real) return Real is
112 begin
113 if A1 = 0.0
114 and then A2 = 0.0
115 then
116 raise Argument_Error;
118 elsif A1 < 0.0 then
119 raise Argument_Error;
121 elsif A2 = 0.0 then
122 return 1.0;
124 elsif A1 = 0.0 then
125 if A2 < 0.0 then
126 raise Constraint_Error;
127 else
128 return 0.0;
129 end if;
131 elsif A1 = 1.0 then
132 return 1.0;
134 elsif A2 = 1.0 then
135 return A1;
137 else
138 begin
139 if A2 = 2.0 then
140 return A1 * A1;
141 else
142 return
143 Real (Aux.pow (Double (A1), Double (A2)));
144 end if;
146 exception
147 when others =>
148 raise Constraint_Error;
149 end;
150 end if;
151 end "**";
153 ------------
154 -- Arccos --
155 ------------
157 -- Natural cycle
159 function Arccos (A : Real) return Real is
160 Temp : Real'Base;
162 begin
163 if abs A > 1.0 then
164 raise Argument_Error;
166 elsif abs A < Square_Root_Epsilon then
167 return Pi / 2.0 - A;
169 elsif A = 1.0 then
170 return 0.0;
172 elsif A = -1.0 then
173 return Pi;
174 end if;
176 Temp := Real (Aux.acos (Double (A)));
178 if Temp < 0.0 then
179 Temp := Pi + Temp;
180 end if;
182 return Temp;
183 end Arccos;
185 -- Arbitrary cycle
187 function Arccos (A, Cycle : Real) return Real is
188 Temp : Real'Base;
190 begin
191 if Cycle <= 0.0 then
192 raise Argument_Error;
194 elsif abs A > 1.0 then
195 raise Argument_Error;
197 elsif abs A < Square_Root_Epsilon then
198 return Cycle / 4.0;
200 elsif A = 1.0 then
201 return 0.0;
203 elsif A = -1.0 then
204 return Cycle / 2.0;
205 end if;
207 Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
209 if Temp < 0.0 then
210 Temp := Cycle / 2.0 + Temp;
211 end if;
213 return Temp;
214 end Arccos;
216 -------------
217 -- Arccosh --
218 -------------
220 function Arccosh (A : Real) return Real is
221 begin
222 -- Return Log (A - Sqrt (A * A - 1.0)); double valued,
223 -- only positive value returned
224 -- What is this comment ???
226 if A < 1.0 then
227 raise Argument_Error;
229 elsif A < 1.0 + Square_Root_Epsilon then
230 return A - 1.0;
232 elsif abs A > 1.0 / Square_Root_Epsilon then
233 return Log (A) + Log_Two;
235 else
236 return Log (A + Sqrt (A * A - 1.0));
237 end if;
238 end Arccosh;
240 ------------
241 -- Arccot --
242 ------------
244 -- Natural cycle
246 function Arccot
247 (A : Real;
248 Y : Real := 1.0)
249 return Real
250 is
251 begin
252 -- Just reverse arguments
254 return Arctan (Y, A);
255 end Arccot;
257 -- Arbitrary cycle
259 function Arccot
260 (A : Real;
261 Y : Real := 1.0;
262 Cycle : Real)
263 return Real
264 is
265 begin
266 -- Just reverse arguments
268 return Arctan (Y, A, Cycle);
269 end Arccot;
271 -------------
272 -- Arccoth --
273 -------------
275 function Arccoth (A : Real) return Real is
276 begin
277 if abs A = 1.0 then
278 raise Constraint_Error;
280 elsif abs A < 1.0 then
281 raise Argument_Error;
283 elsif abs A > 1.0 / Epsilon then
284 return 0.0;
286 else
287 return 0.5 * Log ((1.0 + A) / (A - 1.0));
288 end if;
289 end Arccoth;
291 ------------
292 -- Arcsin --
293 ------------
295 -- Natural cycle
297 function Arcsin (A : Real) return Real is
298 begin
299 if abs A > 1.0 then
300 raise Argument_Error;
302 elsif abs A < Square_Root_Epsilon then
303 return A;
305 elsif A = 1.0 then
306 return Pi / 2.0;
308 elsif A = -1.0 then
309 return -Pi / 2.0;
310 end if;
312 return Real (Aux.asin (Double (A)));
313 end Arcsin;
315 -- Arbitrary cycle
317 function Arcsin (A, Cycle : Real) return Real is
318 begin
319 if Cycle <= 0.0 then
320 raise Argument_Error;
322 elsif abs A > 1.0 then
323 raise Argument_Error;
325 elsif A = 0.0 then
326 return A;
328 elsif A = 1.0 then
329 return Cycle / 4.0;
331 elsif A = -1.0 then
332 return -Cycle / 4.0;
333 end if;
335 return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
336 end Arcsin;
338 -------------
339 -- Arcsinh --
340 -------------
342 function Arcsinh (A : Real) return Real is
343 begin
344 if abs A < Square_Root_Epsilon then
345 return A;
347 elsif A > 1.0 / Square_Root_Epsilon then
348 return Log (A) + Log_Two;
350 elsif A < -1.0 / Square_Root_Epsilon then
351 return -(Log (-A) + Log_Two);
353 elsif A < 0.0 then
354 return -Log (abs A + Sqrt (A * A + 1.0));
356 else
357 return Log (A + Sqrt (A * A + 1.0));
358 end if;
359 end Arcsinh;
361 ------------
362 -- Arctan --
363 ------------
365 -- Natural cycle
367 function Arctan
368 (Y : Real;
369 A : Real := 1.0)
370 return Real
371 is
372 begin
373 if A = 0.0
374 and then Y = 0.0
375 then
376 raise Argument_Error;
378 elsif Y = 0.0 then
379 if A > 0.0 then
380 return 0.0;
381 else -- A < 0.0
382 return Pi;
383 end if;
385 elsif A = 0.0 then
386 if Y > 0.0 then
387 return Half_Pi;
388 else -- Y < 0.0
389 return -Half_Pi;
390 end if;
392 else
393 return Local_Atan (Y, A);
394 end if;
395 end Arctan;
397 -- Arbitrary cycle
399 function Arctan
400 (Y : Real;
401 A : Real := 1.0;
402 Cycle : Real)
403 return Real
404 is
405 begin
406 if Cycle <= 0.0 then
407 raise Argument_Error;
409 elsif A = 0.0
410 and then Y = 0.0
411 then
412 raise Argument_Error;
414 elsif Y = 0.0 then
415 if A > 0.0 then
416 return 0.0;
417 else -- A < 0.0
418 return Cycle / 2.0;
419 end if;
421 elsif A = 0.0 then
422 if Y > 0.0 then
423 return Cycle / 4.0;
424 else -- Y < 0.0
425 return -Cycle / 4.0;
426 end if;
428 else
429 return Local_Atan (Y, A) * Cycle / Two_Pi;
430 end if;
431 end Arctan;
433 -------------
434 -- Arctanh --
435 -------------
437 function Arctanh (A : Real) return Real is
438 begin
439 if abs A = 1.0 then
440 raise Constraint_Error;
442 elsif abs A > 1.0 then
443 raise Argument_Error;
445 elsif abs A < Square_Root_Epsilon then
446 return A;
448 else
449 return 0.5 * Log ((1.0 + A) / (1.0 - A));
450 end if;
451 end Arctanh;
453 ---------
454 -- Cos --
455 ---------
457 -- Natural cycle
459 function Cos (A : Real) return Real is
460 begin
461 if A = 0.0 then
462 return 1.0;
464 elsif abs A < Square_Root_Epsilon then
465 return 1.0;
467 end if;
469 return Real (Aux.Cos (Double (A)));
470 end Cos;
472 -- Arbitrary cycle
474 function Cos (A, Cycle : Real) return Real is
475 T : Real'Base;
477 begin
478 if Cycle <= 0.0 then
479 raise Argument_Error;
481 elsif A = 0.0 then
482 return 1.0;
483 end if;
485 T := Exact_Remainder (abs (A), Cycle) / Cycle;
487 if T = 0.25
488 or else T = 0.75
489 or else T = -0.25
490 or else T = -0.75
491 then
492 return 0.0;
494 elsif T = 0.5 or T = -0.5 then
495 return -1.0;
496 end if;
498 return Real (Aux.Cos (Double (T * Two_Pi)));
499 end Cos;
501 ----------
502 -- Cosh --
503 ----------
505 function Cosh (A : Real) return Real is
506 begin
507 if abs A < Square_Root_Epsilon then
508 return 1.0;
510 elsif abs A > Log_Inverse_Epsilon then
511 return Exp ((abs A) - Log_Two);
512 end if;
514 return Real (Aux.cosh (Double (A)));
516 exception
517 when others =>
518 raise Constraint_Error;
519 end Cosh;
521 ---------
522 -- Cot --
523 ---------
525 -- Natural cycle
527 function Cot (A : Real) return Real is
528 begin
529 if A = 0.0 then
530 raise Constraint_Error;
532 elsif abs A < Square_Root_Epsilon then
533 return 1.0 / A;
534 end if;
536 return Real (1.0 / Real'Base (Aux.tan (Double (A))));
537 end Cot;
539 -- Arbitrary cycle
541 function Cot (A, Cycle : Real) return Real is
542 T : Real'Base;
544 begin
545 if Cycle <= 0.0 then
546 raise Argument_Error;
548 elsif A = 0.0 then
549 raise Constraint_Error;
551 elsif abs A < Square_Root_Epsilon then
552 return 1.0 / A;
553 end if;
555 T := Exact_Remainder (A, Cycle) / Cycle;
557 if T = 0.0 or T = 0.5 or T = -0.5 then
558 raise Constraint_Error;
559 else
560 return Cos (T * Two_Pi) / Sin (T * Two_Pi);
561 end if;
562 end Cot;
564 ----------
565 -- Coth --
566 ----------
568 function Coth (A : Real) return Real is
569 begin
570 if A = 0.0 then
571 raise Constraint_Error;
573 elsif A < Half_Log_Epsilon then
574 return -1.0;
576 elsif A > -Half_Log_Epsilon then
577 return 1.0;
579 elsif abs A < Square_Root_Epsilon then
580 return 1.0 / A;
581 end if;
583 return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
584 end Coth;
586 ---------------------
587 -- Exact_Remainder --
588 ---------------------
590 function Exact_Remainder
591 (A : Real;
592 Y : Real)
593 return Real
594 is
595 Denominator : Real'Base := abs A;
596 Divisor : Real'Base := abs Y;
597 Reducer : Real'Base;
598 Sign : Real'Base := 1.0;
600 begin
601 if Y = 0.0 then
602 raise Constraint_Error;
604 elsif A = 0.0 then
605 return 0.0;
607 elsif A = Y then
608 return 0.0;
610 elsif Denominator < Divisor then
611 return A;
612 end if;
614 while Denominator >= Divisor loop
616 -- Put divisors mantissa with denominators exponent to make reducer
618 Reducer := Divisor;
620 begin
621 while Reducer * 1_048_576.0 < Denominator loop
622 Reducer := Reducer * 1_048_576.0;
623 end loop;
625 exception
626 when others => null;
627 end;
629 begin
630 while Reducer * 1_024.0 < Denominator loop
631 Reducer := Reducer * 1_024.0;
632 end loop;
634 exception
635 when others => null;
636 end;
638 begin
639 while Reducer * 2.0 < Denominator loop
640 Reducer := Reducer * 2.0;
641 end loop;
643 exception
644 when others => null;
645 end;
647 Denominator := Denominator - Reducer;
648 end loop;
650 if A < 0.0 then
651 return -Denominator;
652 else
653 return Denominator;
654 end if;
655 end Exact_Remainder;
657 ---------
658 -- Exp --
659 ---------
661 function Exp (A : Real) return Real is
662 Result : Real'Base;
664 begin
665 if A = 0.0 then
666 return 1.0;
668 else
669 Result := Real (Aux.Exp (Double (A)));
671 -- The check here catches the case of Exp returning IEEE infinity
673 if Result > Real'Last then
674 raise Constraint_Error;
675 else
676 return Result;
677 end if;
678 end if;
679 end Exp;
681 ----------------------
682 -- Half_Log_Epsilon --
683 ----------------------
685 -- Cannot precompute this constant, because this is required to be a
686 -- pure package, which allows no state. A pity, but no way around it!
688 function Half_Log_Epsilon return Real is
689 begin
690 return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
691 end Half_Log_Epsilon;
693 ----------------
694 -- Local_Atan --
695 ----------------
697 function Local_Atan
698 (Y : Real;
699 A : Real := 1.0)
700 return Real
701 is
702 Z : Real'Base;
703 Raw_Atan : Real'Base;
705 begin
706 if abs Y > abs A then
707 Z := abs (A / Y);
708 else
709 Z := abs (Y / A);
710 end if;
712 if Z < Square_Root_Epsilon then
713 Raw_Atan := Z;
715 elsif Z = 1.0 then
716 Raw_Atan := Pi / 4.0;
718 elsif Z < Square_Root_Epsilon then
719 Raw_Atan := Z;
721 else
722 Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
723 end if;
725 if abs Y > abs A then
726 Raw_Atan := Half_Pi - Raw_Atan;
727 end if;
729 if A > 0.0 then
730 if Y > 0.0 then
731 return Raw_Atan;
732 else -- Y < 0.0
733 return -Raw_Atan;
734 end if;
736 else -- A < 0.0
737 if Y > 0.0 then
738 return Pi - Raw_Atan;
739 else -- Y < 0.0
740 return -(Pi - Raw_Atan);
741 end if;
742 end if;
743 end Local_Atan;
745 ---------
746 -- Log --
747 ---------
749 -- Natural base
751 function Log (A : Real) return Real is
752 begin
753 if A < 0.0 then
754 raise Argument_Error;
756 elsif A = 0.0 then
757 raise Constraint_Error;
759 elsif A = 1.0 then
760 return 0.0;
761 end if;
763 return Real (Aux.Log (Double (A)));
764 end Log;
766 -- Arbitrary base
768 function Log (A, Base : Real) return Real is
769 begin
770 if A < 0.0 then
771 raise Argument_Error;
773 elsif Base <= 0.0 or else Base = 1.0 then
774 raise Argument_Error;
776 elsif A = 0.0 then
777 raise Constraint_Error;
779 elsif A = 1.0 then
780 return 0.0;
781 end if;
783 return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
784 end Log;
786 -------------------------
787 -- Log_Inverse_Epsilon --
788 -------------------------
790 -- Cannot precompute this constant, because this is required to be a
791 -- pure package, which allows no state. A pity, but no way around it!
793 function Log_Inverse_Epsilon return Real is
794 begin
795 return Real (Aux.Log (DIEpsilon));
796 end Log_Inverse_Epsilon;
798 ---------
799 -- Sin --
800 ---------
802 -- Natural cycle
804 function Sin (A : Real) return Real is
805 begin
806 if abs A < Square_Root_Epsilon then
807 return A;
808 end if;
810 return Real (Aux.Sin (Double (A)));
811 end Sin;
813 -- Arbitrary cycle
815 function Sin (A, Cycle : Real) return Real is
816 T : Real'Base;
818 begin
819 if Cycle <= 0.0 then
820 raise Argument_Error;
822 elsif A = 0.0 then
823 return A;
824 end if;
826 T := Exact_Remainder (A, Cycle) / Cycle;
828 if T = 0.0 or T = 0.5 or T = -0.5 then
829 return 0.0;
831 elsif T = 0.25 or T = -0.75 then
832 return 1.0;
834 elsif T = -0.25 or T = 0.75 then
835 return -1.0;
837 end if;
839 return Real (Aux.Sin (Double (T * Two_Pi)));
840 end Sin;
842 ----------
843 -- Sinh --
844 ----------
846 function Sinh (A : Real) return Real is
847 begin
848 if abs A < Square_Root_Epsilon then
849 return A;
851 elsif A > Log_Inverse_Epsilon then
852 return Exp (A - Log_Two);
854 elsif A < -Log_Inverse_Epsilon then
855 return -Exp ((-A) - Log_Two);
856 end if;
858 return Real (Aux.Sinh (Double (A)));
860 exception
861 when others =>
862 raise Constraint_Error;
863 end Sinh;
865 -------------------------
866 -- Square_Root_Epsilon --
867 -------------------------
869 -- Cannot precompute this constant, because this is required to be a
870 -- pure package, which allows no state. A pity, but no way around it!
872 function Square_Root_Epsilon return Real is
873 begin
874 return Real (Aux.Sqrt (DEpsilon));
875 end Square_Root_Epsilon;
877 ----------
878 -- Sqrt --
879 ----------
881 function Sqrt (A : Real) return Real is
882 begin
883 if A < 0.0 then
884 raise Argument_Error;
886 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
888 elsif A = 0.0 then
889 return A;
891 -- Sqrt (1.0) must be exact for good complex accuracy
893 elsif A = 1.0 then
894 return 1.0;
896 end if;
898 return Real (Aux.Sqrt (Double (A)));
899 end Sqrt;
901 ---------
902 -- Tan --
903 ---------
905 -- Natural cycle
907 function Tan (A : Real) return Real is
908 begin
909 if abs A < Square_Root_Epsilon then
910 return A;
912 elsif abs A = Pi / 2.0 then
913 raise Constraint_Error;
914 end if;
916 return Real (Aux.tan (Double (A)));
917 end Tan;
919 -- Arbitrary cycle
921 function Tan (A, Cycle : Real) return Real is
922 T : Real'Base;
924 begin
925 if Cycle <= 0.0 then
926 raise Argument_Error;
928 elsif A = 0.0 then
929 return A;
930 end if;
932 T := Exact_Remainder (A, Cycle) / Cycle;
934 if T = 0.25
935 or else T = 0.75
936 or else T = -0.25
937 or else T = -0.75
938 then
939 raise Constraint_Error;
941 else
942 return Sin (T * Two_Pi) / Cos (T * Two_Pi);
943 end if;
944 end Tan;
946 ----------
947 -- Tanh --
948 ----------
950 function Tanh (A : Real) return Real is
951 begin
952 if A < Half_Log_Epsilon then
953 return -1.0;
955 elsif A > -Half_Log_Epsilon then
956 return 1.0;
958 elsif abs A < Square_Root_Epsilon then
959 return A;
960 end if;
962 return Real (Aux.tanh (Double (A)));
963 end Tanh;
965 ----------------------------
966 -- DEC-Specific functions --
967 ----------------------------
969 function LOG10 (A : REAL) return REAL is
970 begin
971 return Log (A, 10.0);
972 end LOG10;
974 function LOG2 (A : REAL) return REAL is
975 begin
976 return Log (A, 2.0);
977 end LOG2;
979 function ASIN (A : REAL) return REAL renames Arcsin;
980 function ACOS (A : REAL) return REAL renames Arccos;
982 function ATAN (A : REAL) return REAL is
983 begin
984 return Arctan (A, 1.0);
985 end ATAN;
987 function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
989 function SIND (A : REAL) return REAL is
990 begin
991 return Sin (A, 360.0);
992 end SIND;
994 function COSD (A : REAL) return REAL is
995 begin
996 return Cos (A, 360.0);
997 end COSD;
999 function TAND (A : REAL) return REAL is
1000 begin
1001 return Tan (A, 360.0);
1002 end TAND;
1004 function ASIND (A : REAL) return REAL is
1005 begin
1006 return Arcsin (A, 360.0);
1007 end ASIND;
1009 function ACOSD (A : REAL) return REAL is
1010 begin
1011 return Arccos (A, 360.0);
1012 end ACOSD;
1014 function Arctan (A : REAL) return REAL is
1015 begin
1016 return Arctan (A, 1.0, 360.0);
1017 end Arctan;
1019 function ATAND (A : REAL) return REAL is
1020 begin
1021 return Arctan (A, 1.0, 360.0);
1022 end ATAND;
1024 function ATAN2D (A1, A2 : REAL) return REAL is
1025 begin
1026 return Arctan (A1, A2, 360.0);
1027 end ATAN2D;
1029 end Math_Lib;