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