Dead
[official-gcc.git] / gomp-20050608-branch / gcc / ada / a-ngelfu.adb
blob97d170dd474762a4c6d581ab8ecb8326fc36091c
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
5 -- ADA.NUMERICS.GENERIC_ELEMENTARY_FUNCTIONS --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 1992-2005, 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, 51 Franklin Street, Fifth Floor, --
20 -- Boston, MA 02110-1301, 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 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan,
41 -- sinh, cosh, tanh from C library via math.h
43 with Ada.Numerics.Aux;
45 package body Ada.Numerics.Generic_Elementary_Functions is
47 use type Ada.Numerics.Aux.Double;
49 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
50 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
51 Half_Log_Two : constant := Log_Two / 2;
53 subtype T is Float_Type'Base;
54 subtype Double is Aux.Double;
56 Two_Pi : constant T := 2.0 * Pi;
57 Half_Pi : constant T := Pi / 2.0;
59 Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two;
60 Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two;
61 Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
63 -----------------------
64 -- Local Subprograms --
65 -----------------------
67 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base;
68 -- Cody/Waite routine, supposedly more precise than the library
69 -- version. Currently only needed for Sinh/Cosh on X86 with the largest
70 -- FP type.
72 function Local_Atan
73 (Y : Float_Type'Base;
74 X : Float_Type'Base := 1.0)
75 return Float_Type'Base;
76 -- Common code for arc tangent after cyele reduction
78 ----------
79 -- "**" --
80 ----------
82 function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
83 A_Right : Float_Type'Base;
84 Int_Part : Integer;
85 Result : Float_Type'Base;
86 R1 : Float_Type'Base;
87 Rest : Float_Type'Base;
89 begin
90 if Left = 0.0
91 and then Right = 0.0
92 then
93 raise Argument_Error;
95 elsif Left < 0.0 then
96 raise Argument_Error;
98 elsif Right = 0.0 then
99 return 1.0;
101 elsif Left = 0.0 then
102 if Right < 0.0 then
103 raise Constraint_Error;
104 else
105 return 0.0;
106 end if;
108 elsif Left = 1.0 then
109 return 1.0;
111 elsif Right = 1.0 then
112 return Left;
114 else
115 begin
116 if Right = 2.0 then
117 return Left * Left;
119 elsif Right = 0.5 then
120 return Sqrt (Left);
122 else
123 A_Right := abs (Right);
125 -- If exponent is larger than one, compute integer exponen-
126 -- tiation if possible, and evaluate fractional part with
127 -- more precision. The relative error is now proportional
128 -- to the fractional part of the exponent only.
130 if A_Right > 1.0
131 and then A_Right < Float_Type'Base (Integer'Last)
132 then
133 Int_Part := Integer (Float_Type'Base'Truncation (A_Right));
134 Result := Left ** Int_Part;
135 Rest := A_Right - Float_Type'Base (Int_Part);
137 -- Compute with two leading bits of the mantissa using
138 -- square roots. Bound to be better than logarithms, and
139 -- easily extended to greater precision.
141 if Rest >= 0.5 then
142 R1 := Sqrt (Left);
143 Result := Result * R1;
144 Rest := Rest - 0.5;
146 if Rest >= 0.25 then
147 Result := Result * Sqrt (R1);
148 Rest := Rest - 0.25;
149 end if;
151 elsif Rest >= 0.25 then
152 Result := Result * Sqrt (Sqrt (Left));
153 Rest := Rest - 0.25;
154 end if;
156 Result := Result *
157 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));
159 if Right >= 0.0 then
160 return Result;
161 else
162 return (1.0 / Result);
163 end if;
164 else
165 return
166 Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
167 end if;
168 end if;
170 exception
171 when others =>
172 raise Constraint_Error;
173 end;
174 end if;
175 end "**";
177 ------------
178 -- Arccos --
179 ------------
181 -- Natural cycle
183 function Arccos (X : Float_Type'Base) return Float_Type'Base is
184 Temp : Float_Type'Base;
186 begin
187 if abs X > 1.0 then
188 raise Argument_Error;
190 elsif abs X < Sqrt_Epsilon then
191 return Pi / 2.0 - X;
193 elsif X = 1.0 then
194 return 0.0;
196 elsif X = -1.0 then
197 return Pi;
198 end if;
200 Temp := Float_Type'Base (Aux.Acos (Double (X)));
202 if Temp < 0.0 then
203 Temp := Pi + Temp;
204 end if;
206 return Temp;
207 end Arccos;
209 -- Arbitrary cycle
211 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
212 Temp : Float_Type'Base;
214 begin
215 if Cycle <= 0.0 then
216 raise Argument_Error;
218 elsif abs X > 1.0 then
219 raise Argument_Error;
221 elsif abs X < Sqrt_Epsilon then
222 return Cycle / 4.0;
224 elsif X = 1.0 then
225 return 0.0;
227 elsif X = -1.0 then
228 return Cycle / 2.0;
229 end if;
231 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);
233 if Temp < 0.0 then
234 Temp := Cycle / 2.0 + Temp;
235 end if;
237 return Temp;
238 end Arccos;
240 -------------
241 -- Arccosh --
242 -------------
244 function Arccosh (X : Float_Type'Base) return Float_Type'Base is
245 begin
246 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or
247 -- the proper approximation for X close to 1 or >> 1.
249 if X < 1.0 then
250 raise Argument_Error;
252 elsif X < 1.0 + Sqrt_Epsilon then
253 return Sqrt (2.0 * (X - 1.0));
255 elsif X > 1.0 / Sqrt_Epsilon then
256 return Log (X) + Log_Two;
258 else
259 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
260 end if;
261 end Arccosh;
263 ------------
264 -- Arccot --
265 ------------
267 -- Natural cycle
269 function Arccot
270 (X : Float_Type'Base;
271 Y : Float_Type'Base := 1.0)
272 return Float_Type'Base
274 begin
275 -- Just reverse arguments
277 return Arctan (Y, X);
278 end Arccot;
280 -- Arbitrary cycle
282 function Arccot
283 (X : Float_Type'Base;
284 Y : Float_Type'Base := 1.0;
285 Cycle : Float_Type'Base)
286 return Float_Type'Base
288 begin
289 -- Just reverse arguments
291 return Arctan (Y, X, Cycle);
292 end Arccot;
294 -------------
295 -- Arccoth --
296 -------------
298 function Arccoth (X : Float_Type'Base) return Float_Type'Base is
299 begin
300 if abs X > 2.0 then
301 return Arctanh (1.0 / X);
303 elsif abs X = 1.0 then
304 raise Constraint_Error;
306 elsif abs X < 1.0 then
307 raise Argument_Error;
309 else
310 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the
311 -- other has error 0 or Epsilon.
313 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
314 end if;
315 end Arccoth;
317 ------------
318 -- Arcsin --
319 ------------
321 -- Natural cycle
323 function Arcsin (X : Float_Type'Base) return Float_Type'Base is
324 begin
325 if abs X > 1.0 then
326 raise Argument_Error;
328 elsif abs X < Sqrt_Epsilon then
329 return X;
331 elsif X = 1.0 then
332 return Pi / 2.0;
334 elsif X = -1.0 then
335 return -Pi / 2.0;
336 end if;
338 return Float_Type'Base (Aux.Asin (Double (X)));
339 end Arcsin;
341 -- Arbitrary cycle
343 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
344 begin
345 if Cycle <= 0.0 then
346 raise Argument_Error;
348 elsif abs X > 1.0 then
349 raise Argument_Error;
351 elsif X = 0.0 then
352 return X;
354 elsif X = 1.0 then
355 return Cycle / 4.0;
357 elsif X = -1.0 then
358 return -Cycle / 4.0;
359 end if;
361 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
362 end Arcsin;
364 -------------
365 -- Arcsinh --
366 -------------
368 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
369 begin
370 if abs X < Sqrt_Epsilon then
371 return X;
373 elsif X > 1.0 / Sqrt_Epsilon then
374 return Log (X) + Log_Two;
376 elsif X < -1.0 / Sqrt_Epsilon then
377 return -(Log (-X) + Log_Two);
379 elsif X < 0.0 then
380 return -Log (abs X + Sqrt (X * X + 1.0));
382 else
383 return Log (X + Sqrt (X * X + 1.0));
384 end if;
385 end Arcsinh;
387 ------------
388 -- Arctan --
389 ------------
391 -- Natural cycle
393 function Arctan
394 (Y : Float_Type'Base;
395 X : Float_Type'Base := 1.0)
396 return Float_Type'Base
398 begin
399 if X = 0.0
400 and then Y = 0.0
401 then
402 raise Argument_Error;
404 elsif Y = 0.0 then
405 if X > 0.0 then
406 return 0.0;
407 else -- X < 0.0
408 return Pi * Float_Type'Copy_Sign (1.0, Y);
409 end if;
411 elsif X = 0.0 then
412 if Y > 0.0 then
413 return Half_Pi;
414 else -- Y < 0.0
415 return -Half_Pi;
416 end if;
418 else
419 return Local_Atan (Y, X);
420 end if;
421 end Arctan;
423 -- Arbitrary cycle
425 function Arctan
426 (Y : Float_Type'Base;
427 X : Float_Type'Base := 1.0;
428 Cycle : Float_Type'Base)
429 return Float_Type'Base
431 begin
432 if Cycle <= 0.0 then
433 raise Argument_Error;
435 elsif X = 0.0
436 and then Y = 0.0
437 then
438 raise Argument_Error;
440 elsif Y = 0.0 then
441 if X > 0.0 then
442 return 0.0;
443 else -- X < 0.0
444 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
445 end if;
447 elsif X = 0.0 then
448 if Y > 0.0 then
449 return Cycle / 4.0;
450 else -- Y < 0.0
451 return -Cycle / 4.0;
452 end if;
454 else
455 return Local_Atan (Y, X) * Cycle / Two_Pi;
456 end if;
457 end Arctan;
459 -------------
460 -- Arctanh --
461 -------------
463 function Arctanh (X : Float_Type'Base) return Float_Type'Base is
464 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;
465 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa;
467 begin
468 -- The naive formula:
470 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X)
472 -- is not well-behaved numerically when X < 0.5 and when X is close
473 -- to one. The following is accurate but probably not optimal.
475 if abs X = 1.0 then
476 raise Constraint_Error;
478 elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then
480 if abs X >= 1.0 then
481 raise Argument_Error;
482 else
484 -- The one case that overflows if put through the method below:
485 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is
486 -- accurate. This simplifies to:
488 return Float_Type'Copy_Sign (
489 Half_Log_Two * Float_Type'Base (Mantissa + 1), X);
490 end if;
492 -- elsif abs X <= 0.5 then
493 -- why is above line commented out ???
495 else
496 -- Use several piecewise linear approximations.
497 -- A is close to X, chosen so 1.0 + A, 1.0 - A, and X - A are exact.
498 -- The two scalings remove the low-order bits of X.
500 A := Float_Type'Base'Scaling (
501 Float_Type'Base (Long_Long_Integer
502 (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa);
504 B := X - A; -- This is exact; abs B <= 2**(-Mantissa).
505 A_Plus_1 := 1.0 + A; -- This is exact.
506 A_From_1 := 1.0 - A; -- Ditto.
507 D := A_Plus_1 * A_From_1; -- 1 - A*A.
509 -- use one term of the series expansion:
510 -- f (x + e) = f(x) + e * f'(x) + ..
512 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is
513 -- A*(B/D)**2 (if a quadratic approximation is ever needed).
515 return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D;
517 -- else
518 -- return 0.5 * Log ((X + 1.0) / (1.0 - X));
519 -- why are above lines commented out ???
520 end if;
521 end Arctanh;
523 ---------
524 -- Cos --
525 ---------
527 -- Natural cycle
529 function Cos (X : Float_Type'Base) return Float_Type'Base is
530 begin
531 if X = 0.0 then
532 return 1.0;
534 elsif abs X < Sqrt_Epsilon then
535 return 1.0;
537 end if;
539 return Float_Type'Base (Aux.Cos (Double (X)));
540 end Cos;
542 -- Arbitrary cycle
544 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
545 begin
546 -- Just reuse the code for Sin. The potential small
547 -- loss of speed is negligible with proper (front-end) inlining.
549 return -Sin (abs X - Cycle * 0.25, Cycle);
550 end Cos;
552 ----------
553 -- Cosh --
554 ----------
556 function Cosh (X : Float_Type'Base) return Float_Type'Base is
557 Lnv : constant Float_Type'Base := 8#0.542714#;
558 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
559 Y : constant Float_Type'Base := abs X;
560 Z : Float_Type'Base;
562 begin
563 if Y < Sqrt_Epsilon then
564 return 1.0;
566 elsif Y > Log_Inverse_Epsilon then
567 Z := Exp_Strict (Y - Lnv);
568 return (Z + V2minus1 * Z);
570 else
571 Z := Exp_Strict (Y);
572 return 0.5 * (Z + 1.0 / Z);
573 end if;
575 end Cosh;
577 ---------
578 -- Cot --
579 ---------
581 -- Natural cycle
583 function Cot (X : Float_Type'Base) return Float_Type'Base is
584 begin
585 if X = 0.0 then
586 raise Constraint_Error;
588 elsif abs X < Sqrt_Epsilon then
589 return 1.0 / X;
590 end if;
592 return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
593 end Cot;
595 -- Arbitrary cycle
597 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
598 T : Float_Type'Base;
600 begin
601 if Cycle <= 0.0 then
602 raise Argument_Error;
603 end if;
605 T := Float_Type'Base'Remainder (X, Cycle);
607 if T = 0.0 or abs T = 0.5 * Cycle then
608 raise Constraint_Error;
610 elsif abs T < Sqrt_Epsilon then
611 return 1.0 / T;
613 elsif abs T = 0.25 * Cycle then
614 return 0.0;
616 else
617 T := T / Cycle * Two_Pi;
618 return Cos (T) / Sin (T);
619 end if;
620 end Cot;
622 ----------
623 -- Coth --
624 ----------
626 function Coth (X : Float_Type'Base) return Float_Type'Base is
627 begin
628 if X = 0.0 then
629 raise Constraint_Error;
631 elsif X < Half_Log_Epsilon then
632 return -1.0;
634 elsif X > -Half_Log_Epsilon then
635 return 1.0;
637 elsif abs X < Sqrt_Epsilon then
638 return 1.0 / X;
639 end if;
641 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
642 end Coth;
644 ---------
645 -- Exp --
646 ---------
648 function Exp (X : Float_Type'Base) return Float_Type'Base is
649 Result : Float_Type'Base;
651 begin
652 if X = 0.0 then
653 return 1.0;
654 end if;
656 Result := Float_Type'Base (Aux.Exp (Double (X)));
658 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
659 -- is False, then we can just leave it as an infinity (and indeed we
660 -- prefer to do so). But if Machine_Overflows is True, then we have
661 -- to raise a Constraint_Error exception as required by the RM.
663 if Float_Type'Machine_Overflows and then not Result'Valid then
664 raise Constraint_Error;
665 end if;
667 return Result;
668 end Exp;
670 ----------------
671 -- Exp_Strict --
672 ----------------
674 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
675 G : Float_Type'Base;
676 Z : Float_Type'Base;
678 P0 : constant := 0.25000_00000_00000_00000;
679 P1 : constant := 0.75753_18015_94227_76666E-2;
680 P2 : constant := 0.31555_19276_56846_46356E-4;
682 Q0 : constant := 0.5;
683 Q1 : constant := 0.56817_30269_85512_21787E-1;
684 Q2 : constant := 0.63121_89437_43985_02557E-3;
685 Q3 : constant := 0.75104_02839_98700_46114E-6;
687 C1 : constant := 8#0.543#;
688 C2 : constant := -2.1219_44400_54690_58277E-4;
689 Le : constant := 1.4426_95040_88896_34074;
691 XN : Float_Type'Base;
692 P, Q, R : Float_Type'Base;
694 begin
695 if X = 0.0 then
696 return 1.0;
697 end if;
699 XN := Float_Type'Base'Rounding (X * Le);
700 G := (X - XN * C1) - XN * C2;
701 Z := G * G;
702 P := G * ((P2 * Z + P1) * Z + P0);
703 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
704 R := 0.5 + P / (Q - P);
706 R := Float_Type'Base'Scaling (R, Integer (XN) + 1);
708 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
709 -- is False, then we can just leave it as an infinity (and indeed we
710 -- prefer to do so). But if Machine_Overflows is True, then we have
711 -- to raise a Constraint_Error exception as required by the RM.
713 if Float_Type'Machine_Overflows and then not R'Valid then
714 raise Constraint_Error;
715 else
716 return R;
717 end if;
719 end Exp_Strict;
721 ----------------
722 -- Local_Atan --
723 ----------------
725 function Local_Atan
726 (Y : Float_Type'Base;
727 X : Float_Type'Base := 1.0)
728 return Float_Type'Base
730 Z : Float_Type'Base;
731 Raw_Atan : Float_Type'Base;
733 begin
734 if abs Y > abs X then
735 Z := abs (X / Y);
736 else
737 Z := abs (Y / X);
738 end if;
740 if Z < Sqrt_Epsilon then
741 Raw_Atan := Z;
743 elsif Z = 1.0 then
744 Raw_Atan := Pi / 4.0;
746 else
747 Raw_Atan := Float_Type'Base (Aux.Atan (Double (Z)));
748 end if;
750 if abs Y > abs X then
751 Raw_Atan := Half_Pi - Raw_Atan;
752 end if;
754 if X > 0.0 then
755 if Y > 0.0 then
756 return Raw_Atan;
757 else -- Y < 0.0
758 return -Raw_Atan;
759 end if;
761 else -- X < 0.0
762 if Y > 0.0 then
763 return Pi - Raw_Atan;
764 else -- Y < 0.0
765 return -(Pi - Raw_Atan);
766 end if;
767 end if;
768 end Local_Atan;
770 ---------
771 -- Log --
772 ---------
774 -- Natural base
776 function Log (X : Float_Type'Base) return Float_Type'Base is
777 begin
778 if X < 0.0 then
779 raise Argument_Error;
781 elsif X = 0.0 then
782 raise Constraint_Error;
784 elsif X = 1.0 then
785 return 0.0;
786 end if;
788 return Float_Type'Base (Aux.Log (Double (X)));
789 end Log;
791 -- Arbitrary base
793 function Log (X, Base : Float_Type'Base) return Float_Type'Base is
794 begin
795 if X < 0.0 then
796 raise Argument_Error;
798 elsif Base <= 0.0 or else Base = 1.0 then
799 raise Argument_Error;
801 elsif X = 0.0 then
802 raise Constraint_Error;
804 elsif X = 1.0 then
805 return 0.0;
806 end if;
808 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
809 end Log;
811 ---------
812 -- Sin --
813 ---------
815 -- Natural cycle
817 function Sin (X : Float_Type'Base) return Float_Type'Base is
818 begin
819 if abs X < Sqrt_Epsilon then
820 return X;
821 end if;
823 return Float_Type'Base (Aux.Sin (Double (X)));
824 end Sin;
826 -- Arbitrary cycle
828 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
829 T : Float_Type'Base;
831 begin
832 if Cycle <= 0.0 then
833 raise Argument_Error;
835 elsif X = 0.0 then
836 -- Is this test really needed on any machine ???
837 return X;
838 end if;
840 T := Float_Type'Base'Remainder (X, Cycle);
842 -- The following two reductions reduce the argument
843 -- to the interval [-0.25 * Cycle, 0.25 * Cycle].
844 -- This reduction is exact and is needed to prevent
845 -- inaccuracy that may result if the sinus function
846 -- a different (more accurate) value of Pi in its
847 -- reduction than is used in the multiplication with Two_Pi.
849 if abs T > 0.25 * Cycle then
850 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
851 end if;
853 -- Could test for 12.0 * abs T = Cycle, and return
854 -- an exact value in those cases. It is not clear that
855 -- this is worth the extra test though.
857 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
858 end Sin;
860 ----------
861 -- Sinh --
862 ----------
864 function Sinh (X : Float_Type'Base) return Float_Type'Base is
865 Lnv : constant Float_Type'Base := 8#0.542714#;
866 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
867 Y : constant Float_Type'Base := abs X;
868 F : constant Float_Type'Base := Y * Y;
869 Z : Float_Type'Base;
871 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;
873 begin
874 if Y < Sqrt_Epsilon then
875 return X;
877 elsif Y > Log_Inverse_Epsilon then
878 Z := Exp_Strict (Y - Lnv);
879 Z := Z + V2minus1 * Z;
881 elsif Y < 1.0 then
883 if Float_Digits_1_6 then
885 -- Use expansion provided by Cody and Waite, p. 226. Note that
886 -- leading term of the polynomial in Q is exactly 1.0.
888 declare
889 P0 : constant := -0.71379_3159E+1;
890 P1 : constant := -0.19033_3399E+0;
891 Q0 : constant := -0.42827_7109E+2;
893 begin
894 Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
895 end;
897 else
898 declare
899 P0 : constant := -0.35181_28343_01771_17881E+6;
900 P1 : constant := -0.11563_52119_68517_68270E+5;
901 P2 : constant := -0.16375_79820_26307_51372E+3;
902 P3 : constant := -0.78966_12741_73570_99479E+0;
903 Q0 : constant := -0.21108_77005_81062_71242E+7;
904 Q1 : constant := 0.36162_72310_94218_36460E+5;
905 Q2 : constant := -0.27773_52311_96507_01667E+3;
907 begin
908 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
909 / (((F + Q2) * F + Q1) * F + Q0);
910 end;
911 end if;
913 else
914 Z := Exp_Strict (Y);
915 Z := 0.5 * (Z - 1.0 / Z);
916 end if;
918 if X > 0.0 then
919 return Z;
920 else
921 return -Z;
922 end if;
923 end Sinh;
925 ----------
926 -- Sqrt --
927 ----------
929 function Sqrt (X : Float_Type'Base) return Float_Type'Base is
930 begin
931 if X < 0.0 then
932 raise Argument_Error;
934 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
936 elsif X = 0.0 then
937 return X;
939 end if;
941 return Float_Type'Base (Aux.Sqrt (Double (X)));
942 end Sqrt;
944 ---------
945 -- Tan --
946 ---------
948 -- Natural cycle
950 function Tan (X : Float_Type'Base) return Float_Type'Base is
951 begin
952 if abs X < Sqrt_Epsilon then
953 return X;
955 elsif abs X = Pi / 2.0 then
956 raise Constraint_Error;
957 end if;
959 return Float_Type'Base (Aux.Tan (Double (X)));
960 end Tan;
962 -- Arbitrary cycle
964 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
965 T : Float_Type'Base;
967 begin
968 if Cycle <= 0.0 then
969 raise Argument_Error;
971 elsif X = 0.0 then
972 return X;
973 end if;
975 T := Float_Type'Base'Remainder (X, Cycle);
977 if abs T = 0.25 * Cycle then
978 raise Constraint_Error;
980 elsif abs T = 0.5 * Cycle then
981 return 0.0;
983 else
984 T := T / Cycle * Two_Pi;
985 return Sin (T) / Cos (T);
986 end if;
988 end Tan;
990 ----------
991 -- Tanh --
992 ----------
994 function Tanh (X : Float_Type'Base) return Float_Type'Base is
995 P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
996 P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
997 P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;
999 Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4;
1000 Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4;
1001 Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3;
1002 Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1;
1004 Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
1006 P, Q, R : Float_Type'Base;
1007 Y : constant Float_Type'Base := abs X;
1008 G : constant Float_Type'Base := Y * Y;
1010 Float_Type_Digits_15_Or_More : constant Boolean :=
1011 Float_Type'Digits > 14;
1013 begin
1014 if X < Half_Log_Epsilon then
1015 return -1.0;
1017 elsif X > -Half_Log_Epsilon then
1018 return 1.0;
1020 elsif Y < Sqrt_Epsilon then
1021 return X;
1023 elsif Y < Half_Ln3
1024 and then Float_Type_Digits_15_Or_More
1025 then
1026 P := (P2 * G + P1) * G + P0;
1027 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
1028 R := G * (P / Q);
1029 return X + X * R;
1031 else
1032 return Float_Type'Base (Aux.Tanh (Double (X)));
1033 end if;
1034 end Tanh;
1036 end Ada.Numerics.Generic_Elementary_Functions;