File : a-ngelfu.adb
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-2015, 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 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. --
17 -- --
18 -- --
19 -- --
20 -- --
21 -- --
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/>. --
26 -- --
27 -- GNAT was originally developed by the GNAT team at New York University. --
28 -- Extensive contributions were provided by Ada Core Technologies Inc. --
29 -- --
30 ------------------------------------------------------------------------------
31
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 or to meet Annex G strict mode requirements.
35
36 -- Uses functions sqrt, exp, log, pow, sin, asin, cos, acos, tan, atan, sinh,
37 -- cosh, tanh from C library via math.h
38
39 with Ada.Numerics.Aux;
40
41 package body Ada.Numerics.Generic_Elementary_Functions is
42
43 use type Ada.Numerics.Aux.Double;
44
45 Sqrt_Two : constant := 1.41421_35623_73095_04880_16887_24209_69807_85696;
46 Log_Two : constant := 0.69314_71805_59945_30941_72321_21458_17656_80755;
47
48 Half_Log_Two : constant := Log_Two / 2;
49
50 subtype T is Float_Type'Base;
51 subtype Double is Aux.Double;
52
53 Two_Pi : constant T := 2.0 * Pi;
54 Half_Pi : constant T := Pi / 2.0;
55
56 Half_Log_Epsilon : constant T := T (1 - T'Model_Mantissa) * Half_Log_Two;
57 Log_Inverse_Epsilon : constant T := T (T'Model_Mantissa - 1) * Log_Two;
58 Sqrt_Epsilon : constant T := Sqrt_Two ** (1 - T'Model_Mantissa);
59
60 -----------------------
61 -- Local Subprograms --
62 -----------------------
63
64 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base;
65 -- Cody/Waite routine, supposedly more precise than the library version.
66 -- Currently only needed for Sinh/Cosh on X86 with the largest FP type.
67
68 function Local_Atan
69 (Y : Float_Type'Base;
70 X : Float_Type'Base := 1.0) return Float_Type'Base;
71 -- Common code for arc tangent after cycle reduction
72
73 ----------
74 -- "**" --
75 ----------
76
77 function "**" (Left, Right : Float_Type'Base) return Float_Type'Base is
78 A_Right : Float_Type'Base;
79 Int_Part : Integer;
80 Result : Float_Type'Base;
81 R1 : Float_Type'Base;
82 Rest : Float_Type'Base;
83
84 begin
85 if Left = 0.0
86 and then Right = 0.0
87 then
88 raise Argument_Error;
89
90 elsif Left < 0.0 then
91 raise Argument_Error;
92
93 elsif Right = 0.0 then
94 return 1.0;
95
96 elsif Left = 0.0 then
97 if Right < 0.0 then
98 raise Constraint_Error;
99 else
100 return 0.0;
101 end if;
102
103 elsif Left = 1.0 then
104 return 1.0;
105
106 elsif Right = 1.0 then
107 return Left;
108
109 else
110 begin
111 if Right = 2.0 then
112 return Left * Left;
113
114 elsif Right = 0.5 then
115 return Sqrt (Left);
116
117 else
118 A_Right := abs (Right);
119
120 -- If exponent is larger than one, compute integer exponen-
121 -- tiation if possible, and evaluate fractional part with more
122 -- precision. The relative error is now proportional to the
123 -- fractional part of the exponent only.
124
125 if A_Right > 1.0
126 and then A_Right < Float_Type'Base (Integer'Last)
127 then
128 Int_Part := Integer (Float_Type'Base'Truncation (A_Right));
129 Result := Left ** Int_Part;
130 Rest := A_Right - Float_Type'Base (Int_Part);
131
132 -- Compute with two leading bits of the mantissa using
133 -- square roots. Bound to be better than logarithms, and
134 -- easily extended to greater precision.
135
136 if Rest >= 0.5 then
137 R1 := Sqrt (Left);
138 Result := Result * R1;
139 Rest := Rest - 0.5;
140
141 if Rest >= 0.25 then
142 Result := Result * Sqrt (R1);
143 Rest := Rest - 0.25;
144 end if;
145
146 elsif Rest >= 0.25 then
147 Result := Result * Sqrt (Sqrt (Left));
148 Rest := Rest - 0.25;
149 end if;
150
151 Result := Result *
152 Float_Type'Base (Aux.Pow (Double (Left), Double (Rest)));
153
154 if Right >= 0.0 then
155 return Result;
156 else
157 return (1.0 / Result);
158 end if;
159 else
160 return
161 Float_Type'Base (Aux.Pow (Double (Left), Double (Right)));
162 end if;
163 end if;
164
165 exception
166 when others =>
167 raise Constraint_Error;
168 end;
169 end if;
170 end "**";
171
172 ------------
173 -- Arccos --
174 ------------
175
176 -- Natural cycle
177
178 function Arccos (X : Float_Type'Base) return Float_Type'Base is
179 Temp : Float_Type'Base;
180
181 begin
182 if abs X > 1.0 then
183 raise Argument_Error;
184
185 elsif abs X < Sqrt_Epsilon then
186 return Pi / 2.0 - X;
187
188 elsif X = 1.0 then
189 return 0.0;
190
191 elsif X = -1.0 then
192 return Pi;
193 end if;
194
195 Temp := Float_Type'Base (Aux.Acos (Double (X)));
196
197 if Temp < 0.0 then
198 Temp := Pi + Temp;
199 end if;
200
201 return Temp;
202 end Arccos;
203
204 -- Arbitrary cycle
205
206 function Arccos (X, Cycle : Float_Type'Base) return Float_Type'Base is
207 Temp : Float_Type'Base;
208
209 begin
210 if Cycle <= 0.0 then
211 raise Argument_Error;
212
213 elsif abs X > 1.0 then
214 raise Argument_Error;
215
216 elsif abs X < Sqrt_Epsilon then
217 return Cycle / 4.0;
218
219 elsif X = 1.0 then
220 return 0.0;
221
222 elsif X = -1.0 then
223 return Cycle / 2.0;
224 end if;
225
226 Temp := Arctan (Sqrt ((1.0 - X) * (1.0 + X)) / X, 1.0, Cycle);
227
228 if Temp < 0.0 then
229 Temp := Cycle / 2.0 + Temp;
230 end if;
231
232 return Temp;
233 end Arccos;
234
235 -------------
236 -- Arccosh --
237 -------------
238
239 function Arccosh (X : Float_Type'Base) return Float_Type'Base is
240 begin
241 -- Return positive branch of Log (X - Sqrt (X * X - 1.0)), or the proper
242 -- approximation for X close to 1 or >> 1.
243
244 if X < 1.0 then
245 raise Argument_Error;
246
247 elsif X < 1.0 + Sqrt_Epsilon then
248 return Sqrt (2.0 * (X - 1.0));
249
250 elsif X > 1.0 / Sqrt_Epsilon then
251 return Log (X) + Log_Two;
252
253 else
254 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
255 end if;
256 end Arccosh;
257
258 ------------
259 -- Arccot --
260 ------------
261
262 -- Natural cycle
263
264 function Arccot
265 (X : Float_Type'Base;
266 Y : Float_Type'Base := 1.0)
267 return Float_Type'Base
268 is
269 begin
270 -- Just reverse arguments
271
272 return Arctan (Y, X);
273 end Arccot;
274
275 -- Arbitrary cycle
276
277 function Arccot
278 (X : Float_Type'Base;
279 Y : Float_Type'Base := 1.0;
280 Cycle : Float_Type'Base)
281 return Float_Type'Base
282 is
283 begin
284 -- Just reverse arguments
285
286 return Arctan (Y, X, Cycle);
287 end Arccot;
288
289 -------------
290 -- Arccoth --
291 -------------
292
293 function Arccoth (X : Float_Type'Base) return Float_Type'Base is
294 begin
295 if abs X > 2.0 then
296 return Arctanh (1.0 / X);
297
298 elsif abs X = 1.0 then
299 raise Constraint_Error;
300
301 elsif abs X < 1.0 then
302 raise Argument_Error;
303
304 else
305 -- 1.0 < abs X <= 2.0. One of X + 1.0 and X - 1.0 is exact, the other
306 -- has error 0 or Epsilon.
307
308 return 0.5 * (Log (abs (X + 1.0)) - Log (abs (X - 1.0)));
309 end if;
310 end Arccoth;
311
312 ------------
313 -- Arcsin --
314 ------------
315
316 -- Natural cycle
317
318 function Arcsin (X : Float_Type'Base) return Float_Type'Base is
319 begin
320 if abs X > 1.0 then
321 raise Argument_Error;
322
323 elsif abs X < Sqrt_Epsilon then
324 return X;
325
326 elsif X = 1.0 then
327 return Pi / 2.0;
328
329 elsif X = -1.0 then
330 return -(Pi / 2.0);
331 end if;
332
333 return Float_Type'Base (Aux.Asin (Double (X)));
334 end Arcsin;
335
336 -- Arbitrary cycle
337
338 function Arcsin (X, Cycle : Float_Type'Base) return Float_Type'Base is
339 begin
340 if Cycle <= 0.0 then
341 raise Argument_Error;
342
343 elsif abs X > 1.0 then
344 raise Argument_Error;
345
346 elsif X = 0.0 then
347 return X;
348
349 elsif X = 1.0 then
350 return Cycle / 4.0;
351
352 elsif X = -1.0 then
353 return -(Cycle / 4.0);
354 end if;
355
356 return Arctan (X / Sqrt ((1.0 - X) * (1.0 + X)), 1.0, Cycle);
357 end Arcsin;
358
359 -------------
360 -- Arcsinh --
361 -------------
362
363 function Arcsinh (X : Float_Type'Base) return Float_Type'Base is
364 begin
365 if abs X < Sqrt_Epsilon then
366 return X;
367
368 elsif X > 1.0 / Sqrt_Epsilon then
369 return Log (X) + Log_Two;
370
371 elsif X < -(1.0 / Sqrt_Epsilon) then
372 return -(Log (-X) + Log_Two);
373
374 elsif X < 0.0 then
375 return -Log (abs X + Sqrt (X * X + 1.0));
376
377 else
378 return Log (X + Sqrt (X * X + 1.0));
379 end if;
380 end Arcsinh;
381
382 ------------
383 -- Arctan --
384 ------------
385
386 -- Natural cycle
387
388 function Arctan
389 (Y : Float_Type'Base;
390 X : Float_Type'Base := 1.0)
391 return Float_Type'Base
392 is
393 begin
394 if X = 0.0 and then Y = 0.0 then
395 raise Argument_Error;
396
397 elsif Y = 0.0 then
398 if X > 0.0 then
399 return 0.0;
400 else -- X < 0.0
401 return Pi * Float_Type'Copy_Sign (1.0, Y);
402 end if;
403
404 elsif X = 0.0 then
405 return Float_Type'Copy_Sign (Half_Pi, Y);
406
407 else
408 return Local_Atan (Y, X);
409 end if;
410 end Arctan;
411
412 -- Arbitrary cycle
413
414 function Arctan
415 (Y : Float_Type'Base;
416 X : Float_Type'Base := 1.0;
417 Cycle : Float_Type'Base)
418 return Float_Type'Base
419 is
420 begin
421 if Cycle <= 0.0 then
422 raise Argument_Error;
423
424 elsif X = 0.0 and then Y = 0.0 then
425 raise Argument_Error;
426
427 elsif Y = 0.0 then
428 if X > 0.0 then
429 return 0.0;
430 else -- X < 0.0
431 return Cycle / 2.0 * Float_Type'Copy_Sign (1.0, Y);
432 end if;
433
434 elsif X = 0.0 then
435 return Float_Type'Copy_Sign (Cycle / 4.0, Y);
436
437 else
438 return Local_Atan (Y, X) * Cycle / Two_Pi;
439 end if;
440 end Arctan;
441
442 -------------
443 -- Arctanh --
444 -------------
445
446 function Arctanh (X : Float_Type'Base) return Float_Type'Base is
447 A, B, D, A_Plus_1, A_From_1 : Float_Type'Base;
448
449 Mantissa : constant Integer := Float_Type'Base'Machine_Mantissa;
450
451 begin
452 -- The naive formula:
453
454 -- Arctanh (X) := (1/2) * Log (1 + X) / (1 - X)
455
456 -- is not well-behaved numerically when X < 0.5 and when X is close
457 -- to one. The following is accurate but probably not optimal.
458
459 if abs X = 1.0 then
460 raise Constraint_Error;
461
462 elsif abs X >= 1.0 - 2.0 ** (-Mantissa) then
463
464 if abs X >= 1.0 then
465 raise Argument_Error;
466 else
467
468 -- The one case that overflows if put through the method below:
469 -- abs X = 1.0 - Epsilon. In this case (1/2) log (2/Epsilon) is
470 -- accurate. This simplifies to:
471
472 return Float_Type'Copy_Sign (
473 Half_Log_Two * Float_Type'Base (Mantissa + 1), X);
474 end if;
475
476 -- elsif abs X <= 0.5 then
477 -- why is above line commented out ???
478
479 else
480 -- Use several piecewise linear approximations. A is close to X,
481 -- chosen so 1.0 + A, 1.0 - A, and X - A are exact. The two scalings
482 -- remove the low-order bits of X.
483
484 A := Float_Type'Base'Scaling (
485 Float_Type'Base (Long_Long_Integer
486 (Float_Type'Base'Scaling (X, Mantissa - 1))), 1 - Mantissa);
487
488 B := X - A; -- This is exact; abs B <= 2**(-Mantissa).
489 A_Plus_1 := 1.0 + A; -- This is exact.
490 A_From_1 := 1.0 - A; -- Ditto.
491 D := A_Plus_1 * A_From_1; -- 1 - A*A.
492
493 -- use one term of the series expansion:
494
495 -- f (x + e) = f(x) + e * f'(x) + ..
496
497 -- The derivative of Arctanh at A is 1/(1-A*A). Next term is
498 -- A*(B/D)**2 (if a quadratic approximation is ever needed).
499
500 return 0.5 * (Log (A_Plus_1) - Log (A_From_1)) + B / D;
501 end if;
502 end Arctanh;
503
504 ---------
505 -- Cos --
506 ---------
507
508 -- Natural cycle
509
510 function Cos (X : Float_Type'Base) return Float_Type'Base is
511 begin
512 if abs X < Sqrt_Epsilon then
513 return 1.0;
514 end if;
515
516 return Float_Type'Base (Aux.Cos (Double (X)));
517 end Cos;
518
519 -- Arbitrary cycle
520
521 function Cos (X, Cycle : Float_Type'Base) return Float_Type'Base is
522 begin
523 -- Just reuse the code for Sin. The potential small loss of speed is
524 -- negligible with proper (front-end) inlining.
525
526 return -Sin (abs X - Cycle * 0.25, Cycle);
527 end Cos;
528
529 ----------
530 -- Cosh --
531 ----------
532
533 function Cosh (X : Float_Type'Base) return Float_Type'Base is
534 Lnv : constant Float_Type'Base := 8#0.542714#;
535 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
536 Y : constant Float_Type'Base := abs X;
537 Z : Float_Type'Base;
538
539 begin
540 if Y < Sqrt_Epsilon then
541 return 1.0;
542
543 elsif Y > Log_Inverse_Epsilon then
544 Z := Exp_Strict (Y - Lnv);
545 return (Z + V2minus1 * Z);
546
547 else
548 Z := Exp_Strict (Y);
549 return 0.5 * (Z + 1.0 / Z);
550 end if;
551
552 end Cosh;
553
554 ---------
555 -- Cot --
556 ---------
557
558 -- Natural cycle
559
560 function Cot (X : Float_Type'Base) return Float_Type'Base is
561 begin
562 if X = 0.0 then
563 raise Constraint_Error;
564
565 elsif abs X < Sqrt_Epsilon then
566 return 1.0 / X;
567 end if;
568
569 return 1.0 / Float_Type'Base (Aux.Tan (Double (X)));
570 end Cot;
571
572 -- Arbitrary cycle
573
574 function Cot (X, Cycle : Float_Type'Base) return Float_Type'Base is
575 T : Float_Type'Base;
576
577 begin
578 if Cycle <= 0.0 then
579 raise Argument_Error;
580 end if;
581
582 T := Float_Type'Base'Remainder (X, Cycle);
583
584 if T = 0.0 or else abs T = 0.5 * Cycle then
585 raise Constraint_Error;
586
587 elsif abs T < Sqrt_Epsilon then
588 return 1.0 / T;
589
590 elsif abs T = 0.25 * Cycle then
591 return 0.0;
592
593 else
594 T := T / Cycle * Two_Pi;
595 return Cos (T) / Sin (T);
596 end if;
597 end Cot;
598
599 ----------
600 -- Coth --
601 ----------
602
603 function Coth (X : Float_Type'Base) return Float_Type'Base is
604 begin
605 if X = 0.0 then
606 raise Constraint_Error;
607
608 elsif X < Half_Log_Epsilon then
609 return -1.0;
610
611 elsif X > -Half_Log_Epsilon then
612 return 1.0;
613
614 elsif abs X < Sqrt_Epsilon then
615 return 1.0 / X;
616 end if;
617
618 return 1.0 / Float_Type'Base (Aux.Tanh (Double (X)));
619 end Coth;
620
621 ---------
622 -- Exp --
623 ---------
624
625 function Exp (X : Float_Type'Base) return Float_Type'Base is
626 Result : Float_Type'Base;
627
628 begin
629 if X = 0.0 then
630 return 1.0;
631 end if;
632
633 Result := Float_Type'Base (Aux.Exp (Double (X)));
634
635 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
636 -- is False, then we can just leave it as an infinity (and indeed we
637 -- prefer to do so). But if Machine_Overflows is True, then we have
638 -- to raise a Constraint_Error exception as required by the RM.
639
640 if Float_Type'Machine_Overflows and then not Result'Valid then
641 raise Constraint_Error;
642 end if;
643
644 return Result;
645 end Exp;
646
647 ----------------
648 -- Exp_Strict --
649 ----------------
650
651 function Exp_Strict (X : Float_Type'Base) return Float_Type'Base is
652 G : Float_Type'Base;
653 Z : Float_Type'Base;
654
655 P0 : constant := 0.25000_00000_00000_00000;
656 P1 : constant := 0.75753_18015_94227_76666E-2;
657 P2 : constant := 0.31555_19276_56846_46356E-4;
658
659 Q0 : constant := 0.5;
660 Q1 : constant := 0.56817_30269_85512_21787E-1;
661 Q2 : constant := 0.63121_89437_43985_02557E-3;
662 Q3 : constant := 0.75104_02839_98700_46114E-6;
663
664 C1 : constant := 8#0.543#;
665 C2 : constant := -2.1219_44400_54690_58277E-4;
666 Le : constant := 1.4426_95040_88896_34074;
667
668 XN : Float_Type'Base;
669 P, Q, R : Float_Type'Base;
670
671 begin
672 if X = 0.0 then
673 return 1.0;
674 end if;
675
676 XN := Float_Type'Base'Rounding (X * Le);
677 G := (X - XN * C1) - XN * C2;
678 Z := G * G;
679 P := G * ((P2 * Z + P1) * Z + P0);
680 Q := ((Q3 * Z + Q2) * Z + Q1) * Z + Q0;
681 R := 0.5 + P / (Q - P);
682
683 R := Float_Type'Base'Scaling (R, Integer (XN) + 1);
684
685 -- Deal with case of Exp returning IEEE infinity. If Machine_Overflows
686 -- is False, then we can just leave it as an infinity (and indeed we
687 -- prefer to do so). But if Machine_Overflows is True, then we have to
688 -- raise a Constraint_Error exception as required by the RM.
689
690 if Float_Type'Machine_Overflows and then not R'Valid then
691 raise Constraint_Error;
692 else
693 return R;
694 end if;
695
696 end Exp_Strict;
697
698 ----------------
699 -- Local_Atan --
700 ----------------
701
702 function Local_Atan
703 (Y : Float_Type'Base;
704 X : Float_Type'Base := 1.0) return Float_Type'Base
705 is
706 Z : Float_Type'Base;
707 Raw_Atan : Float_Type'Base;
708
709 begin
710 Z := (if abs Y > abs X then abs (X / Y) else abs (Y / X));
711
712 Raw_Atan :=
713 (if Z < Sqrt_Epsilon then Z
714 elsif Z = 1.0 then Pi / 4.0
715 else Float_Type'Base (Aux.Atan (Double (Z))));
716
717 if abs Y > abs X then
718 Raw_Atan := Half_Pi - Raw_Atan;
719 end if;
720
721 if X > 0.0 then
722 return Float_Type'Copy_Sign (Raw_Atan, Y);
723 else
724 return Float_Type'Copy_Sign (Pi - Raw_Atan, Y);
725 end if;
726 end Local_Atan;
727
728 ---------
729 -- Log --
730 ---------
731
732 -- Natural base
733
734 function Log (X : Float_Type'Base) return Float_Type'Base is
735 begin
736 if X < 0.0 then
737 raise Argument_Error;
738
739 elsif X = 0.0 then
740 raise Constraint_Error;
741
742 elsif X = 1.0 then
743 return 0.0;
744 end if;
745
746 return Float_Type'Base (Aux.Log (Double (X)));
747 end Log;
748
749 -- Arbitrary base
750
751 function Log (X, Base : Float_Type'Base) return Float_Type'Base is
752 begin
753 if X < 0.0 then
754 raise Argument_Error;
755
756 elsif Base <= 0.0 or else Base = 1.0 then
757 raise Argument_Error;
758
759 elsif X = 0.0 then
760 raise Constraint_Error;
761
762 elsif X = 1.0 then
763 return 0.0;
764 end if;
765
766 return Float_Type'Base (Aux.Log (Double (X)) / Aux.Log (Double (Base)));
767 end Log;
768
769 ---------
770 -- Sin --
771 ---------
772
773 -- Natural cycle
774
775 function Sin (X : Float_Type'Base) return Float_Type'Base is
776 begin
777 if abs X < Sqrt_Epsilon then
778 return X;
779 end if;
780
781 return Float_Type'Base (Aux.Sin (Double (X)));
782 end Sin;
783
784 -- Arbitrary cycle
785
786 function Sin (X, Cycle : Float_Type'Base) return Float_Type'Base is
787 T : Float_Type'Base;
788
789 begin
790 if Cycle <= 0.0 then
791 raise Argument_Error;
792
793 -- If X is zero, return it as the result, preserving the argument sign.
794 -- Is this test really needed on any machine ???
795
796 elsif X = 0.0 then
797 return X;
798 end if;
799
800 T := Float_Type'Base'Remainder (X, Cycle);
801
802 -- The following two reductions reduce the argument to the interval
803 -- [-0.25 * Cycle, 0.25 * Cycle]. This reduction is exact and is needed
804 -- to prevent inaccuracy that may result if the sine function uses a
805 -- different (more accurate) value of Pi in its reduction than is used
806 -- in the multiplication with Two_Pi.
807
808 if abs T > 0.25 * Cycle then
809 T := 0.5 * Float_Type'Copy_Sign (Cycle, T) - T;
810 end if;
811
812 -- Could test for 12.0 * abs T = Cycle, and return an exact value in
813 -- those cases. It is not clear this is worth the extra test though.
814
815 return Float_Type'Base (Aux.Sin (Double (T / Cycle * Two_Pi)));
816 end Sin;
817
818 ----------
819 -- Sinh --
820 ----------
821
822 function Sinh (X : Float_Type'Base) return Float_Type'Base is
823 Lnv : constant Float_Type'Base := 8#0.542714#;
824 V2minus1 : constant Float_Type'Base := 0.13830_27787_96019_02638E-4;
825 Y : constant Float_Type'Base := abs X;
826 F : constant Float_Type'Base := Y * Y;
827 Z : Float_Type'Base;
828
829 Float_Digits_1_6 : constant Boolean := Float_Type'Digits < 7;
830
831 begin
832 if Y < Sqrt_Epsilon then
833 return X;
834
835 elsif Y > Log_Inverse_Epsilon then
836 Z := Exp_Strict (Y - Lnv);
837 Z := Z + V2minus1 * Z;
838
839 elsif Y < 1.0 then
840
841 if Float_Digits_1_6 then
842
843 -- Use expansion provided by Cody and Waite, p. 226. Note that
844 -- leading term of the polynomial in Q is exactly 1.0.
845
846 declare
847 P0 : constant := -0.71379_3159E+1;
848 P1 : constant := -0.19033_3399E+0;
849 Q0 : constant := -0.42827_7109E+2;
850
851 begin
852 Z := Y + Y * F * (P1 * F + P0) / (F + Q0);
853 end;
854
855 else
856 declare
857 P0 : constant := -0.35181_28343_01771_17881E+6;
858 P1 : constant := -0.11563_52119_68517_68270E+5;
859 P2 : constant := -0.16375_79820_26307_51372E+3;
860 P3 : constant := -0.78966_12741_73570_99479E+0;
861 Q0 : constant := -0.21108_77005_81062_71242E+7;
862 Q1 : constant := 0.36162_72310_94218_36460E+5;
863 Q2 : constant := -0.27773_52311_96507_01667E+3;
864
865 begin
866 Z := Y + Y * F * (((P3 * F + P2) * F + P1) * F + P0)
867 / (((F + Q2) * F + Q1) * F + Q0);
868 end;
869 end if;
870
871 else
872 Z := Exp_Strict (Y);
873 Z := 0.5 * (Z - 1.0 / Z);
874 end if;
875
876 if X > 0.0 then
877 return Z;
878 else
879 return -Z;
880 end if;
881 end Sinh;
882
883 ----------
884 -- Sqrt --
885 ----------
886
887 function Sqrt (X : Float_Type'Base) return Float_Type'Base is
888 begin
889 if X < 0.0 then
890 raise Argument_Error;
891
892 -- Special case Sqrt (0.0) to preserve possible minus sign per IEEE
893
894 elsif X = 0.0 then
895 return X;
896 end if;
897
898 return Float_Type'Base (Aux.Sqrt (Double (X)));
899 end Sqrt;
900
901 ---------
902 -- Tan --
903 ---------
904
905 -- Natural cycle
906
907 function Tan (X : Float_Type'Base) return Float_Type'Base is
908 begin
909 if abs X < Sqrt_Epsilon then
910 return X;
911 end if;
912
913 -- Note: if X is exactly pi/2, then we should raise an exception, since
914 -- the result would overflow. But for all floating-point formats we deal
915 -- with, it is impossible for X to be exactly pi/2, and the result is
916 -- always in range.
917
918 return Float_Type'Base (Aux.Tan (Double (X)));
919 end Tan;
920
921 -- Arbitrary cycle
922
923 function Tan (X, Cycle : Float_Type'Base) return Float_Type'Base is
924 T : Float_Type'Base;
925
926 begin
927 if Cycle <= 0.0 then
928 raise Argument_Error;
929
930 elsif X = 0.0 then
931 return X;
932 end if;
933
934 T := Float_Type'Base'Remainder (X, Cycle);
935
936 if abs T = 0.25 * Cycle then
937 raise Constraint_Error;
938
939 elsif abs T = 0.5 * Cycle then
940 return 0.0;
941
942 else
943 T := T / Cycle * Two_Pi;
944 return Sin (T) / Cos (T);
945 end if;
946
947 end Tan;
948
949 ----------
950 -- Tanh --
951 ----------
952
953 function Tanh (X : Float_Type'Base) return Float_Type'Base is
954 P0 : constant Float_Type'Base := -0.16134_11902_39962_28053E+4;
955 P1 : constant Float_Type'Base := -0.99225_92967_22360_83313E+2;
956 P2 : constant Float_Type'Base := -0.96437_49277_72254_69787E+0;
957
958 Q0 : constant Float_Type'Base := 0.48402_35707_19886_88686E+4;
959 Q1 : constant Float_Type'Base := 0.22337_72071_89623_12926E+4;
960 Q2 : constant Float_Type'Base := 0.11274_47438_05349_49335E+3;
961 Q3 : constant Float_Type'Base := 0.10000_00000_00000_00000E+1;
962
963 Half_Ln3 : constant Float_Type'Base := 0.54930_61443_34054_84570;
964
965 P, Q, R : Float_Type'Base;
966 Y : constant Float_Type'Base := abs X;
967 G : constant Float_Type'Base := Y * Y;
968
969 Float_Type_Digits_15_Or_More : constant Boolean :=
970 Float_Type'Digits > 14;
971
972 begin
973 if X < Half_Log_Epsilon then
974 return -1.0;
975
976 elsif X > -Half_Log_Epsilon then
977 return 1.0;
978
979 elsif Y < Sqrt_Epsilon then
980 return X;
981
982 elsif Y < Half_Ln3
983 and then Float_Type_Digits_15_Or_More
984 then
985 P := (P2 * G + P1) * G + P0;
986 Q := ((Q3 * G + Q2) * G + Q1) * G + Q0;
987 R := G * (P / Q);
988 return X + X * R;
989
990 else
991 return Float_Type'Base (Aux.Tanh (Double (X)));
992 end if;
993 end Tanh;
994
995 end Ada.Numerics.Generic_Elementary_Functions;