File : s-libm-ada.adb
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT COMPILER COMPONENTS --
4 -- --
5 -- S Y S T E M . L I B M --
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-libm.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
45 package body System.Libm is
46 type Unsigned_64 is mod 2**64;
47
48 generic
49 type T is private;
50 with function Multiply_Add (X, Y, Z : T) return T is <>;
51 -- The Multiply_Add function returns the value of X * Y + Z, ideally
52 -- (but not necessarily) using a wider intermediate type, or a fused
53 -- multiply-add operation with only a single rounding. They are used
54 -- for evaluating polynomials.
55 package Generic_Polynomials is
56
57 type Polynomial is array (Natural range <>) of T;
58 -- A value P of type PolynomialRepresents the polynomial
59 -- P (X) = P_0 + P_1 * X + ... + P_(n-1) * X**(n-1) + P_n * X**n,
60 --
61 -- where n = P'Length - 1, P_0 is P (P'First) and P_n is P (P'Last)
62
63 -- P (X) = P_0 + X * (P_1 + X * (P_2 + X * (... + X * P_n)))
64
65 function Compute_Horner (P : Polynomial; X : T) return T with Inline;
66 -- Computes the polynomial P using the Horner scheme:
67 -- P (X) = P_0 + X * (P_1 + X * (P_2 + X * (... + X * P_n)))
68 end Generic_Polynomials;
69
70 ------------------------
71 -- Generic_Polynomial --
72 ------------------------
73
74 package body Generic_Polynomials is
75
76 --------------------
77 -- Compute_Horner --
78 ---------------------
79
80 function Compute_Horner (P : Polynomial; X : T) return T is
81 Result : T := P (P'Last);
82
83 begin
84 for P_j of reverse P (P'First .. P'Last - 1) loop
85 Result := Multiply_Add (Result, X, P_j);
86 end loop;
87
88 return Result;
89 end Compute_Horner;
90 end Generic_Polynomials;
91
92 ----------------------------------
93 -- Generic_Float_Approximations --
94 ----------------------------------
95
96 package body Generic_Approximations is
97
98 function Multiply_Add (X, Y, Z : T) return T is (X * Y + Z);
99
100 package Float_Polynomials is new Generic_Polynomials (T);
101 use Float_Polynomials;
102
103 -----------------
104 -- Approx_Asin --
105 -----------------
106
107 function Approx_Asin (X : T) return T is
108 P : T;
109 Q : T;
110
111 begin
112 if Mantissa <= 24 then
113 declare
114
115 -- Approximation MRE = 6.0128E-9
116
117 P1 : constant T := Exact (0.93393_5835);
118 P2 : constant T := Exact (-0.50440_0557);
119
120 Q0 : constant T := Exact (5.6036_3004);
121 Q1 : constant T := Exact (-5.5484_6723);
122
123 begin
124 P := Compute_Horner ((P1, P2), X);
125 Q := Compute_Horner ((Q0, Q1 + X), X);
126 end;
127
128 else
129 declare
130 -- Approximation MRE = 2.0975E-18
131
132 P1 : constant T := Exact (-0.27368_49452_41642_55994E+2);
133 P2 : constant T := Exact (+0.57208_22787_78917_31407E+2);
134 P3 : constant T := Exact (-0.39688_86299_75048_77339E+2);
135 P4 : constant T := Exact (+0.10152_52223_38064_63645E+2);
136 P5 : constant T := Exact (-0.69674_57344_73506_46411);
137
138 Q0 : constant T := Exact (-0.16421_09671_44985_60795E+3);
139 Q1 : constant T := Exact (+0.41714_43024_82604_12556E+3);
140 Q2 : constant T := Exact (-0.38186_30336_17501_49284E+3);
141 Q3 : constant T := Exact (+0.15095_27084_10306_04719E+3);
142 Q4 : constant T := Exact (-0.23823_85915_36702_38830E+2);
143
144 begin
145 P := Compute_Horner ((P1, P2, P3, P4, P5), X);
146 Q := Compute_Horner ((Q0, Q1, Q2, Q3, Q4 + X), X);
147 end;
148 end if;
149
150 return X * P / Q;
151 end Approx_Asin;
152
153 -----------------
154 -- Approx_Atan --
155 -----------------
156
157 function Approx_Atan (X : T) return T is
158 G : constant T := X * X;
159 P, Q : T;
160
161 begin
162 if Mantissa <= 24 then
163 declare
164 -- Approximation MRE = 3.2002E-9
165
166 P0 : constant T := Exact (-0.47083_25141);
167 P1 : constant T := Exact (-0.50909_58253E-1);
168
169 Q0 : constant T := Exact (0.14125_00740E1);
170
171 begin
172 P := Compute_Horner ((P0, P1), G);
173 Q := Q0 + G;
174 end;
175
176 else
177 declare
178 -- Approximation MRE = 1.8154E-18
179
180 P0 : constant T := Exact (-0.13688_76889_41919_26929E2);
181 P1 : constant T := Exact (-0.20505_85519_58616_51981E2);
182 P2 : constant T := Exact (-0.84946_24035_13206_83534E1);
183 P3 : constant T := Exact (-0.83758_29936_81500_59274);
184
185 Q0 : constant T := Exact (0.41066_30668_25757_81263E2);
186 Q1 : constant T := Exact (0.86157_34959_71302_42515E2);
187 Q2 : constant T := Exact (0.59578_43614_25973_44465E2);
188 Q3 : constant T := Exact (0.15024_00116_00285_76121E2);
189
190 begin
191 P := Compute_Horner ((P0, P1, P2, P3), G);
192 Q := Compute_Horner ((Q0, Q1, Q2, Q3 + G), G);
193 end;
194 end if;
195
196 return Multiply_Add (X, (G * P / Q), X);
197 end Approx_Atan;
198
199 function Approx_Cos (X : T) return T is
200 -- Note: The reference tables named below for cosine lists
201 -- constants for cos((pi/4) * x) ~= P (x^2), in order to get
202 -- cos (x), the constants have been adjusted by division of
203 -- appropriate powers of (pi/4) ^ n, for n 0,2,4,6 etc.
204 Cos_P : constant Polynomial :=
205 (if Mantissa <= 24
206 then
207 -- Hart's constants : #COS 3821# (p. 209)
208 -- Approximation MRE = 8.1948E-9 ???
209
210 (0 => Exact (0.99999_99999),
211 1 => Exact (-0.49999_99957),
212 2 => Exact (0.41666_61323E-1),
213 3 => Exact (-0.13888_52915E-2),
214 4 => Exact (0.24372_67909E-4))
215
216 elsif Mantissa <= 53
217 then
218 -- Hart's constants : #COS 3824# (p. 209)
219 -- Approximation MRE = 1.2548E-18
220
221 (0 => Exact (0.99999_99999_99999_99995),
222 1 => Exact (-0.49999_99999_99999_99279),
223 2 => Exact (+0.04166_66666_66666_430254),
224 3 => Exact (-0.13888_88888_88589_60433E-2),
225 4 => Exact (+0.24801_58728_28994_63022E-4),
226 5 => Exact (-0.27557_31286_56960_82219E-6),
227 6 => Exact (+0.20875_55514_56778_82872E-8),
228 7 => Exact (-0.11352_12320_57839_39582E-10))
229 else
230 -- Hart's constants : #COS 3825# (p. 209-210)
231 -- Approximation MRE = ???
232
233 (0 => Exact (+1.0),
234 1 => Exact (-0.49999_99999_99999_99994_57899),
235 2 => Exact (+0.41666_66666_66666_66467_89627E-1),
236 3 => Exact (-0.13888_88888_88888_57508_03579E-2),
237 4 => Exact (+0.24801_58730_15616_31808_80662E-4),
238 5 => Exact (-0.27557_31921_21557_14660_22522E-6),
239 6 => Exact (+0.20876_75377_75228_35357_18906E-8),
240 7 => Exact (-0.11470_23678_56189_18819_10735E-10),
241 8 => Exact (+0.47358_93914_72413_21156_01793E-13)));
242 begin
243 return Compute_Horner (Cos_P, X * X);
244 end Approx_Cos;
245
246 ----------------
247 -- Approx_Exp --
248 ----------------
249
250 function Approx_Exp (X : T) return T is
251 -- Cody and Waite implementation (page 69)
252 Exp_P : constant Polynomial :=
253 (if Mantissa <= 24
254 then -- Approximation MRE = 8.1529E-10
255 (0 => Exact (0.24999_99995_0),
256 1 => Exact (0.41602_88626_0E-2))
257 elsif Mantissa <= 53
258 then -- Approximation MRE = 1.0259E-17
259 (0 => Exact (0.24999_99999_99999_993),
260 1 => Exact (0.69436_00015_11792_852E-2),
261 2 => Exact (0.16520_33002_68279_130E-4))
262 else
263 (0 => Exact (0.25),
264 1 => Exact (0.75753_18015_94227_76666E-2),
265 2 => Exact (0.31555_19276_56846_46356E-4)));
266
267 Exp_Q : constant Polynomial :=
268 (if Mantissa <= 24
269 then
270 (0 => Exact (0.5),
271 1 => Exact (0.49987_17877_8E-1))
272 elsif Mantissa <= 53
273 then
274 (0 => Exact (0.5),
275 1 => Exact (0.55553_86669_69001_188E-1),
276 2 => Exact (0.49586_28849_05441_294E-3))
277 else
278 (0 => Exact (0.5),
279 1 => Exact (0.56817_30269_85512_21787E-1),
280 2 => Exact (0.63121_89437_43985_03557E-3),
281 3 => Exact (0.75104_02839_98700_46114E-6)));
282
283 G : constant T := X * X;
284 P : T;
285 Q : T;
286
287 begin
288 P := Compute_Horner (Exp_P, G);
289 Q := Compute_Horner (Exp_Q, G);
290
291 return Exact (2.0) * Multiply_Add (X, P / (Multiply_Add (-X, P, Q)),
292 Exact (0.5));
293 end Approx_Exp;
294
295 ----------------
296 -- Approx_Log --
297 ----------------
298
299 function Approx_Log (X : T) return T is
300
301 Log_P : constant Polynomial :=
302 (if Mantissa <= 24
303 then -- Approximation MRE = 1.0368E-10
304 (0 => Exact (-0.46490_62303_464),
305 1 => Exact (0.013600_95468_621))
306 else -- Approximation MRE = 4.7849E-19
307 (0 => Exact (-0.64124_94342_37455_81147E+2),
308 1 => Exact (0.16383_94356_30215_34222E+2),
309 2 => Exact (-0.78956_11288_74912_57267)));
310
311 Log_Q : constant Polynomial :=
312 (if Mantissa <= 24
313 then
314 (0 => Exact (-5.5788_73750_242),
315 1 => Exact (1.0))
316 else
317 (0 => Exact (-0.76949_93210_84948_79777E+3),
318 1 => Exact (0.31203_22209_19245_32844E+3),
319 2 => Exact (-0.35667_97773_90346_46171E+2),
320 3 => Exact (1.0)));
321
322 G : T;
323 P : T;
324 Q : T;
325
326 ZNum, ZDen, Z : T;
327
328 begin
329 ZNum := (X + Exact (-0.5)) + Exact (-0.5);
330 ZDen := X * Exact (0.5) + Exact (0.5);
331 Z := ZNum / ZDen;
332 G := Z * Z;
333 P := Compute_Horner (Log_P, G);
334 Q := Compute_Horner (Log_Q, G);
335 return Multiply_Add (Z, G * (P / Q), Z);
336 end Approx_Log;
337
338 ----------------------
339 -- Approx_Power Log --
340 ----------------------
341
342 function Approx_Power_Log (X : T) return T is
343 Power_Log_P : constant Polynomial :=
344 (if Mantissa <= 24
345 then -- Approximation MRE = 7.9529E-4
346 (1 => Exact (0.83357_541E-1))
347 else -- Approximation MRE = 8.7973E-8
348 (1 => Exact (0.83333_33333_33332_11405E-1),
349 2 => Exact (0.12500_00000_05037_99174E-1),
350 3 => Exact (0.22321_42128_59242_58967E-2),
351 4 => Exact (0.43445_77567_21631_19635E-3)));
352
353 K : constant T := Exact (0.44269_50408_88963_40736);
354 G : constant T := X * X;
355 P : T;
356
357 begin
358 P := Compute_Horner (Power_Log_P, G);
359 P := (P * G) * X;
360 P := Multiply_Add (P, K, P);
361
362 return Multiply_Add (X, K, P) + X;
363 end Approx_Power_Log;
364
365 -----------------
366 -- Approx_Exp2 --
367 -----------------
368
369 function Approx_Exp2 (X : T) return T is
370 Exp2_P : constant Polynomial :=
371 (if Mantissa > 24
372 then -- Approximation MRE = 1.7418E-17
373 (1 => Exact (0.69314_71805_59945_29629),
374 2 => Exact (0.24022_65069_59095_37056),
375 3 => Exact (0.55504_10866_40855_95326E-1),
376 4 => Exact (0.96181_29059_51724_16964E-2),
377 5 => Exact (0.13333_54131_35857_84703E-2),
378 6 => Exact (0.15400_29044_09897_64601E-3),
379 7 => Exact (0.14928_85268_05956_08186E-4))
380 else -- Approximation MRE = 3.3642E-9
381 (1 => Exact (0.69314_675),
382 2 => Exact (0.24018_510),
383 3 => Exact (0.54360_383E-1)));
384 begin
385 return Exact (1.0) + Compute_Horner (Exp2_P, X) * X;
386 end Approx_Exp2;
387
388 ----------------
389 -- Approx_Sin --
390 ----------------
391
392 function Approx_Sin (X : T) return T is
393 -- Note: The reference tables named below for sine lists constants
394 -- for sin((pi/4) * x) ~= x * P (x^2), in order to get sin (x),
395 -- the constants have been adjusted by division of appropriate
396 -- powers of (pi/4) ^ n, for n 1,3,5, etc.
397 Sin_P : constant Polynomial :=
398 (if Mantissa <= 24
399 then
400 -- Hart's constants: #SIN 3040# (p. 199)
401
402 (1 => Exact (-0.16666_65022),
403 2 => Exact (0.83320_16396E-2),
404 3 => Exact (-0.19501_81843E-3))
405 else
406 -- Hart's constants: #SIN 3044# (p. 199)
407 -- Approximation MRE = 2.4262E-18 ???
408
409 (1 => Exact (-0.16666_66666_66666_66628),
410 2 => Exact (0.83333_33333_33332_03335E-2),
411 3 => Exact (-0.19841_26984_12531_05860E-3),
412 4 => Exact (0.27557_31921_33901_68712E-5),
413 5 => Exact (-0.25052_10473_82673_30950E-7),
414 6 => Exact (0.16058_34762_32246_06553E-9),
415 7 => Exact (-0.75778_67884_01271_15607E-12)));
416
417 Sqrt_Epsilon_LLF : constant Long_Long_Float :=
418 Sqrt_2 ** (1 - Long_Long_Float'Machine_Mantissa);
419
420 G : constant T := X * X;
421
422 begin
423 if abs X <= Exact (Sqrt_Epsilon_LLF) then
424 return X;
425 end if;
426
427 return Multiply_Add (X, Compute_Horner (Sin_P, G) * G, X);
428 end Approx_Sin;
429
430 -----------------
431 -- Approx_Sinh --
432 -----------------
433
434 function Approx_Sinh (X : T) return T is
435 Sinh_P : constant Polynomial :=
436 (if Mantissa <= 24
437 then -- Approximation MRE = 2.6841E-8
438 (0 => Exact (-0.71379_3159E1),
439 1 => Exact (-0.19033_3300))
440 else -- Approximation MRE = 4.6429E-18
441 (0 => Exact (-0.35181_28343_01771_17881E6),
442 1 => Exact (-0.11563_52119_68517_68270E5),
443 2 => Exact (-0.16375_79820_26307_51372E3),
444 3 => Exact (-0.78966_12741_73570_99479)));
445
446 Sinh_Q : constant Polynomial :=
447 (if Mantissa <= 24
448 then
449 (0 => Exact (-0.42827_7109E2),
450 1 => Exact (1.0))
451 else
452 (0 => Exact (-0.21108_77005_81062_71242E7),
453 1 => Exact (0.36162_72310_94218_36460E5),
454 2 => Exact (-0.27773_52311_96507_01667E3),
455 3 => Exact (1.0)));
456
457 G : constant T := X * X;
458 P : T;
459 Q : T;
460
461 begin
462 P := Compute_Horner (Sinh_P, G);
463 Q := Compute_Horner (Sinh_Q, G);
464
465 return Multiply_Add (X, (G * P / Q), X);
466 end Approx_Sinh;
467
468 ----------------
469 -- Approx_Tan --
470 ----------------
471
472 function Approx_Tan (X : T) return T is
473 Tan_P : constant Polynomial :=
474 (if Mantissa <= 24
475 then -- Approximation MRE = 2.7824E-8
476 (1 => Exact (-0.95801_7723E-1))
477 else -- Approximation MRE = 3.5167E-18
478 (1 => Exact (-0.13338_35000_64219_60681),
479 2 => Exact (0.34248_87823_58905_89960E-2),
480 3 => Exact (-0.17861_70734_22544_26711E-4)));
481
482 Tan_Q : constant Polynomial :=
483 (if Mantissa <= 24
484 then
485 (0 => Exact (1.0),
486 1 => Exact (-0.42913_5777),
487 2 => Exact (0.97168_5835E-2))
488 else
489 (0 => Exact (1.0),
490 1 => Exact (-0.46671_68333_97552_94240),
491 2 => Exact (0.25663_83228_94401_12864E-1),
492 3 => Exact (-0.31181_53190_70100_27307E-3),
493 4 => Exact (0.49819_43399_37865_12270E-6)));
494
495 G : constant T := X * X;
496 P : constant T := Multiply_Add (X, G * Compute_Horner (Tan_P, G), X);
497 Q : constant T := Compute_Horner (Tan_Q, G);
498
499 begin
500 return P / Q;
501 end Approx_Tan;
502
503 ----------------
504 -- Approx_Cot --
505 ----------------
506
507 function Approx_Cot (X : T) return T is
508 Tan_P : constant Polynomial :=
509 (if Mantissa <= 24
510 then -- Approxmiation MRE = 1.5113E-17
511 (1 => Exact (-0.95801_7723E-1))
512 else
513 (1 => Exact (-0.13338_35000_64219_60681),
514 2 => Exact (0.34248_87823_58905_89960E-2),
515 3 => Exact (-0.17861_70734_22544_26711E-4)));
516
517 Tan_Q : constant Polynomial :=
518 (if Mantissa <= 24
519 then
520 (0 => Exact (1.0),
521 1 => Exact (-0.42913_5777),
522 2 => Exact (0.97168_5835E-2))
523 else
524 (0 => Exact (1.0),
525 1 => Exact (-0.46671_68333_97552_94240),
526 2 => Exact (0.25663_83228_94401_12864E-1),
527 3 => Exact (-0.31181_53190_70100_27307E-3),
528 4 => Exact (0.49819_43399_37865_12270E-6)));
529 G : constant T := X * X;
530 P : constant T := Multiply_Add (X, G * Compute_Horner (Tan_P, G), X);
531 Q : constant T := Compute_Horner (Tan_Q, G);
532
533 begin
534 return -Q / P;
535 end Approx_Cot;
536
537 -----------------
538 -- Approx_Tanh --
539 -----------------
540
541 function Approx_Tanh (X : T) return T is
542 Tanh_P : constant Polynomial :=
543 (if Mantissa <= 24
544 then -- Approximation MRE = 2.7166E-9
545 (0 => Exact (-0.82377_28127),
546 1 => Exact (-0.38310_10665E-2))
547 else -- Approximation MRE = 3.2436E-18
548 (0 => Exact (-0.16134_11902_39962_28053E4),
549 1 => Exact (-0.99225_92967_22360_83313E2),
550 2 => Exact (-0.96437_49277_72254_69787)));
551
552 Tanh_Q : constant Polynomial :=
553 (if Mantissa <= 24
554 then
555 (0 => Exact (2.4713_19654),
556 1 => Exact (1.0))
557 else
558 (0 => Exact (0.48402_35707_19886_88686E4),
559 1 => Exact (0.22337_72071_89623_12926E4),
560 2 => Exact (0.11274_47438_05349_49335E3),
561 3 => Exact (1.0)));
562
563 G : constant T := X * X;
564 P, Q : T;
565
566 begin
567 P := Compute_Horner (Tanh_P, G);
568 Q := Compute_Horner (Tanh_Q, G);
569
570 return Multiply_Add (X, G * P / Q, X);
571 end Approx_Tanh;
572
573 ----------
574 -- Asin --
575 ----------
576
577 function Asin (X : T) return T is
578
579 -- Cody and Waite implementation (page 174)
580
581 Y : T := abs X;
582 G : T;
583 Result : T;
584
585 begin
586 if Y <= Exact (0.5) then
587 Result := X + X * Approx_Asin (X * X);
588
589 else
590 G := (Exact (1.0) + (-Y)) * Exact (0.5);
591 Y := Sqrt (G);
592
593 Result :=
594 Exact (Pi / 2.0) - Exact (2.0) * (Y + Y * Approx_Asin (G));
595
596 if not (Exact (0.0) <= X) then
597 Result := -Result;
598 end if;
599 end if;
600
601 return Result;
602 end Asin;
603
604 end Generic_Approximations;
605
606 ------------------
607 -- Generic_Acos --
608 ------------------
609
610 function Generic_Acos (X : T) return T is
611
612 -- Cody and Waite implementation (page 174)
613
614 Y : T := abs (X);
615 G : T;
616 Result : T;
617
618 begin
619 if Y <= 0.5 then
620
621 -- No reduction needed
622
623 G := Y * Y;
624 Result := T'Copy_Sign (Y + Y * Approx_Asin (G), X);
625 return 0.5 * Pi - Result;
626 end if;
627
628 -- In the reduction step that follows, it is not Y, but rather G that
629 -- is reduced. The reduced G is in 0.0 .. 0.25.
630
631 G := (1.0 - Y) / 2.0;
632 Y := -2.0 * Sqrt (G);
633
634 Result := Y + Y * Approx_Asin (G);
635
636 return (if X < 0.0 then Pi + Result else -Result);
637 end Generic_Acos;
638
639 -------------------
640 -- Generic_Atan2 --
641 -------------------
642
643 function Generic_Atan2 (Y, X : T) return T is
644
645 -- Cody and Waite implementation (page 194)
646
647 F : T;
648 N : Integer := -1;
649
650 -- Default value for N is -1 so that if X=0 or over/underflow
651 -- tests on N are all false.
652
653 Result : T;
654
655 begin
656 if Y = 0.0 then
657 if T'Copy_Sign (1.0, X) < 0.0 then
658 return T'Copy_Sign (Pi, Y);
659 else
660 return T'Copy_Sign (0.0, Y);
661 end if;
662
663 elsif X = 0.0 then
664 return T'Copy_Sign (Half_Pi, Y);
665
666 elsif abs (Y) > T'Last * abs (X) then -- overflow
667 Result := T (Half_Pi);
668
669 elsif abs (X) > T'Last * abs (Y) then -- underflow
670 Result := 0.0;
671
672 elsif abs (X) > T'Last and then abs (Y) > T'Last then
673
674 -- NaN
675
676 if X < 0.0 then
677 return T'Copy_Sign (3.0 * Pi / 4.0, Y);
678 else
679 return T'Copy_Sign (Pi / 4.0, Y);
680 end if;
681
682 else
683 F := abs (Y / X);
684
685 if F > 1.0 then
686 F := 1.0 / F;
687 N := 2;
688 else
689 N := 0;
690 end if;
691
692 if F > 2.0 - Sqrt_3 then
693 F := (((Sqrt_3 - 1.0) * F - 1.0) + F) / (Sqrt_3 + F);
694 N := N + 1;
695 end if;
696
697 Result := Approx_Atan (F);
698 end if;
699
700 if N > 1 then
701 Result := -Result;
702 end if;
703
704 case N is
705 when 1 => Result := Result + Sixth_Pi;
706 when 2 => Result := Result + Half_Pi;
707 when 3 => Result := Result + Third_Pi;
708 when others => null;
709 end case;
710
711 if T'Copy_Sign (1.0, X) < 0.0 then
712 Result := Pi - Result;
713 end if;
714
715 return T'Copy_Sign (Result, Y);
716 end Generic_Atan2;
717
718 procedure Generic_Pow_Special_Cases
719 (Left : T;
720 Right : T;
721 Is_Special : out Boolean;
722 Result : out T)
723 is
724 ------------
725 -- Is_Even --
726 ------------
727
728 function Is_Even (X : T) return Boolean is
729 (abs X >= 2.0**T'Machine_Mantissa
730 or else Unsigned_64 (abs X) mod 2 = 0);
731 pragma Assert (T'Machine_Mantissa <= 64);
732 -- If X is large enough, then X is a multiple of 2. Otherwise,
733 -- conversion to Unsigned_64 is safe, assuming a mantissa of at
734 -- most 64 bits.
735
736 begin
737 Is_Special := True;
738 Result := 0.0;
739
740 -- value 'Result' is not used if the input is
741 -- not a couple of special values
742
743 if Right = 0.0 or else not (Left /= 1.0) then
744 Result := (if Right = 0.0 then 1.0 else Left);
745
746 elsif Left = 0.0 then
747 if Right < 0.0 then
748 if Right = T'Rounding (Right) and then not Is_Even (Right) then
749 Result := 1.0 / Left; -- Infinity with sign of Left
750 else
751 Result := 1.0 / abs Left; -- +Infinity
752 end if;
753
754 else
755 if Right = T'Rounding (Right)
756 and then not Is_Even (Right)
757 then
758 Result := Left;
759 else
760 Result := +0.0;
761 end if;
762 end if;
763
764 elsif abs (Right) > T'Last and then Left = -1.0 then
765 Result := 1.0;
766
767 elsif Left < 0.0
768 and then Left >= T'First
769 and then abs (Right) <= T'Last
770 and then Right /= T'Rounding (Right)
771 then
772 Result := 0.0 / (Left - Left); -- NaN
773
774 elsif Right < T'First then
775 if abs (Left) < 1.0 then
776 Result := -Right; -- Infinity
777 else
778 Result := 0.0; -- Cases where Left=+-1 are dealt with above
779 end if;
780
781 elsif Right > T'Last then
782 if abs (Left) < 1.0 then
783 Result := 0.0;
784 else
785 Result := Right;
786 end if;
787
788 elsif Left > T'Last then
789 if Right < 0.0 then
790 Result := 0.0;
791 else
792 Result := Left;
793 end if;
794
795 elsif Left < T'First then
796 if Right > 0.0 then
797 if Right = T'Rounding (Right)
798 and then not Is_Even (Right)
799 then
800 Result := Left;
801 else
802 Result := -Left; -- -Left = +INF
803 end if;
804 else
805 if Right = T'Rounding (Right)
806 and then not Is_Even (Right)
807 then
808 Result := -0.0;
809 else
810 Result := +0.0;
811 end if;
812 end if;
813
814 else
815 Is_Special := False;
816 end if;
817 end Generic_Pow_Special_Cases;
818
819 end System.Libm;