| blob | 64fed359913ec7a082f90ff53ff3a2e19f18ac38 |
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUNTIME COMPONENTS --
4 -- --
5 -- M A T H _ L I B --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 1992-2000 Free Software Foundation, Inc. --
10 -- --
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, 59 Temple Place - Suite 330, Boston, --
20 -- MA 02111-1307, USA. --
21 -- --
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. --
28 -- --
29 -- GNAT was originally developed by the GNAT team at New York University. --
30 -- Extensive contributions were provided by Ada Core Technologies Inc. --
31 -- --
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.
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 -----------------------
72 function Arctan
73 (Y : Real;
74 A : Real := 1.0)
75 return Real;
77 function Arctan
78 (Y : Real;
79 A : Real := 1.0;
80 Cycle : Real)
81 return Real;
83 function Exact_Remainder
84 (A : Real;
85 Y : Real)
86 return Real;
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)
92 function Local_Atan
93 (Y : Real;
94 A : Real := 1.0)
95 return Real;
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)
104 ----------
105 -- "**" --
106 ----------
108 function "**" (A1, A2 : Real) return Real is
110 begin
111 if A1 = 0.0
112 and then A2 = 0.0
113 then
114 raise Argument_Error;
116 elsif A1 < 0.0 then
117 raise Argument_Error;
119 elsif A2 = 0.0 then
120 return 1.0;
122 elsif A1 = 0.0 then
123 if A2 < 0.0 then
124 raise Constraint_Error;
125 else
126 return 0.0;
127 end if;
129 elsif A1 = 1.0 then
130 return 1.0;
132 elsif A2 = 1.0 then
133 return A1;
135 else
136 begin
137 if A2 = 2.0 then
138 return A1 * A1;
139 else
140 return
141 Real (Aux.pow (Double (A1), Double (A2)));
142 end if;
144 exception
145 when others =>
146 raise Constraint_Error;
147 end;
148 end if;
149 end "**";
151 ------------
152 -- Arccos --
153 ------------
155 -- Natural cycle
157 function Arccos (A : Real) return Real is
158 Temp : Real'Base;
160 begin
161 if abs A > 1.0 then
162 raise Argument_Error;
164 elsif abs A < Square_Root_Epsilon then
165 return Pi / 2.0 - A;
167 elsif A = 1.0 then
168 return 0.0;
170 elsif A = -1.0 then
171 return Pi;
172 end if;
174 Temp := Real (Aux.acos (Double (A)));
176 if Temp < 0.0 then
177 Temp := Pi + Temp;
178 end if;
180 return Temp;
181 end Arccos;
183 -- Arbitrary cycle
185 function Arccos (A, Cycle : Real) return Real is
186 Temp : Real'Base;
188 begin
189 if Cycle <= 0.0 then
190 raise Argument_Error;
192 elsif abs A > 1.0 then
193 raise Argument_Error;
195 elsif abs A < Square_Root_Epsilon then
196 return Cycle / 4.0;
198 elsif A = 1.0 then
199 return 0.0;
201 elsif A = -1.0 then
202 return Cycle / 2.0;
203 end if;
205 Temp := Arctan (Sqrt (1.0 - A * A) / A, 1.0, Cycle);
207 if Temp < 0.0 then
208 Temp := Cycle / 2.0 + Temp;
209 end if;
211 return Temp;
212 end Arccos;
214 -------------
215 -- Arccosh --
216 -------------
218 function Arccosh (A : Real) return Real is
219 begin
220 -- Return Log (A - Sqrt (A * A - 1.0)); double valued,
221 -- only positive value returned
222 -- What is this comment ???
224 if A < 1.0 then
225 raise Argument_Error;
227 elsif A < 1.0 + Square_Root_Epsilon then
228 return A - 1.0;
230 elsif abs A > 1.0 / Square_Root_Epsilon then
231 return Log (A) + Log_Two;
233 else
234 return Log (A + Sqrt (A * A - 1.0));
235 end if;
236 end Arccosh;
238 ------------
239 -- Arccot --
240 ------------
242 -- Natural cycle
244 function Arccot
245 (A : Real;
246 Y : Real := 1.0)
247 return Real
248 is
249 begin
250 -- Just reverse arguments
252 return Arctan (Y, A);
253 end Arccot;
255 -- Arbitrary cycle
257 function Arccot
258 (A : Real;
259 Y : Real := 1.0;
260 Cycle : Real)
261 return Real
262 is
263 begin
264 -- Just reverse arguments
266 return Arctan (Y, A, Cycle);
267 end Arccot;
269 -------------
270 -- Arccoth --
271 -------------
273 function Arccoth (A : Real) return Real is
274 begin
275 if abs A = 1.0 then
276 raise Constraint_Error;
278 elsif abs A < 1.0 then
279 raise Argument_Error;
281 elsif abs A > 1.0 / Epsilon then
282 return 0.0;
284 else
285 return 0.5 * Log ((1.0 + A) / (A - 1.0));
286 end if;
287 end Arccoth;
289 ------------
290 -- Arcsin --
291 ------------
293 -- Natural cycle
295 function Arcsin (A : Real) return Real is
296 begin
297 if abs A > 1.0 then
298 raise Argument_Error;
300 elsif abs A < Square_Root_Epsilon then
301 return A;
303 elsif A = 1.0 then
304 return Pi / 2.0;
306 elsif A = -1.0 then
307 return -Pi / 2.0;
308 end if;
310 return Real (Aux.asin (Double (A)));
311 end Arcsin;
313 -- Arbitrary cycle
315 function Arcsin (A, Cycle : Real) return Real is
316 begin
317 if Cycle <= 0.0 then
318 raise Argument_Error;
320 elsif abs A > 1.0 then
321 raise Argument_Error;
323 elsif A = 0.0 then
324 return A;
326 elsif A = 1.0 then
327 return Cycle / 4.0;
329 elsif A = -1.0 then
330 return -Cycle / 4.0;
331 end if;
333 return Arctan (A / Sqrt (1.0 - A * A), 1.0, Cycle);
334 end Arcsin;
336 -------------
337 -- Arcsinh --
338 -------------
340 function Arcsinh (A : Real) return Real is
341 begin
342 if abs A < Square_Root_Epsilon then
343 return A;
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);
351 elsif A < 0.0 then
352 return -Log (abs A + Sqrt (A * A + 1.0));
354 else
355 return Log (A + Sqrt (A * A + 1.0));
356 end if;
357 end Arcsinh;
359 ------------
360 -- Arctan --
361 ------------
363 -- Natural cycle
365 function Arctan
366 (Y : Real;
367 A : Real := 1.0)
368 return Real
369 is
370 begin
371 if A = 0.0
372 and then Y = 0.0
373 then
374 raise Argument_Error;
376 elsif Y = 0.0 then
377 if A > 0.0 then
378 return 0.0;
379 else -- A < 0.0
380 return Pi;
381 end if;
383 elsif A = 0.0 then
384 if Y > 0.0 then
385 return Half_Pi;
386 else -- Y < 0.0
387 return -Half_Pi;
388 end if;
390 else
391 return Local_Atan (Y, A);
392 end if;
393 end Arctan;
395 -- Arbitrary cycle
397 function Arctan
398 (Y : Real;
399 A : Real := 1.0;
400 Cycle : Real)
401 return Real
402 is
403 begin
404 if Cycle <= 0.0 then
405 raise Argument_Error;
407 elsif A = 0.0
408 and then Y = 0.0
409 then
410 raise Argument_Error;
412 elsif Y = 0.0 then
413 if A > 0.0 then
414 return 0.0;
415 else -- A < 0.0
416 return Cycle / 2.0;
417 end if;
419 elsif A = 0.0 then
420 if Y > 0.0 then
421 return Cycle / 4.0;
422 else -- Y < 0.0
423 return -Cycle / 4.0;
424 end if;
426 else
427 return Local_Atan (Y, A) * Cycle / Two_Pi;
428 end if;
429 end Arctan;
431 -------------
432 -- Arctanh --
433 -------------
435 function Arctanh (A : Real) return Real is
436 begin
437 if abs A = 1.0 then
438 raise Constraint_Error;
440 elsif abs A > 1.0 then
441 raise Argument_Error;
443 elsif abs A < Square_Root_Epsilon then
444 return A;
446 else
447 return 0.5 * Log ((1.0 + A) / (1.0 - A));
448 end if;
449 end Arctanh;
451 ---------
452 -- Cos --
453 ---------
455 -- Natural cycle
457 function Cos (A : Real) return Real is
458 begin
459 if A = 0.0 then
460 return 1.0;
462 elsif abs A < Square_Root_Epsilon then
463 return 1.0;
465 end if;
467 return Real (Aux.Cos (Double (A)));
468 end Cos;
470 -- Arbitrary cycle
472 function Cos (A, Cycle : Real) return Real is
473 T : Real'Base;
475 begin
476 if Cycle <= 0.0 then
477 raise Argument_Error;
479 elsif A = 0.0 then
480 return 1.0;
481 end if;
483 T := Exact_Remainder (abs (A), Cycle) / Cycle;
485 if T = 0.25
486 or else T = 0.75
487 or else T = -0.25
488 or else T = -0.75
489 then
490 return 0.0;
492 elsif T = 0.5 or T = -0.5 then
493 return -1.0;
494 end if;
496 return Real (Aux.Cos (Double (T * Two_Pi)));
497 end Cos;
499 ----------
500 -- Cosh --
501 ----------
503 function Cosh (A : Real) return Real is
504 begin
505 if abs A < Square_Root_Epsilon then
506 return 1.0;
508 elsif abs A > Log_Inverse_Epsilon then
509 return Exp ((abs A) - Log_Two);
510 end if;
512 return Real (Aux.cosh (Double (A)));
514 exception
515 when others =>
516 raise Constraint_Error;
517 end Cosh;
519 ---------
520 -- Cot --
521 ---------
523 -- Natural cycle
525 function Cot (A : Real) return Real is
526 begin
527 if A = 0.0 then
528 raise Constraint_Error;
530 elsif abs A < Square_Root_Epsilon then
531 return 1.0 / A;
532 end if;
534 return Real (1.0 / Real'Base (Aux.tan (Double (A))));
535 end Cot;
537 -- Arbitrary cycle
539 function Cot (A, Cycle : Real) return Real is
540 T : Real'Base;
542 begin
543 if Cycle <= 0.0 then
544 raise Argument_Error;
546 elsif A = 0.0 then
547 raise Constraint_Error;
549 elsif abs A < Square_Root_Epsilon then
550 return 1.0 / A;
551 end if;
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;
557 else
558 return Cos (T * Two_Pi) / Sin (T * Two_Pi);
559 end if;
560 end Cot;
562 ----------
563 -- Coth --
564 ----------
566 function Coth (A : Real) return Real is
567 begin
568 if A = 0.0 then
569 raise Constraint_Error;
571 elsif A < Half_Log_Epsilon then
572 return -1.0;
574 elsif A > -Half_Log_Epsilon then
575 return 1.0;
577 elsif abs A < Square_Root_Epsilon then
578 return 1.0 / A;
579 end if;
581 return Real (1.0 / Real'Base (Aux.tanh (Double (A))));
582 end Coth;
584 ---------------------
585 -- Exact_Remainder --
586 ---------------------
588 function Exact_Remainder
589 (A : Real;
590 Y : Real)
591 return Real
592 is
593 Denominator : Real'Base := abs A;
594 Divisor : Real'Base := abs Y;
595 Reducer : Real'Base;
596 Sign : Real'Base := 1.0;
598 begin
599 if Y = 0.0 then
600 raise Constraint_Error;
602 elsif A = 0.0 then
603 return 0.0;
605 elsif A = Y then
606 return 0.0;
608 elsif Denominator < Divisor then
609 return A;
610 end if;
612 while Denominator >= Divisor loop
614 -- Put divisors mantissa with denominators exponent to make reducer
616 Reducer := Divisor;
618 begin
619 while Reducer * 1_048_576.0 < Denominator loop
620 Reducer := Reducer * 1_048_576.0;
621 end loop;
623 exception
624 when others => null;
625 end;
627 begin
628 while Reducer * 1_024.0 < Denominator loop
629 Reducer := Reducer * 1_024.0;
630 end loop;
632 exception
633 when others => null;
634 end;
636 begin
637 while Reducer * 2.0 < Denominator loop
638 Reducer := Reducer * 2.0;
639 end loop;
641 exception
642 when others => null;
643 end;
645 Denominator := Denominator - Reducer;
646 end loop;
648 if A < 0.0 then
649 return -Denominator;
650 else
651 return Denominator;
652 end if;
653 end Exact_Remainder;
655 ---------
656 -- Exp --
657 ---------
659 function Exp (A : Real) return Real is
660 Result : Real'Base;
662 begin
663 if A = 0.0 then
664 return 1.0;
666 else
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;
673 else
674 return Result;
675 end if;
676 end if;
677 end Exp;
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
687 begin
688 return Real (0.5 * Real'Base (Aux.Log (DEpsilon)));
689 end Half_Log_Epsilon;
691 ----------------
692 -- Local_Atan --
693 ----------------
695 function Local_Atan
696 (Y : Real;
697 A : Real := 1.0)
698 return Real
699 is
700 Z : Real'Base;
701 Raw_Atan : Real'Base;
703 begin
704 if abs Y > abs A then
705 Z := abs (A / Y);
706 else
707 Z := abs (Y / A);
708 end if;
710 if Z < Square_Root_Epsilon then
711 Raw_Atan := Z;
713 elsif Z = 1.0 then
714 Raw_Atan := Pi / 4.0;
716 elsif Z < Square_Root_Epsilon then
717 Raw_Atan := Z;
719 else
720 Raw_Atan := Real'Base (Aux.Atan (Double (Z)));
721 end if;
723 if abs Y > abs A then
724 Raw_Atan := Half_Pi - Raw_Atan;
725 end if;
727 if A > 0.0 then
728 if Y > 0.0 then
729 return Raw_Atan;
730 else -- Y < 0.0
731 return -Raw_Atan;
732 end if;
734 else -- A < 0.0
735 if Y > 0.0 then
736 return Pi - Raw_Atan;
737 else -- Y < 0.0
738 return -(Pi - Raw_Atan);
739 end if;
740 end if;
741 end Local_Atan;
743 ---------
744 -- Log --
745 ---------
747 -- Natural base
749 function Log (A : Real) return Real is
750 begin
751 if A < 0.0 then
752 raise Argument_Error;
754 elsif A = 0.0 then
755 raise Constraint_Error;
757 elsif A = 1.0 then
758 return 0.0;
759 end if;
761 return Real (Aux.Log (Double (A)));
762 end Log;
764 -- Arbitrary base
766 function Log (A, Base : Real) return Real is
767 begin
768 if A < 0.0 then
769 raise Argument_Error;
771 elsif Base <= 0.0 or else Base = 1.0 then
772 raise Argument_Error;
774 elsif A = 0.0 then
775 raise Constraint_Error;
777 elsif A = 1.0 then
778 return 0.0;
779 end if;
781 return Real (Aux.Log (Double (A)) / Aux.Log (Double (Base)));
782 end Log;
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
792 begin
793 return Real (Aux.Log (DIEpsilon));
794 end Log_Inverse_Epsilon;
796 ---------
797 -- Sin --
798 ---------
800 -- Natural cycle
802 function Sin (A : Real) return Real is
803 begin
804 if abs A < Square_Root_Epsilon then
805 return A;
806 end if;
808 return Real (Aux.Sin (Double (A)));
809 end Sin;
811 -- Arbitrary cycle
813 function Sin (A, Cycle : Real) return Real is
814 T : Real'Base;
816 begin
817 if Cycle <= 0.0 then
818 raise Argument_Error;
820 elsif A = 0.0 then
821 return A;
822 end if;
824 T := Exact_Remainder (A, Cycle) / Cycle;
826 if T = 0.0 or T = 0.5 or T = -0.5 then
827 return 0.0;
829 elsif T = 0.25 or T = -0.75 then
830 return 1.0;
832 elsif T = -0.25 or T = 0.75 then
833 return -1.0;
835 end if;
837 return Real (Aux.Sin (Double (T * Two_Pi)));
838 end Sin;
840 ----------
841 -- Sinh --
842 ----------
844 function Sinh (A : Real) return Real is
845 begin
846 if abs A < Square_Root_Epsilon then
847 return A;
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);
854 end if;
856 return Real (Aux.Sinh (Double (A)));
858 exception
859 when others =>
860 raise Constraint_Error;
861 end Sinh;
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
871 begin
872 return Real (Aux.Sqrt (DEpsilon));
873 end Square_Root_Epsilon;
875 ----------
876 -- Sqrt --
877 ----------
879 function Sqrt (A : Real) return Real is
880 begin
881 if A < 0.0 then
882 raise Argument_Error;
884 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
886 elsif A = 0.0 then
887 return A;
889 -- Sqrt (1.0) must be exact for good complex accuracy
891 elsif A = 1.0 then
892 return 1.0;
894 end if;
896 return Real (Aux.Sqrt (Double (A)));
897 end Sqrt;
899 ---------
900 -- Tan --
901 ---------
903 -- Natural cycle
905 function Tan (A : Real) return Real is
906 begin
907 if abs A < Square_Root_Epsilon then
908 return A;
910 elsif abs A = Pi / 2.0 then
911 raise Constraint_Error;
912 end if;
914 return Real (Aux.tan (Double (A)));
915 end Tan;
917 -- Arbitrary cycle
919 function Tan (A, Cycle : Real) return Real is
920 T : Real'Base;
922 begin
923 if Cycle <= 0.0 then
924 raise Argument_Error;
926 elsif A = 0.0 then
927 return A;
928 end if;
930 T := Exact_Remainder (A, Cycle) / Cycle;
932 if T = 0.25
933 or else T = 0.75
934 or else T = -0.25
935 or else T = -0.75
936 then
937 raise Constraint_Error;
939 else
940 return Sin (T * Two_Pi) / Cos (T * Two_Pi);
941 end if;
942 end Tan;
944 ----------
945 -- Tanh --
946 ----------
948 function Tanh (A : Real) return Real is
949 begin
950 if A < Half_Log_Epsilon then
951 return -1.0;
953 elsif A > -Half_Log_Epsilon then
954 return 1.0;
956 elsif abs A < Square_Root_Epsilon then
957 return A;
958 end if;
960 return Real (Aux.tanh (Double (A)));
961 end Tanh;
963 ----------------------------
964 -- DEC-Specific functions --
965 ----------------------------
967 function LOG10 (A : REAL) return REAL is
968 begin
969 return Log (A, 10.0);
970 end LOG10;
972 function LOG2 (A : REAL) return REAL is
973 begin
974 return Log (A, 2.0);
975 end LOG2;
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
981 begin
982 return Arctan (A, 1.0);
983 end ATAN;
985 function ATAN2 (A1, A2 : REAL) return REAL renames Arctan;
987 function SIND (A : REAL) return REAL is
988 begin
989 return Sin (A, 360.0);
990 end SIND;
992 function COSD (A : REAL) return REAL is
993 begin
994 return Cos (A, 360.0);
995 end COSD;
997 function TAND (A : REAL) return REAL is
998 begin
999 return Tan (A, 360.0);
1000 end TAND;
1002 function ASIND (A : REAL) return REAL is
1003 begin
1004 return Arcsin (A, 360.0);
1005 end ASIND;
1007 function ACOSD (A : REAL) return REAL is
1008 begin
1009 return Arccos (A, 360.0);
1010 end ACOSD;
1012 function Arctan (A : REAL) return REAL is
1013 begin
1014 return Arctan (A, 1.0, 360.0);
1015 end Arctan;
1017 function ATAND (A : REAL) return REAL is
1018 begin
1019 return Arctan (A, 1.0, 360.0);
1020 end ATAND;
1022 function ATAN2D (A1, A2 : REAL) return REAL is
1023 begin
1024 return Arctan (A1, A2, 360.0);
1025 end ATAN2D;
1027 end Math_Lib;