1 ------------------------------------------------------------------------------
3 -- GNAT RUN-TIME COMPONENTS --
5 -- A D A . N U M E R I C S . A U X --
8 -- (Machine Version for x86) --
10 -- Copyright (C) 1998-2006, Free Software Foundation, Inc. --
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 2, 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. See the GNU General Public License --
18 -- for more details. You should have received a copy of the GNU General --
19 -- Public License distributed with GNAT; see file COPYING. If not, write --
20 -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, --
21 -- Boston, MA 02110-1301, USA. --
23 -- As a special exception, if other files instantiate generics from this --
24 -- unit, or you link this unit with other files to produce an executable, --
25 -- this unit does not by itself cause the resulting executable to be --
26 -- covered by the GNU General Public License. This exception does not --
27 -- however invalidate any other reasons why the executable file might be --
28 -- covered by the GNU Public License. --
30 -- GNAT was originally developed by the GNAT team at New York University. --
31 -- Extensive contributions were provided by Ada Core Technologies Inc. --
33 ------------------------------------------------------------------------------
35 -- File a-numaux.adb <- 86numaux.adb
37 -- This version of Numerics.Aux is for the IEEE Double Extended floating
38 -- point format on x86.
40 with System.Machine_Code; use System.Machine_Code;
42 package body Ada.Numerics.Aux is
44 NL : constant String := ASCII.LF & ASCII.HT;
46 -----------------------
47 -- Local subprograms --
48 -----------------------
50 function Is_Nan (X : Double) return Boolean;
51 -- Return True iff X is a IEEE NaN value
53 function Logarithmic_Pow (X, Y : Double) return Double;
54 -- Implementation of X**Y using Exp and Log functions (binary base)
55 -- to calculate the exponentiation. This is used by Pow for values
56 -- for values of Y in the open interval (-0.25, 0.25)
58 procedure Reduce (X : in out Double; Q : out Natural);
59 -- Implements reduction of X by Pi/2. Q is the quadrant of the final
60 -- result in the range 0 .. 3. The absolute value of X is at most Pi.
62 pragma Inline (Is_Nan);
63 pragma Inline (Reduce);
65 --------------------------------
66 -- Basic Elementary Functions --
67 --------------------------------
69 -- This section implements a few elementary functions that are used to
70 -- build the more complex ones. This ordering enables better inlining.
76 function Atan (X : Double) return Double is
83 Outputs => Double'Asm_Output ("=t", Result),
84 Inputs => Double'Asm_Input ("0", X));
86 -- The result value is NaN iff input was invalid
88 if not (Result = Result) then
99 function Exp (X : Double) return Double is
104 & "fmulp %%st, %%st(1)" & NL -- X * log2 (E)
105 & "fld %%st(0) " & NL
106 & "frndint " & NL -- Integer (X * Log2 (E))
107 & "fsubr %%st, %%st(1)" & NL -- Fraction (X * Log2 (E))
109 & "f2xm1 " & NL -- 2**(...) - 1
111 & "faddp %%st, %%st(1)" & NL -- 2**(Fraction (X * Log2 (E)))
112 & "fscale " & NL -- E ** X
114 Outputs => Double'Asm_Output ("=t", Result),
115 Inputs => Double'Asm_Input ("0", X));
123 function Is_Nan (X : Double) return Boolean is
125 -- The IEEE NaN values are the only ones that do not equal themselves
134 function Log (X : Double) return Double is
142 Outputs => Double'Asm_Output ("=t", Result),
143 Inputs => Double'Asm_Input ("0", X));
151 procedure Reduce (X : in out Double; Q : out Natural) is
152 Half_Pi : constant := Pi / 2.0;
153 Two_Over_Pi : constant := 2.0 / Pi;
155 HM : constant := Integer'Min (Double'Machine_Mantissa / 2, Natural'Size);
156 M : constant Double := 0.5 + 2.0**(1 - HM); -- Splitting constant
157 P1 : constant Double := Double'Leading_Part (Half_Pi, HM);
158 P2 : constant Double := Double'Leading_Part (Half_Pi - P1, HM);
159 P3 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2, HM);
160 P4 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3, HM);
161 P5 : constant Double := Double'Leading_Part (Half_Pi - P1 - P2 - P3
163 P6 : constant Double := Double'Model (Half_Pi - P1 - P2 - P3 - P4 - P5);
164 K : Double := X * Two_Over_Pi;
166 -- For X < 2.0**32, all products below are computed exactly.
167 -- Due to cancellation effects all subtractions are exact as well.
168 -- As no double extended floating-point number has more than 75
169 -- zeros after the binary point, the result will be the correctly
170 -- rounded result of X - K * (Pi / 2.0).
172 while abs K >= 2.0**HM loop
173 K := K * M - (K * M - K);
174 X := (((((X - K * P1) - K * P2) - K * P3)
175 - K * P4) - K * P5) - K * P6;
176 K := X * Two_Over_Pi;
181 -- K is not a number, because X was not finite
183 raise Constraint_Error;
186 K := Double'Rounding (K);
187 Q := Integer (K) mod 4;
188 X := (((((X - K * P1) - K * P2) - K * P3)
189 - K * P4) - K * P5) - K * P6;
196 function Sqrt (X : Double) return Double is
201 raise Argument_Error;
204 Asm (Template => "fsqrt",
205 Outputs => Double'Asm_Output ("=t", Result),
206 Inputs => Double'Asm_Input ("0", X));
211 --------------------------------
212 -- Other Elementary Functions --
213 --------------------------------
215 -- These are built using the previously implemented basic functions
221 function Acos (X : Double) return Double is
225 Result := 2.0 * Atan (Sqrt ((1.0 - X) / (1.0 + X)));
227 -- The result value is NaN iff input was invalid
229 if Is_Nan (Result) then
230 raise Argument_Error;
240 function Asin (X : Double) return Double is
244 Result := Atan (X / Sqrt ((1.0 - X) * (1.0 + X)));
246 -- The result value is NaN iff input was invalid
248 if Is_Nan (Result) then
249 raise Argument_Error;
259 function Cos (X : Double) return Double is
260 Reduced_X : Double := abs X;
262 Quadrant : Natural range 0 .. 3;
265 if Reduced_X > Pi / 4.0 then
266 Reduce (Reduced_X, Quadrant);
270 Asm (Template => "fcos",
271 Outputs => Double'Asm_Output ("=t", Result),
272 Inputs => Double'Asm_Input ("0", Reduced_X));
274 Asm (Template => "fsin",
275 Outputs => Double'Asm_Output ("=t", Result),
276 Inputs => Double'Asm_Input ("0", -Reduced_X));
278 Asm (Template => "fcos ; fchs",
279 Outputs => Double'Asm_Output ("=t", Result),
280 Inputs => Double'Asm_Input ("0", Reduced_X));
282 Asm (Template => "fsin",
283 Outputs => Double'Asm_Output ("=t", Result),
284 Inputs => Double'Asm_Input ("0", Reduced_X));
288 Asm (Template => "fcos",
289 Outputs => Double'Asm_Output ("=t", Result),
290 Inputs => Double'Asm_Input ("0", Reduced_X));
296 ---------------------
297 -- Logarithmic_Pow --
298 ---------------------
300 function Logarithmic_Pow (X, Y : Double) return Double is
303 Asm (Template => "" -- X : Y
304 & "fyl2x " & NL -- Y * Log2 (X)
305 & "fst %%st(1) " & NL -- Y * Log2 (X) : Y * Log2 (X)
306 & "frndint " & NL -- Int (...) : Y * Log2 (X)
307 & "fsubr %%st, %%st(1)" & NL -- Int (...) : Fract (...)
308 & "fxch " & NL -- Fract (...) : Int (...)
309 & "f2xm1 " & NL -- 2**Fract (...) - 1 : Int (...)
310 & "fld1 " & NL -- 1 : 2**Fract (...) - 1 : Int (...)
311 & "faddp %%st, %%st(1)" & NL -- 2**Fract (...) : Int (...)
312 & "fscale " & NL -- 2**(Fract (...) + Int (...))
314 Outputs => Double'Asm_Output ("=t", Result),
316 (Double'Asm_Input ("0", X),
317 Double'Asm_Input ("u", Y)));
325 function Pow (X, Y : Double) return Double is
326 type Mantissa_Type is mod 2**Double'Machine_Mantissa;
327 -- Modular type that can hold all bits of the mantissa of Double
329 -- For negative exponents, do divide at the end of the processing
331 Negative_Y : constant Boolean := Y < 0.0;
332 Abs_Y : constant Double := abs Y;
334 -- During this function the following invariant is kept:
335 -- X ** (abs Y) = Base**(Exp_High + Exp_Mid + Exp_Low) * Factor
339 Exp_High : Double := Double'Floor (Abs_Y);
342 Exp_Int : Mantissa_Type;
344 Factor : Double := 1.0;
347 -- Select algorithm for calculating Pow (integer cases fall through)
349 if Exp_High >= 2.0**Double'Machine_Mantissa then
351 -- In case of Y that is IEEE infinity, just raise constraint error
353 if Exp_High > Double'Safe_Last then
354 raise Constraint_Error;
357 -- Large values of Y are even integers and will stay integer
358 -- after division by two.
361 -- Exp_Mid and Exp_Low are zero, so
362 -- X**(abs Y) = Base ** Exp_High = (Base**2) ** (Exp_High / 2)
364 Exp_High := Exp_High / 2.0;
366 exit when Exp_High < 2.0**Double'Machine_Mantissa;
369 elsif Exp_High /= Abs_Y then
370 Exp_Low := Abs_Y - Exp_High;
373 if Exp_Low /= 0.0 then
375 -- Exp_Low now is in interval (0.0, 1.0)
376 -- Exp_Mid := Double'Floor (Exp_Low * 4.0) / 4.0;
379 Exp_Low := Exp_Low - Exp_Mid;
381 if Exp_Low >= 0.5 then
383 Exp_Low := Exp_Low - 0.5; -- exact
385 if Exp_Low >= 0.25 then
386 Factor := Factor * Sqrt (Factor);
387 Exp_Low := Exp_Low - 0.25; -- exact
390 elsif Exp_Low >= 0.25 then
391 Factor := Sqrt (Sqrt (X));
392 Exp_Low := Exp_Low - 0.25; -- exact
395 -- Exp_Low now is in interval (0.0, 0.25)
397 -- This means it is safe to call Logarithmic_Pow
398 -- for the remaining part.
400 Factor := Factor * Logarithmic_Pow (X, Exp_Low);
407 -- Exp_High is non-zero integer smaller than 2**Double'Machine_Mantissa
409 Exp_Int := Mantissa_Type (Exp_High);
411 -- Standard way for processing integer powers > 0
413 while Exp_Int > 1 loop
414 if (Exp_Int and 1) = 1 then
416 -- Base**Y = Base**(Exp_Int - 1) * Exp_Int for Exp_Int > 0
418 Factor := Factor * Base;
421 -- Exp_Int is even and Exp_Int > 0, so
422 -- Base**Y = (Base**2)**(Exp_Int / 2)
425 Exp_Int := Exp_Int / 2;
428 -- Exp_Int = 1 or Exp_Int = 0
431 Factor := Base * Factor;
435 Factor := 1.0 / Factor;
445 function Sin (X : Double) return Double is
446 Reduced_X : Double := X;
448 Quadrant : Natural range 0 .. 3;
451 if abs X > Pi / 4.0 then
452 Reduce (Reduced_X, Quadrant);
456 Asm (Template => "fsin",
457 Outputs => Double'Asm_Output ("=t", Result),
458 Inputs => Double'Asm_Input ("0", Reduced_X));
460 Asm (Template => "fcos",
461 Outputs => Double'Asm_Output ("=t", Result),
462 Inputs => Double'Asm_Input ("0", Reduced_X));
464 Asm (Template => "fsin",
465 Outputs => Double'Asm_Output ("=t", Result),
466 Inputs => Double'Asm_Input ("0", -Reduced_X));
468 Asm (Template => "fcos ; fchs",
469 Outputs => Double'Asm_Output ("=t", Result),
470 Inputs => Double'Asm_Input ("0", Reduced_X));
474 Asm (Template => "fsin",
475 Outputs => Double'Asm_Output ("=t", Result),
476 Inputs => Double'Asm_Input ("0", Reduced_X));
486 function Tan (X : Double) return Double is
487 Reduced_X : Double := X;
489 Quadrant : Natural range 0 .. 3;
492 if abs X > Pi / 4.0 then
493 Reduce (Reduced_X, Quadrant);
495 if Quadrant mod 2 = 0 then
496 Asm (Template => "fptan" & NL
497 & "ffree %%st(0)" & NL
499 Outputs => Double'Asm_Output ("=t", Result),
500 Inputs => Double'Asm_Input ("0", Reduced_X));
502 Asm (Template => "fsincos" & NL
503 & "fdivp %%st, %%st(1)" & NL
505 Outputs => Double'Asm_Output ("=t", Result),
506 Inputs => Double'Asm_Input ("0", Reduced_X));
512 & "ffree %%st(0) " & NL
514 Outputs => Double'Asm_Output ("=t", Result),
515 Inputs => Double'Asm_Input ("0", Reduced_X));
525 function Sinh (X : Double) return Double is
527 -- Mathematically Sinh (x) is defined to be (Exp (X) - Exp (-X)) / 2.0
530 return (Exp (X) - Exp (-X)) / 2.0;
532 return Exp (X) / 2.0;
540 function Cosh (X : Double) return Double is
542 -- Mathematically Cosh (X) is defined to be (Exp (X) + Exp (-X)) / 2.0
545 return (Exp (X) + Exp (-X)) / 2.0;
547 return Exp (X) / 2.0;
555 function Tanh (X : Double) return Double is
557 -- Return the Hyperbolic Tangent of x
561 -- Tanh (X) is defined to be ----------- = --------
566 return Double'Copy_Sign (1.0, X);
569 return 1.0 / (1.0 + Exp (-(2.0 * X))) - 1.0 / (1.0 + Exp (2.0 * X));
572 end Ada.Numerics.Aux;