Daily bump.
[official-gcc.git] / gcc / ada / a-numaux-x86.adb
blob811485d859b866fb651dfc208ae2fbdedc5e0bb4
1 ------------------------------------------------------------------------------
2 -- --
3 -- GNAT RUN-TIME COMPONENTS --
4 -- --
5 -- A D A . N U M E R I C S . A U X --
6 -- --
7 -- B o d y --
8 -- (Machine Version for x86) --
9 -- --
10 -- Copyright (C) 1998-2009, Free Software Foundation, Inc. --
11 -- --
12 -- GNAT is free software; you can redistribute it and/or modify it under --
13 -- terms of the GNU General Public License as published by the Free Soft- --
14 -- ware Foundation; either version 3, or (at your option) any later ver- --
15 -- sion. GNAT is distributed in the hope that it will be useful, but WITH- --
16 -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY --
17 -- or FITNESS FOR A PARTICULAR PURPOSE. --
18 -- --
19 -- As a special exception under Section 7 of GPL version 3, you are granted --
20 -- additional permissions described in the GCC Runtime Library Exception, --
21 -- version 3.1, as published by the Free Software Foundation. --
22 -- --
23 -- You should have received a copy of the GNU General Public License and --
24 -- a copy of the GCC Runtime Library Exception along with this program; --
25 -- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see --
26 -- <http://www.gnu.org/licenses/>. --
27 -- --
28 -- GNAT was originally developed by the GNAT team at New York University. --
29 -- Extensive contributions were provided by Ada Core Technologies Inc. --
30 -- --
31 ------------------------------------------------------------------------------
33 -- File a-numaux.adb <- 86numaux.adb
35 -- This version of Numerics.Aux is for the IEEE Double Extended floating
36 -- point format on x86.
38 with System.Machine_Code; use System.Machine_Code;
40 package body Ada.Numerics.Aux is
42 NL : constant String := ASCII.LF & ASCII.HT;
44 -----------------------
45 -- Local subprograms --
46 -----------------------
48 function Is_Nan (X : Double) return Boolean;
49 -- Return True iff X is a IEEE NaN value
51 function Logarithmic_Pow (X, Y : Double) return Double;
52 -- Implementation of X**Y using Exp and Log functions (binary base)
53 -- to calculate the exponentiation. This is used by Pow for values
54 -- for values of Y in the open interval (-0.25, 0.25)
56 procedure Reduce (X : in out Double; Q : out Natural);
57 -- Implements reduction of X by Pi/2. Q is the quadrant of the final
58 -- result in the range 0 .. 3. The absolute value of X is at most Pi.
60 pragma Inline (Is_Nan);
61 pragma Inline (Reduce);
63 --------------------------------
64 -- Basic Elementary Functions --
65 --------------------------------
67 -- This section implements a few elementary functions that are used to
68 -- build the more complex ones. This ordering enables better inlining.
70 ----------
71 -- Atan --
72 ----------
74 function Atan (X : Double) return Double is
75 Result : Double;
77 begin
78 Asm (Template =>
79 "fld1" & NL
80 & "fpatan",
81 Outputs => Double'Asm_Output ("=t", Result),
82 Inputs => Double'Asm_Input ("0", X));
84 -- The result value is NaN iff input was invalid
86 if not (Result = Result) then
87 raise Argument_Error;
88 end if;
90 return Result;
91 end Atan;
93 ---------
94 -- Exp --
95 ---------
97 function Exp (X : Double) return Double is
98 Result : Double;
99 begin
100 Asm (Template =>
101 "fldl2e " & NL
102 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E)
103 & "fld %%st(0) " & NL
104 & "frndint " & NL -- Integer (X * Log2 (E))
105 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
106 & "fxch " & NL
107 & "f2xm1 " & NL -- 2**(...) - 1
108 & "fld1 " & NL
109 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
110 & "fscale " & NL -- E ** X
111 & "fstp %%st(1) ",
112 Outputs => Double'Asm_Output ("=t", Result),
113 Inputs => Double'Asm_Input ("0", X));
114 return Result;
115 end Exp;
117 ------------
118 -- Is_Nan --
119 ------------
121 function Is_Nan (X : Double) return Boolean is
122 begin
123 -- The IEEE NaN values are the only ones that do not equal themselves
125 return not (X = X);
126 end Is_Nan;
128 ---------
129 -- Log --
130 ---------
132 function Log (X : Double) return Double is
133 Result : Double;
135 begin
136 Asm (Template =>
137 "fldln2 " & NL
138 & "fxch " & NL
139 & "fyl2x " & NL,
140 Outputs => Double'Asm_Output ("=t", Result),
141 Inputs => Double'Asm_Input ("0", X));
142 return Result;
143 end Log;
145 ------------
146 -- Reduce --
147 ------------
149 procedure Reduce (X : in out Double; Q : out Natural) is
150 Half_Pi : constant := Pi / 2.0;
151 Two_Over_Pi : constant := 2.0 / Pi;
153 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
154 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
155 P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
156 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
157 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
158 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
159 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
160 - P4, HM);
161 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
162 K : Double := X * Two_Over_Pi;
163 begin
164 -- For X < 2.0**32, all products below are computed exactly.
165 -- Due to cancellation effects all subtractions are exact as well.
166 -- As no double extended floating-point number has more than 75
167 -- zeros after the binary point, the result will be the correctly
168 -- rounded result of X - K * (Pi / 2.0).
170 while abs K >= 2.0**HM loop
171 K := K * M - (K * M - K);
172 X := (((((X - K * P1) - K * P2) - K * P3)
173 - K * P4) - K * P5) - K * P6;
174 K := X * Two_Over_Pi;
175 end loop;
177 if K /= K then
179 -- K is not a number, because X was not finite
181 raise Constraint_Error;
182 end if;
184 K := Double'Rounding (K);
185 Q := Integer (K) mod 4;
186 X := (((((X - K * P1) - K * P2) - K * P3)
187 - K * P4) - K * P5) - K * P6;
188 end Reduce;
190 ----------
191 -- Sqrt --
192 ----------
194 function Sqrt (X : Double) return Double is
195 Result : Double;
197 begin
198 if X < 0.0 then
199 raise Argument_Error;
200 end if;
202 Asm (Template => "fsqrt",
203 Outputs => Double'Asm_Output ("=t", Result),
204 Inputs => Double'Asm_Input ("0", X));
206 return Result;
207 end Sqrt;
209 --------------------------------
210 -- Other Elementary Functions --
211 --------------------------------
213 -- These are built using the previously implemented basic functions
215 ----------
216 -- Acos --
217 ----------
219 function Acos (X : Double) return Double is
220 Result : Double;
222 begin
223 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
225 -- The result value is NaN iff input was invalid
227 if Is_Nan (Result) then
228 raise Argument_Error;
229 end if;
231 return Result;
232 end Acos;
234 ----------
235 -- Asin --
236 ----------
238 function Asin (X : Double) return Double is
239 Result : Double;
241 begin
242 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
244 -- The result value is NaN iff input was invalid
246 if Is_Nan (Result) then
247 raise Argument_Error;
248 end if;
250 return Result;
251 end Asin;
253 ---------
254 -- Cos --
255 ---------
257 function Cos (X : Double) return Double is
258 Reduced_X : Double := abs X;
259 Result : Double;
260 Quadrant : Natural range 0 .. 3;
262 begin
263 if Reduced_X > Pi / 4.0 then
264 Reduce (Reduced_X, Quadrant);
266 case Quadrant is
267 when 0 =>
268 Asm (Template => "fcos",
269 Outputs => Double'Asm_Output ("=t", Result),
270 Inputs => Double'Asm_Input ("0", Reduced_X));
271 when 1 =>
272 Asm (Template => "fsin",
273 Outputs => Double'Asm_Output ("=t", Result),
274 Inputs => Double'Asm_Input ("0", -Reduced_X));
275 when 2 =>
276 Asm (Template => "fcos ; fchs",
277 Outputs => Double'Asm_Output ("=t", Result),
278 Inputs => Double'Asm_Input ("0", Reduced_X));
279 when 3 =>
280 Asm (Template => "fsin",
281 Outputs => Double'Asm_Output ("=t", Result),
282 Inputs => Double'Asm_Input ("0", Reduced_X));
283 end case;
285 else
286 Asm (Template => "fcos",
287 Outputs => Double'Asm_Output ("=t", Result),
288 Inputs => Double'Asm_Input ("0", Reduced_X));
289 end if;
291 return Result;
292 end Cos;
294 ---------------------
295 -- Logarithmic_Pow --
296 ---------------------
298 function Logarithmic_Pow (X, Y : Double) return Double is
299 Result : Double;
300 begin
301 Asm (Template => "" -- X : Y
302 & "fyl2x " & NL -- Y * Log2 (X)
303 & "fld %%st(0) " & NL -- Y * Log2 (X) : Y * Log2 (X)
304 & "frndint " & NL -- Int (...) : Y * Log2 (X)
305 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...)
306 & "fxch " & NL -- Fract (...) : Int (...)
307 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...)
308 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...)
309 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...)
310 & "fscale ", -- 2**(Fract (...) + Int (...))
311 Outputs => Double'Asm_Output ("=t", Result),
312 Inputs =>
313 (Double'Asm_Input ("0", X),
314 Double'Asm_Input ("u", Y)));
315 return Result;
316 end Logarithmic_Pow;
318 ---------
319 -- Pow --
320 ---------
322 function Pow (X, Y : Double) return Double is
323 type Mantissa_Type is mod 2**Double'Machine_Mantissa;
324 -- Modular type that can hold all bits of the mantissa of Double
326 -- For negative exponents, do divide at the end of the processing
328 Negative_Y : constant Boolean := Y < 0.0;
329 Abs_Y : constant Double := abs Y;
331 -- During this function the following invariant is kept:
332 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
334 Base : Double := X;
336 Exp_High : Double := Double'Floor (Abs_Y);
337 Exp_Mid : Double;
338 Exp_Low : Double;
339 Exp_Int : Mantissa_Type;
341 Factor : Double := 1.0;
343 begin
344 -- Select algorithm for calculating Pow (integer cases fall through)
346 if Exp_High >= 2.0**Double'Machine_Mantissa then
348 -- In case of Y that is IEEE infinity, just raise constraint error
350 if Exp_High > Double'Safe_Last then
351 raise Constraint_Error;
352 end if;
354 -- Large values of Y are even integers and will stay integer
355 -- after division by two.
357 loop
358 -- Exp_Mid and Exp_Low are zero, so
359 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
361 Exp_High := Exp_High / 2.0;
362 Base := Base * Base;
363 exit when Exp_High < 2.0**Double'Machine_Mantissa;
364 end loop;
366 elsif Exp_High /= Abs_Y then
367 Exp_Low := Abs_Y - Exp_High;
368 Factor := 1.0;
370 if Exp_Low /= 0.0 then
372 -- Exp_Low now is in interval (0.0, 1.0)
373 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
375 Exp_Mid := 0.0;
376 Exp_Low := Exp_Low - Exp_Mid;
378 if Exp_Low >= 0.5 then
379 Factor := Sqrt (X);
380 Exp_Low := Exp_Low - 0.5; -- exact
382 if Exp_Low >= 0.25 then
383 Factor := Factor * Sqrt (Factor);
384 Exp_Low := Exp_Low - 0.25; -- exact
385 end if;
387 elsif Exp_Low >= 0.25 then
388 Factor := Sqrt (Sqrt (X));
389 Exp_Low := Exp_Low - 0.25; -- exact
390 end if;
392 -- Exp_Low now is in interval (0.0, 0.25)
394 -- This means it is safe to call Logarithmic_Pow
395 -- for the remaining part.
397 Factor := Factor * Logarithmic_Pow (X, Exp_Low);
398 end if;
400 elsif X = 0.0 then
401 return 0.0;
402 end if;
404 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
406 Exp_Int := Mantissa_Type (Exp_High);
408 -- Standard way for processing integer powers > 0
410 while Exp_Int > 1 loop
411 if (Exp_Int and 1) = 1 then
413 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
415 Factor := Factor * Base;
416 end if;
418 -- Exp_Int is even and Exp_Int > 0, so
419 -- Base**Y = (Base**2)**(Exp_Int / 2)
421 Base := Base * Base;
422 Exp_Int := Exp_Int / 2;
423 end loop;
425 -- Exp_Int = 1 or Exp_Int = 0
427 if Exp_Int = 1 then
428 Factor := Base * Factor;
429 end if;
431 if Negative_Y then
432 Factor := 1.0 / Factor;
433 end if;
435 return Factor;
436 end Pow;
438 ---------
439 -- Sin --
440 ---------
442 function Sin (X : Double) return Double is
443 Reduced_X : Double := X;
444 Result : Double;
445 Quadrant : Natural range 0 .. 3;
447 begin
448 if abs X > Pi / 4.0 then
449 Reduce (Reduced_X, Quadrant);
451 case Quadrant is
452 when 0 =>
453 Asm (Template => "fsin",
454 Outputs => Double'Asm_Output ("=t", Result),
455 Inputs => Double'Asm_Input ("0", Reduced_X));
456 when 1 =>
457 Asm (Template => "fcos",
458 Outputs => Double'Asm_Output ("=t", Result),
459 Inputs => Double'Asm_Input ("0", Reduced_X));
460 when 2 =>
461 Asm (Template => "fsin",
462 Outputs => Double'Asm_Output ("=t", Result),
463 Inputs => Double'Asm_Input ("0", -Reduced_X));
464 when 3 =>
465 Asm (Template => "fcos ; fchs",
466 Outputs => Double'Asm_Output ("=t", Result),
467 Inputs => Double'Asm_Input ("0", Reduced_X));
468 end case;
470 else
471 Asm (Template => "fsin",
472 Outputs => Double'Asm_Output ("=t", Result),
473 Inputs => Double'Asm_Input ("0", Reduced_X));
474 end if;
476 return Result;
477 end Sin;
479 ---------
480 -- Tan --
481 ---------
483 function Tan (X : Double) return Double is
484 Reduced_X : Double := X;
485 Result : Double;
486 Quadrant : Natural range 0 .. 3;
488 begin
489 if abs X > Pi / 4.0 then
490 Reduce (Reduced_X, Quadrant);
492 if Quadrant mod 2 = 0 then
493 Asm (Template => "fptan" & NL
494 & "ffree %%st(0)" & NL
495 & "fincstp",
496 Outputs => Double'Asm_Output ("=t", Result),
497 Inputs => Double'Asm_Input ("0", Reduced_X));
498 else
499 Asm (Template => "fsincos" & NL
500 & "fdivp %%st, %%st(1)" & NL
501 & "fchs",
502 Outputs => Double'Asm_Output ("=t", Result),
503 Inputs => Double'Asm_Input ("0", Reduced_X));
504 end if;
506 else
507 Asm (Template =>
508 "fptan " & NL
509 & "ffree %%st(0) " & NL
510 & "fincstp ",
511 Outputs => Double'Asm_Output ("=t", Result),
512 Inputs => Double'Asm_Input ("0", Reduced_X));
513 end if;
515 return Result;
516 end Tan;
518 ----------
519 -- Sinh --
520 ----------
522 function Sinh (X : Double) return Double is
523 begin
524 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
526 if abs X < 25.0 then
527 return (Exp (X) - Exp (-X)) / 2.0;
528 else
529 return Exp (X) / 2.0;
530 end if;
531 end Sinh;
533 ----------
534 -- Cosh --
535 ----------
537 function Cosh (X : Double) return Double is
538 begin
539 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
541 if abs X < 22.0 then
542 return (Exp (X) + Exp (-X)) / 2.0;
543 else
544 return Exp (X) / 2.0;
545 end if;
546 end Cosh;
548 ----------
549 -- Tanh --
550 ----------
552 function Tanh (X : Double) return Double is
553 begin
554 -- Return the Hyperbolic Tangent of x
556 -- x -x
557 -- e - e Sinh (X)
558 -- Tanh (X) is defined to be ----------- = --------
559 -- x -x Cosh (X)
560 -- e + e
562 if abs X > 23.0 then
563 return Double'Copy_Sign (1.0, X);
564 end if;
566 return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
567 end Tanh;
569 end Ada.Numerics.Aux;