File : s-libsin-ada.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- S Y S T E M . L I B M _ S I N G L E --
6 -- --
7 -- B o d y --
8 -- --
9 -- Copyright (C) 2014-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 is the Ada Cert Math specific version of s-libsin.adb
33
34 -- When Cody and Waite implementation is cited, it refers to the
35 -- Software Manual for the Elementary Functions by William J. Cody, Jr.
36 -- and William Waite, published by Prentice-Hall Series in Computational
37 -- Mathematics. Copyright 1980. ISBN 0-13-822064-6.
38
39 -- When Hart implementation is cited, it refers to
40 -- "Computer Approximations" by John F. Hart, published by Krieger.
41 -- Copyright 1968, Reprinted 1978 w/ corrections. ISBN 0-88275-642-7.
42
43 with Ada.Numerics; use Ada.Numerics;
44 with System.Libm; use System.Libm;
45 with System.Libm_Single.Squareroot;
46
47 package body System.Libm_Single is
48 subtype F is Float;
49
50 pragma Assert (F'Machine_Radix = 2);
51 pragma Assert (F'Machine_Mantissa = 24);
52
53 F_HM : constant Integer := Float'Machine_Mantissa / 2;
54
55 type Float_Table is array (Positive range <>) of Float;
56
57 -- A1 (i) = Float (2**((1-i)/16))
58
59 A1_Tab_F : constant Float_Table :=
60 (1.0,
61 Root16_Half,
62 Root16_Half**2,
63 Root16_Half**3,
64 Root16_Half**4,
65 Root16_Half**5,
66 Root16_Half**6,
67 Root16_Half**7,
68 Root16_Half**8,
69 Root16_Half**9,
70 Root16_Half**10,
71 Root16_Half**11,
72 Root16_Half**12,
73 Root16_Half**13,
74 Root16_Half**14,
75 Root16_Half**15,
76 0.5);
77
78 -- A2 (i) = 2**((1-2i)/16) - A1(2i)
79
80 A2_Tab_F : constant Float_Table :=
81 (Root16_Half - Float'Machine (Root16_Half),
82 Root16_Half**3 - Float'Machine (Root16_Half**3),
83 Root16_Half**5 - Float'Machine (Root16_Half**5),
84 Root16_Half**7 - Float'Machine (Root16_Half**7),
85 Root16_Half**9 - Float'Machine (Root16_Half**9),
86 Root16_Half**11 - Float'Machine (Root16_Half**11),
87 Root16_Half**13 - Float'Machine (Root16_Half**13),
88 Root16_Half**15 - Float'Machine (Root16_Half**15));
89
90 Sqrt_Epsilon_F : constant Float :=
91 Sqrt_2 ** (1 - Float'Machine_Mantissa);
92
93 -- Intermediary functions
94
95 function Reconstruct_Pow
96 (Z : Float;
97 P : Integer;
98 M : Integer) return Float;
99
100 procedure Reduce_044
101 (X : Float;
102 Reduced_X : out Float;
103 P : out Integer)
104 with Post => abs (Reduced_X) < 0.044;
105
106 function Reduce_1_16 (X : Float) return Float
107 with Post => abs (X - Reduce_1_16'Result) <= 0.0625;
108
109 function Reduce_1_16 (X : Float) return Float is
110 (F'Machine_Rounding (X * 16.0) * (1.0 / 16.0));
111
112 package Float_Approximations is
113 new Generic_Approximations (Float, Mantissa => 24);
114
115 use Float_Approximations;
116
117 -- The following function reduces a positive X into the range
118 -- -ln (2) / 2 .. ln (2) / 2
119
120 -- It returns a reduced X and an integer N such that:
121 -- X'Old = X'New + N * Log (2)
122
123 -- It is used by Exp function
124
125 -- The result should be correctly rounded
126
127 procedure Reduce_Ln_2 (X : in out Float; N : out Integer)
128 with Pre => abs (X) <= Float'Ceiling
129 (Float'Pred (88.72283_90520_7) * Inv_Ln_2);
130 -- @llr Reduce_Ln_2 Float
131 -- The following is postcondition doesn't hold. Suspicious "=" ???
132 -- Post => abs (X) <= Ln_2 / 2.0 and
133 -- X'Old = X + Float (N) * Ln_2;
134
135 -- The reduction is used by the Sin, Cos and Tan functions.
136
137 procedure Reduce_Half_Pi (X : in out Float; Q : out Quadrant)
138 with Pre => X >= 0.0,
139 Post => abs (X) <= Max_Red_Trig_Arg;
140 -- @llr Reduce_Half_Pi Float
141 -- The following functions reduce a positive X into the range
142 -- -(Pi/4 + E) .. Pi/4 + E, with E a small fraction of Pi.
143 --
144 -- The reason the normalization is not strict is that the computation of
145 -- the number of times to subtract half Pi is not exact. The rounding
146 -- error is worst for large arguments, where the number of bits behind
147 -- the radix point reduces to half the mantissa bits.
148
149 -- While it would be possible to correct for this, the polynomial
150 -- approximations work well for values slightly outside the -Pi/4 .. Pi/4
151 -- interval, so it is easier for both error analysis and implementation
152 -- to leave the reduction non-strict, and assume the reduced argument is
153 -- within -0.26 * Pi .. 0.26 * Pi rather than a quarter of pi.
154
155 -- The reduction is guaranteed to be correct to within 0.501 ulp for
156 -- values of X for which Ada's accuracy guarantees apply:
157 -- abs X <= 2.0**(T'Machine_Mantissa / 2)
158 -- For values outside this range, an attempt is made to have significance
159 -- decrease only proportionally with increase of magnitued. In any case,
160 -- for all finite arguments, the reduction will succeed, though the reduced
161 -- value may not agree with the mathematically correct value in even its
162 -- sign.
163
164 procedure Reduce_Half_Pi_Large (X : in out F; N : F; Q : out Quadrant);
165
166 procedure Split_Veltkamp (X : Float; X_Hi, X_Lo : out Float)
167 with Post => X = X_Hi + X_Lo;
168
169 function Multiply_Add (X, Y, Z : F) return F is (X * Y + Z);
170
171 ---------------------
172 -- Reconstruct_Pow --
173 ---------------------
174
175 function Reconstruct_Pow
176 (Z : Float;
177 P : Integer;
178 M : Integer) return Float
179 is
180 -- Cody and Waite implementation (in "**" function page 84)
181
182 Result : Float;
183
184 begin
185 Result := A1_Tab_F (P + 1) * Z;
186 return Float'Scaling (Result, M);
187 end Reconstruct_Pow;
188
189 ----------------
190 -- Reduce_044 --
191 ----------------
192
193 procedure Reduce_044
194 (X : Float;
195 Reduced_X : out Float;
196 P : out Integer)
197 is
198 -- Cody and Waite implementation (in "**" function page 84)
199
200 begin
201 P := 1;
202
203 if X <= A1_Tab_F (9) then
204 P := 9;
205 end if;
206
207 if X <= A1_Tab_F (P + 4) then
208 P := P + 4;
209 end if;
210
211 if X <= A1_Tab_F (P + 2) then
212 P := P + 2;
213 end if;
214
215 Reduced_X := (X - A1_Tab_F (P + 1)) - A2_Tab_F ((P + 1) / 2);
216 Reduced_X := Reduced_X / (X + A1_Tab_F (P + 1));
217 Reduced_X := Reduced_X + Reduced_X;
218 end Reduce_044;
219
220 -----------------
221 -- Reduce_Ln_2 --
222 -----------------
223
224 procedure Reduce_Ln_2 (X : in out Float; N : out Integer) is
225 L1 : constant := Float'Leading_Part (Ln_2, F_HM);
226 L2 : constant := Float'Leading_Part (Ln_2 - L1, F_HM);
227 L3 : constant := Ln_2 - L1 - L2;
228 XN : constant Float := Float'Rounding (X * Inv_Ln_2);
229
230 begin
231 -- The argument passed to the function is smaller than Ymax * 1/log(2)
232 -- No overflow is possible for N (Ymax is the largest machine number
233 -- less than Log (F'Last)).
234
235 N := Integer (XN);
236 X := ((X - XN * L1) - XN * L2) - XN * L3;
237
238 if X < -Ln_2 / 2.0 then
239 X := X + Ln_2;
240 N := N - 1;
241 end if;
242
243 if X > Ln_2 / 2.0 then
244 X := X - Ln_2;
245 N := N + 1;
246 end if;
247 end Reduce_Ln_2;
248
249 --------------------
250 -- Split_Veltkamp --
251 --------------------
252
253 procedure Split_Veltkamp (X : Float; X_Hi, X_Lo : out Float) is
254 M : constant F := 0.5 + 2.0**(1 - F'Machine_Mantissa / 2);
255 begin
256 X_Hi := X * M - (X * M - X);
257 X_Lo := X - X_Hi;
258 end Split_Veltkamp;
259
260 --------------------------
261 -- Reduce_Half_Pi_Large --
262 --------------------------
263
264 procedure Reduce_Half_Pi_Large (X : in out F; N : F; Q : out Quadrant) is
265 type Int_64 is range -2**63 .. 2**63 - 1; -- used for conversions
266
267 HM : constant Positive := F'Machine_Mantissa / 2;
268 C1 : constant F := F'Leading_Part (Half_Pi, HM);
269 C2 : constant F := F'Leading_Part (Half_Pi - C1, HM);
270 C3 : constant F := F'Leading_Part (Half_Pi - C1 - C2, HM);
271 C4 : constant F := Half_Pi - C1 - C2 - C3;
272
273 K : F := N;
274 K_Hi : F;
275 K_Lo : F;
276
277 begin
278 Q := 0;
279
280 loop
281 Split_Veltkamp (X => K, X_Hi => K_Hi, X_Lo => K_Lo);
282
283 X := Multiply_Add (-K_Hi, C1, X);
284 X := Multiply_Add (-K_Hi, C2, Multiply_Add (-K_Lo, C1, X));
285 X := Multiply_Add (-K_Hi, C3, Multiply_Add (-K_Lo, C2, X));
286 X := Multiply_Add (-K_Hi, C4, Multiply_Add (-K_Lo, C3, X));
287 X := Multiply_Add (-K_Lo, C4, X);
288
289 if abs K < 2.0**62 then
290 Q := Quadrant ((Int_64 (Q) + Int_64 (N)) mod 4);
291
292 elsif abs K_Lo <= 2.0**62 then
293 Q := Quadrant ((Int_64 (Q) + Int_64 (N)) mod 4);
294 end if;
295
296 exit when X in -0.26 * Pi .. 0.26 * Pi;
297
298 K := F'Machine_Rounding (X * Half_Pi**(-1));
299 end loop;
300 end Reduce_Half_Pi_Large;
301
302 --------------------
303 -- Reduce_Half_Pi --
304 --------------------
305
306 procedure Reduce_Half_Pi (X : in out Float; Q : out Quadrant) is
307 K : constant := Pi / 2.0;
308 Bits_N : constant := 9;
309 Max_N : constant := 2.0**Bits_N - 1.0;
310 Max_X : constant Float := Float'Pred (K * Max_N);
311 Bits_C : constant := Float'Machine_Mantissa - Bits_N;
312 C1 : constant Float := Float'Leading_Part (K, Bits_C);
313 C2 : constant Float := Float'Leading_Part (K - C1, Bits_C);
314 C3 : constant Float := Float'Leading_Part (K - C1 - C2, Bits_C);
315 C4 : constant Float := K - C1 - C2 - C3;
316 N : constant Float := Float'Machine_Rounding (X * K**(-1));
317
318 begin
319 if abs X < Max_X then
320
321 -- X is in the range for which strict accuracy can be provided
322
323 X := (((X - N * C1) - N * C2) - N * C3) - N * C4;
324 Q := Integer (N) mod 4;
325
326 elsif not X'Valid then
327
328 -- X is an infinity or NaN
329
330 X := X - X;
331 Q := 0;
332
333 else
334 Reduce_Half_Pi_Large (X, N, Q);
335 end if;
336 end Reduce_Half_Pi;
337
338 --------------------
339 -- Instantiations --
340 --------------------
341
342 package Instantiations is
343 function Acos is new Generic_Acos (F);
344 function Atan2 is new Generic_Atan2 (F);
345 end Instantiations;
346
347 ----------
348 -- Acos --
349 ----------
350
351 function Acos (X : F) return F is (Instantiations.Acos (X));
352
353 -----------
354 -- Acosh --
355 -----------
356
357 function Acosh (X : F) return F is
358
359 -- Math based implementation using Log1p: x-> Log (1+x)
360
361 T : constant F := X - 1.0;
362
363 begin
364 if X > 1.0 / Sqrt_Epsilon_F then
365 return Log (X) + Ln_2;
366 elsif X < 2.0 then
367 return Log1p (T + Sqrt (2.0 * T + T * T));
368 else
369 return Log (X + Sqrt ((X - 1.0) * (X + 1.0)));
370 end if;
371 end Acosh;
372
373 ----------
374 -- Asin --
375 ----------
376
377 function Asin (X : F) return F is (Float_Approximations.Asin (X));
378
379 -----------
380 -- Asinh --
381 -----------
382
383 function Asinh (X : F) return F is
384
385 -- Math based implementation using Log1p: x-> Log (1+x)
386
387 Y : constant F := abs X;
388 G : constant F := X * X;
389 Res : F;
390
391 begin
392 if Y < Sqrt_Epsilon_F then
393 Res := Y;
394 elsif Y > 1.0 / Sqrt_Epsilon_F then
395 Res := Log (Y) + Ln_2;
396 elsif Y < 2.0 then
397 Res := Log1p (Y + G / (1.0 + Sqrt (G + 1.0)));
398 else
399 Res := Log (Y + Sqrt (G + 1.0));
400 end if;
401
402 return F'Copy_Sign (Res, X);
403 end Asinh;
404
405 ----------
406 -- Atan --
407 ----------
408
409 function Atan (X : F) return F is (Instantiations.Atan2 (X, 1.0));
410
411 -----------
412 -- Atan2 --
413 -----------
414
415 function Atan2 (Y, X : F) return F is (Instantiations.Atan2 (Y, X));
416
417 -----------
418 -- Atanh --
419 -----------
420
421 function Atanh (X : F) return F is
422
423 -- Math based implementation using Log1p: x-> Log (1+x)
424
425 (if X >= 0.0
426 then Log1p (2.0 * X / (1.0 - X)) / 2.0
427 else -Log1p (-2.0 * X / (1.0 + X)) / 2.0);
428
429 ---------
430 -- Cos --
431 ---------
432
433 function Cos (X : F) return F is
434
435 -- Math based implementation using Hart constants
436
437 Y : F := abs (X);
438 Q : Quadrant;
439 Result : F;
440
441 begin
442 Reduce_Half_Pi (Y, Q);
443
444 if Q mod 2 = 0 then
445 Result := Approx_Cos (Y);
446 else
447 Result := Approx_Sin (Y);
448 end if;
449
450 return (if Q = 1 or else Q = 2 then -Result else Result);
451 end Cos;
452
453 ----------
454 -- Cosh --
455 ----------
456
457 function Cosh (X : F) return F is
458
459 -- Cody and Waite implementation (page 217)
460
461 Y : constant F := abs (X);
462
463 -- Because the overflow threshold for cosh(X) is beyond the overflow
464 -- threshold for exp(X), it appears natural to reformulate the
465 -- computation as:
466
467 -- Cosh (X) = Exp (X - Log (2))
468
469 -- But because Log (2) is not an exact machine number, the finite word
470 -- length of the machine implies that the absolute error in X - Log (2),
471 -- and hence the transmitted error in Cosh (X) is proportional to the
472 -- magnitude of X even when X is error-free.
473
474 -- To avoid this problem, we revise the computation to:
475
476 -- Cosh (X) = V/2 * exp(X - Log (V))
477
478 -- where Log (V) is an exact machine number slightly larger than Log (2)
479 -- with the last few digits of its significand zero.
480
481 Ln_V : constant := 8#0.542714#;
482 -- Machine value slightly above Ln_2
483
484 V_2 : constant := 0.24999_30850_04514_99336;
485 -- V**(-2)
486
487 V_2_1 : constant := 0.13830_27787_96019_02638E-4;
488 -- V / 2 - 1
489
490 Y_Bar : constant := 88.72283_93554_69;
491 -- The last floating point for which exp (Y) does not overflow and
492 -- exp (-Y) does not underflow.
493
494 W : F;
495 Z : F;
496
497 begin
498 if Y >= Y_Bar then
499 W := Y - Ln_V;
500 Z := Exp (W);
501 Z := Z + V_2 / Z;
502
503 return Z + V_2_1 * Z; -- rewriting of V/2 * Z
504
505 else
506 Z := Exp (X);
507 return (Z + 1.0 / Z) / 2.0;
508 end if;
509 end Cosh;
510
511 ---------
512 -- Exp --
513 ---------
514
515 function Exp (X : F) return F is
516
517 -- Cody and Waite implementation (page 60)
518
519 N : Integer;
520 Y : F := X;
521 R : F;
522
523 Ymax : constant Float := F'Pred (88.72283_90520_7);
524 -- The largest machine number less than Log (F'Last)
525
526 begin
527
528 if abs (Y) < 2.0**(-F'Machine_Mantissa - 1) then
529 return 1.0;
530 end if;
531
532 if abs Y > Ymax then
533 return (if Y > 0.0 then Infinity else 0.0);
534 end if;
535
536 Reduce_Ln_2 (Y, N);
537 R := Approx_Exp (Y);
538 return Float'Scaling (R, N);
539 end Exp;
540
541 ----------
542 -- Exp2 --
543 ----------
544
545 function Exp2 (X : F) return F is
546
547 -- Implementation based on Cody and Waite Exp implementation (page 217)
548 -- but using Hart constants
549
550 N : Integer;
551 Y : F := X;
552 R : F;
553
554 Result : F;
555
556 begin
557 if abs (Y) < 2.0**(-F'Machine_Mantissa - 1) then
558 return 1.0;
559 end if;
560
561 if abs Y > F'Pred (F (F'Machine_Emax)) then
562 return (if Y > 0.0 then Infinity else 0.0);
563 end if;
564
565 N := Integer (X);
566 Y := Y - Float (N);
567 R := Approx_Exp2 (Y);
568 Result := Float'Scaling (R, N + 1);
569
570 if Result /= Result then
571 return (if X < F'First then 0.0 else Infinity);
572 else
573 return Result;
574 end if;
575 end Exp2;
576
577 ---------
578 -- Log --
579 ---------
580
581 function Log (X : F) return F is
582
583 -- Cody and Waite implementation (page 35)
584
585 Exponent_X : constant Integer := F'Exponent (X);
586 XN : F := F (Exponent_X);
587 Mantissa_X : F := F'Scaling (X, -Exponent_X);
588 HM_F : constant Integer := Integer (F'Machine_Mantissa / 2);
589 L1 : constant F := F'Leading_Part (Ln_2, HM_F);
590 L2 : constant F := Ln_2 - L1;
591 Result : F;
592
593 begin
594 if X <= 0.0 then
595 if X < 0.0 then
596 return NaN;
597 else
598 return -Infinity;
599 end if;
600
601 elsif X > Float'Last then
602 return X;
603
604 -- Making sure X is in Sqrt(0.5) .. Sqrt (2)
605
606 elsif Mantissa_X <= Sqrt_Half then
607 XN := XN - 1.0;
608 Mantissa_X := Mantissa_X * 2.0;
609 end if;
610
611 Result := Approx_Log (Mantissa_X);
612 Result := (XN * L2 + Result) + XN * L1;
613
614 return Result;
615 end Log;
616
617 -----------
618 -- Log1p --
619 -----------
620
621 function Log1p (X : Float) return Float is
622
623 -- Quick implementation of Log1p not accurate to the Ada regular Log
624 -- requirements, but accurate enough to compute inverse hyperbolic
625 -- functions.
626
627 begin
628 if 1.0 + X = 1.0 then
629 return X;
630 elsif X > Float'Last then
631 return X;
632 else
633 return Log (1.0 + X) * (X / ((1.0 + X) - 1.0));
634 end if;
635 end Log1p;
636
637 ----------
638 -- Log2 --
639 ----------
640
641 function Log2 (X : F) return F is
642
643 -- Quick implementation of Log2 not accurate to the Ada regular Log
644 -- (base e) requirement on the whole definition interval but accurate
645 -- enough on 0 .. 2**(-1/16).
646
647 (Log (X) * (1.0 / Ln_2));
648
649 ----------
650 -- Pow --
651 ----------
652
653 function Pow (Left, Right : F) return F is
654
655 -- Cody and Waite implementation (page 84)
656
657 procedure Pow_Special_Cases is new Generic_Pow_Special_Cases (F);
658
659 One_Over_Sixteen : constant := 0.0625;
660
661 M : constant Integer := F'Exponent (Left);
662 G : constant F := F'Fraction (Left);
663 Y : constant F := Right;
664 Z : F;
665 P : Integer;
666
667 U2, U1, Y1, Y2, W1, W2, W : F;
668 MM, PP, IW1, I : Integer;
669 Is_Special : Boolean;
670 Special_Result : F;
671
672 begin
673 -- Special values
674
675 Pow_Special_Cases (Left, Right, Is_Special, Special_Result);
676
677 if Is_Special then
678 return Special_Result;
679 else
680 -- Left**Right is calculated using the formula
681 -- 2**(Right * Log2 (Left))
682
683 Reduce_044 (G, Z, P);
684
685 -- At this point, Z <= 0.044
686
687 U2 := Approx_Power_Log (Z);
688 U1 := F (M * 16 - P) * 0.0625; -- U2 + U1 = Log2 (Left)
689
690 -- Forming the pseudo extended precision product of U * Right
691
692 Y1 := Reduce_1_16 (Y);
693 Y2 := Y - Y1;
694
695 W := U2 * Y + U1 * Y2;
696 W1 := Reduce_1_16 (W);
697 W2 := W - W1;
698
699 W := W1 + U1 * Y1;
700 W1 := Reduce_1_16 (W);
701 W2 := W2 + (W - W1);
702
703 W := Reduce_1_16 (W2);
704 IW1 := Integer (16.0 * (W1 + W));
705 W2 := W2 - W;
706
707 if W2 > 0.0 then
708 W2 := W2 - One_Over_Sixteen;
709 IW1 := 1 + IW1;
710 end if;
711
712 if IW1 < 0 then
713 I := 0;
714 else
715 I := 1;
716 end if;
717
718 MM := Integer (IW1 / 16) + I;
719 PP := 16 * MM - IW1;
720
721 Z := Approx_Exp2 (W2);
722 return Reconstruct_Pow (Z, PP, MM);
723 end if;
724 end Pow;
725
726 ---------
727 -- Sin --
728 ---------
729
730 function Sin (X : F) return F is
731
732 -- Math based implementation using Hart constants
733
734 Y : F := abs X;
735 Q : Quadrant;
736 Result : F;
737
738 begin
739 Reduce_Half_Pi (Y, Q);
740
741 Result := (if Q mod 2 = 0 then Approx_Sin (Y) else Approx_Cos (Y));
742
743 return F'Copy_Sign (1.0, X) * (if Q >= 2 then -Result else Result);
744 end Sin;
745
746 -----------
747 -- Sinh --
748 -----------
749
750 function Sinh (X : F) return F is
751
752 -- Cody and Waite implementation (page 217)
753
754 Sign : constant F := F'Copy_Sign (1.0, X);
755 Y : constant F := abs (X);
756
757 -- Because the overflow threshold for sinh(X) is beyond the overflow
758 -- threshold for exp(X), it appears natural to reformulate the
759 -- computation as:
760
761 -- Sinh (X) = Exp (X - Log (2))
762
763 -- But because Log (2) is not an exact machine number, the finite word
764 -- length of the machine implies that the absolute error in X - Log (2),
765 -- and hence the transmitted error in Sinh (X), is proportional to the
766 -- magnitude of X even when X is error-free.
767
768 -- To avoid this problem, we revise the computation to:
769
770 -- Sinh (X) = V/2 * exp(X - Log (V))
771
772 -- where Log (V) is an exact machine number slightly larger than Log (2)
773 -- with the last few digits of its significand zero.
774
775 Ln_V : constant := 8#0.542714#;
776 -- Machine value slightly above Ln_2
777
778 V_2 : constant := 0.24999_30850_04514_99336;
779 -- V**(-2)
780
781 V_2_1 : constant := 0.13830_27787_96019_02638E-4;
782 -- V / 2 - 1
783
784 Y_Bar : constant := 88.72283_93554_69;
785 -- The last floating point for which exp (X) does not overflow and
786 -- exp (-x) does not underflow
787
788 W : F;
789 Z : F;
790
791 begin
792 if Y < 1.0 then
793 return Approx_Sinh (X);
794 end if;
795
796 if Y >= Y_Bar then
797 W := Y - Ln_V;
798 Z := Exp (W);
799 Z := Z - V_2 / Z;
800
801 return Sign * (Z + V_2_1 * Z); -- rewriting of V/2 * Z
802
803 else
804 Z := Exp (Y);
805 return Sign * ((Z - 1.0 / Z) / 2.0);
806 end if;
807 end Sinh;
808
809 ----------
810 -- Sqrt --
811 ----------
812
813 function Sqrt (X : Float) return Float renames
814 System.Libm_Single.Squareroot.Sqrt;
815
816 ---------
817 -- Tan --
818 ---------
819
820 function Tan (X : F) return F is
821
822 -- Math based implementation using Hart constants
823
824 Y : F := abs X;
825 N : Integer;
826
827 begin
828 Reduce_Half_Pi (Y, N);
829
830 -- The reconstruction is included in the algebraic fraction in
831 -- Approx_Tan function.
832
833 if N mod 2 = 0 then
834 return Approx_Tan (Y) * F'Copy_Sign (1.0, X);
835
836 else
837 return Approx_Cot (Y) * F'Copy_Sign (1.0, X);
838 end if;
839 end Tan;
840
841 ----------
842 -- Tanh --
843 ----------
844
845 function Tanh (X : F) return F is
846
847 -- Cody and Waite implementation (page 239)
848
849 Y : constant F := abs (X);
850 Xbig : constant := Ln_2 * F (1 + F'Machine_Mantissa);
851 LN_3_2 : constant := 0.54930_61443_34054_84570;
852 Result : F;
853
854 begin
855 if Y > Xbig then
856 Result := 1.0;
857 else
858 if Y > LN_3_2 then
859 Result := 1.0 - 2.0 / (Exp (2.0 * Y) + 1.0);
860 else
861 Result := Approx_Tanh (Y);
862 end if;
863 end if;
864
865 return F'Copy_Sign (Result, X);
866 end Tanh;
867
868 end System.Libm_Single;