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